From 83890dda858a0f6bdd969a79852d280a50d19106 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong @ vulgaris" Date: Fri, 11 Nov 2016 23:22:44 +0100 Subject: [PATCH] Initial commit. Updated material and funcs --- .gitignore | 14 +++ CMakeLists.txt | 139 +++++++++++++++++++++++++++ src/CMakeLists.txt | 69 +++++++++++++ src/funcs/CMakeLists.txt | 7 ++ src/funcs/bessel.cpp | 102 ++++++++++++++++++++ src/funcs/bessel.h | 23 +++++ src/funcs/cbessj.cpp | 187 ++++++++++++++++++++++++++++++++++++ src/funcs/cbessj.h | 22 +++++ src/funcs/rottfuncs.cpp | 126 ++++++++++++++++++++++++ src/funcs/rottfuncs.h | 43 +++++++++ src/funcs/skewsine.cpp | 27 ++++++ src/funcs/skewsine.h | 13 +++ src/funcs/sph_j.cpp | 32 ++++++ src/funcs/sph_j.h | 32 ++++++ src/material/CMakeLists.txt | 7 ++ src/material/air.cpp | 54 +++++++++++ src/material/air.h | 34 +++++++ src/material/gas.cpp | 34 +++++++ src/material/gas.h | 96 ++++++++++++++++++ src/material/helium.cpp | 33 +++++++ src/material/helium.h | 39 ++++++++ src/material/nitrogen.cpp | 49 ++++++++++ src/material/nitrogen.h | 36 +++++++ src/material/perfectgas.cpp | 7 ++ src/material/perfectgas.h | 78 +++++++++++++++ src/material/solid.cpp | 66 +++++++++++++ src/material/solid.h | 78 +++++++++++++++ src/tasmet_consolecolors.h | 27 ++++++ src/tasmet_enum.h | 61 ++++++++++++ src/tasmet_exception.cpp | 18 ++++ src/tasmet_exception.h | 26 +++++ src/tasmet_io.h | 25 +++++ src/tasmet_tracer.cpp | 13 +++ src/tasmet_tracer.h | 89 +++++++++++++++++ src/tasmet_types.h | 81 ++++++++++++++++ src/tasmet_utils.h | 71 ++++++++++++++ 36 files changed, 1858 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 src/CMakeLists.txt create mode 100644 src/funcs/CMakeLists.txt create mode 100644 src/funcs/bessel.cpp create mode 100644 src/funcs/bessel.h create mode 100644 src/funcs/cbessj.cpp create mode 100644 src/funcs/cbessj.h create mode 100644 src/funcs/rottfuncs.cpp create mode 100644 src/funcs/rottfuncs.h create mode 100644 src/funcs/skewsine.cpp create mode 100644 src/funcs/skewsine.h create mode 100644 src/funcs/sph_j.cpp create mode 100644 src/funcs/sph_j.h create mode 100644 src/material/CMakeLists.txt create mode 100644 src/material/air.cpp create mode 100644 src/material/air.h create mode 100644 src/material/gas.cpp create mode 100644 src/material/gas.h create mode 100644 src/material/helium.cpp create mode 100644 src/material/helium.h create mode 100644 src/material/nitrogen.cpp create mode 100644 src/material/nitrogen.h create mode 100644 src/material/perfectgas.cpp create mode 100644 src/material/perfectgas.h create mode 100644 src/material/solid.cpp create mode 100644 src/material/solid.h create mode 100644 src/tasmet_consolecolors.h create mode 100644 src/tasmet_enum.h create mode 100644 src/tasmet_exception.cpp create mode 100644 src/tasmet_exception.h create mode 100644 src/tasmet_io.h create mode 100644 src/tasmet_tracer.cpp create mode 100644 src/tasmet_tracer.h create mode 100644 src/tasmet_types.h create mode 100644 src/tasmet_utils.h diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0e90e0b --- /dev/null +++ b/.gitignore @@ -0,0 +1,14 @@ +CMakeCache.txt +CMakeFiles +Makefile +*.so +nonlinPYTHON_wrap.cxx +*.a +*.dll +*.pyd +__pycache__ +cmake_install.cmake +doc/html +doc/latex +doc/usg.pdf +*.pyc diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..2cdf38a --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,139 @@ +# CMakeList.txt for TaSMET +cmake_minimum_required (VERSION 2.8) +project(TaSMET) +set(PACKAGE_VERSION 0.1) +message("Running Cmake for TaSMET version ${PACKAGE_VERSION}") + +# Set the Python version to link against. If this is not set, Cmake tries to find it automatically. +# set(TaSMET_PY_VERSION "2.7") +set(TaSMET_PY_VERSION "3.5m") + +# Tracer name (name of the variable) +add_definitions(-DTRACER=1) +# add_definitions(-DTRACER_IN_COMMON) + +add_definitions(-DTASMET_FLOAT=64) + +#==================================================== +# Compiler settings ************************************************ +#==================================================== + +# Always required make flags in case of both compilers +set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -pipe -fPIC -Wall -Wextra -Wno-unused-parameter -Wno-unused-variable ") + + +# Stop letting Numpy complain about its API +add_definitions(-DNPY_NO_DEPRECATED_API=NPY_1_4_API_VERSION) + +#================================================== +# Optimized code flags ******************************************* +#================================================== + +#Debug mode +# set(CMAKE_CXX_FLAGS "${CMAKE_CXX)FLAGS} -g -ggdb") + +# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2") +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -march=native -mtune=native") +# set(CMAKE_CLANG "${CMAKE_GCC} -march=native -mtune=native -fopenmp") + +# To disable bound checking on std::vector, and to disable assertions +# add_definitions(-DNDEBUG) + +# To increase speed on Armadillo +# add_definitions(-DARMA_NO_DEBUG) + +# Disable traces +# add_definitions(-DTRACER=0) + +# Pre-allocation size for matrices. Very important setting to tune the +# code in case we know that we are going to run with certain sizes of matrices and vectors. For Nf=6, set this to (2*Nf+1)^2=169 +# For Nf=1: 9 +# For Nf=2: 25 +# For Nf=3: 49 +# For Nf=4: 81 +# For Nf=5: 121 +# For Nf=6: 169 +# For Nf=7: 225 +# For Nf=8: 289 +# For Nf=9: 361 +# For Nf=10: 441 +# For Nf=11: 529 +# For Nf=12: 625 +# Watch out! Setting prealloc too high can give too much overhead for smaller Nf's +# add_definitions(-DARMA_MAT_PREALLOC=625) + + +# ########################## +# Finding the presence of the prerequisites +# ########################## + + +# For importing find directives for Cmake +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/common/cmake_tools) + +# ########################## +# Python ##################### +# ########################## + +if(TaSMET_PY_VERSION) + # Find major version from version string + set(PYTHON_LIBRARY "/usr/lib/libpython${TaSMET_PY_VERSION}.so") + set(PYTHON_INCLUDE_DIR "/usr/include/python${TaSMET_PY_VERSION}") + set(PYTHON_INCLUDE_DIRS "/usr/include/python${TaSMET_PY_VERSION}") + +endif(TaSMET_PY_VERSION) +message("Python include dirs: ${PYTHON_INCLUDE_DIRS}") + +find_package(PythonLibs REQUIRED) +string(REGEX MATCH "^." TaSMET_PY_MAJOR_VERSION ${PYTHONLIBS_VERSION_STRING}) +MESSAGE("Python major version: ${TaSMET_PY_MAJOR_VERSION}") + +# Find the site_packages directory of python +execute_process(COMMAND python${TaSMET_PY_MAJOR_VERSION} -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())" OUTPUT_VARIABLE PYTHON_SITE_PACKAGES OUTPUT_STRIP_TRAILING_WHITESPACE) + +# ################################ +# Initialize swig +# ################################ + +find_package(SWIG REQUIRED) +include(${SWIG_USE_FILE}) + +SET(CMAKE_SWIG_FLAGS -Wall -DSWIG_PYTHON) +if(${TaSMET_PY_MAJOR_VERSION}=="3") + SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} -py${TaSMET_PY_MAJOR_VERSION}) +endif(${TaSMET_PY_MAJOR_VERSION}=="3") + +include_directories(common/src) + +# Armadillo +find_package(Armadillo REQUIRED) +add_definitions(-DARMA_USE_SUPERLU -DARMA_USE_CXX11) + +# ==================== Compile the code in common and src + +# This is the common code (gas and solid libs, etc) +# link_directories(common) + +aux_source_directory(common/src/gas gas) +include_directories( + ${PYTHON_INCLUDE_DIRS} + common/src/swig + common/src/gas + ) +# Add the code subdirectory +add_subdirectory(src) + + +# set_source_files_properties( ${swig_generated_file_fullname} + # PROPERTIES COMPILE_FLAGS "${SWIG_COMMON_COMPILE_FLAGS} ") + + +# ================================== Installation + +# Install common files +install(FILES ${PROJECT_SOURCE_DIR}/common/__init__.py + DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}/common) +install(FILES ${PROJECT_SOURCE_DIR}/common/common.py + DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}/common) + +# Rest of the files is installed from src/CMakeLists.txt diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..267f32f --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,69 @@ +#================ +# The TaSMET code +#================ + +include_directories( + ${PYTHON_INCLUDE_PATH} + ${ARMADILLO_INCLUDE_DIRS} + . + # duct + # duct/cell + # duct/connectors + # duct/eq + # duct/geom + # duct/geom/grid + # duct/drag + # duct/heat + # mech + # python + # seg + # sol + # sys + # var + # volume + ) + +add_subdirectory(funcs) +add_subdirectory(material) +# add_subdirectory(duct) +# add_subdirectory(mech) +# add_subdirectory(seg) +# add_subdirectory(sol) +# add_subdirectory(sys) +# add_subdirectory(var) + +add_library(tasmet_src tasmet_tracer.cpp tasmet_exception.cpp) +target_link_libraries(tasmet_src + # duct + # mech + # seg + # sol + # sys + # var + ) + +# set_source_files_properties(swig/nonlin.i +# PROPERTIES CPLUSPLUS ON) +# swig_add_module(TaSMET python swig/nonlin.i) + +# set(SWIG_MODULE_TaSMET_FILE_nonlin_EXTRA_DEPS +# ${CMAKE_CURRENT_SOURCE_DIR}/common/swig_files/arma_numpy.i +# ) +# swig_link_libraries(TaSMET tasmet_src common ${PYTHON_LIBRARIES} +# ${ARMADILLO_LIBRARIES}) + +# set_source_files_properties( ${swig_generated_file_fullname} + # PROPERTIES COMPILE_FLAGS "${SWIG_COMMON_COMPILE_FLAGS} ") + +# Install swig files to the right place +# install(TARGETS _TaSMET +# DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}) +# install(FILES ${CMAKE_BINARY_DIR}/src/TaSMET.py +# DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}) +# install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/python/__init__.py +# DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}) +# install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/python/post +# DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}) +# install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/python/gui +# DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}) + diff --git a/src/funcs/CMakeLists.txt b/src/funcs/CMakeLists.txt new file mode 100644 index 0000000..1c75ddb --- /dev/null +++ b/src/funcs/CMakeLists.txt @@ -0,0 +1,7 @@ +add_library(funcs + bessel.cpp + cbessj.cpp + rottfuncs.cpp + skewsine.cpp + sph_j.cpp + ) diff --git a/src/funcs/bessel.cpp b/src/funcs/bessel.cpp new file mode 100644 index 0000000..33a04db --- /dev/null +++ b/src/funcs/bessel.cpp @@ -0,0 +1,102 @@ +#include "bessel.h" +#include "cbessj.h" +#include "skewsine.h" + +namespace { + c besselj0_smallarg(c x){ + //Abramowich p. 369, polynomial approximation + //For complex numbers, unfortunately we have to check + //for both imaginary and real part. + /* return 1.0 \ + // -2.2499997*pow(x/3.0,2) \ + // +1.2656208*pow(x/3.0,4) \ + // -0.3163866*pow(x/3.0,6) \ + // +0.0444479*pow(x/3.0,8) \ + // -0.0039444*pow(x/3.0,10) \ + // +0.0002100*pow(x/3.0,12); */ + + // External implementation + d z[2]; + d res[2]; + z[0]=x.real(); + z[1]=x.imag(); + CBESSJ(z,0,res); + return c(res[0],res[1]); + } + c besselj0_largearg(c x){ + // larger than 3" << endl; + // Abramovich, p. 370 + c f0=0.79788456-0.00000077*(3.0/x)- + 0.00552740*pow(3.0/x,2)- + 0.00009512*pow(3.0/x,3)+ + 0.00137237*pow(3.0/x,4)- + 0.00072805*pow(3.0/x,5)+ + 0.00014476*pow(3.0/x,6); + c th0=x-0.78539816- + 0.04166397*(3.0/x)- + 0.00003954*pow(3.0/x,2)+ + 0.00262573*pow(3.0/x,3)- + 0.00054125*pow(3.0/x,4)- + 0.00029333*pow(3.0/x,5)+ + 0.00013558*pow(3.0/x,6); + return pow(x,-0.5)*f0*cos(th0); + } + c besselj1_over_x_smallarg(c x){ + // return 0.5-0.56249985*pow(x/3.0,2)+0.21093573*pow(x/3.0,4)-0.03954289*pow(x/3.0,6)+0.00443319*pow(x/3.0,8)-0.00031761*pow(x/3.0,10)+0.00001109*pow(x/3.0,12); + + // External implementation + d z[2]; + d res[2]; + z[0]=x.real(); + z[1]=x.imag(); + CBESSJ(z,1,res); + return c(res[0],res[1])/x; + } + c besselj1_over_x_largearg(c x){ + c f1= + 0.79788456 + +0.00000156*(3.0/x) + +0.01659667*pow(3.0/x,2) + +0.00017105*pow(3.0/x,3) + -0.00249511*pow(3.0/x,4) + +0.00113653*pow(3.0/x,5) + -0.00020033*pow(3.0/x,6); + + c th1=x + -2.35619449 + +0.12499612*(3.0/x) + +0.00005650*pow(3.0/x,2) + -0.00637879*pow(3.0/x,3) + +0.00074348*pow(3.0/x,4) + +0.00079824*pow(3.0/x,5) + -0.00029166*pow(3.0/x,6); + return pow(x,-1.5)*f1*cos(th1); + + } + +} + +c besselj0(c x){ + d trstart=11.0; // Transition start point + d trdelta=1.0; // Transition deltax + if (abs(x) +#include + +#define MAXK 20 //09/21/2009 (15 before) +#define FALSE 0 +#define TRUE 1 + +double HALF = 0.5, ONE = 1.0, FPF = 5.5; +double PI=4.0*atan(1.0); + +int nu; // order of complex Bessel +double z[2], z1[2]; // Complex numbers +double x,y; + +// Z=Z1/Z2 +void CDIV(double *Z1, double *Z2, double *Z) { + double D; + D=Z2[0]*Z2[0]+Z2[1]*Z2[1]; + if (D > 1e30) { + Z[0]=0.0; Z[1]=0.0; return; + } + if (D == 0) return; + Z[0]=(Z1[0]*Z2[0]+Z1[1]*Z2[1])/D; + Z[1]=(Z1[1]*Z2[0]-Z1[0]*Z2[1])/D; +} + +// Z=Z1*Z2 +void CMUL(double *Z1, double *Z2, double *Z) { + Z[0]=Z1[0]*Z2[0] - Z1[1]*Z2[1]; + Z[1]=Z1[0]*Z2[1] + Z1[1]*Z2[0]; +} + +// compute Z^N +void IZPower(double *z, int n, double *z1) { + double temp[2],temp1[2]; int i; + if (n==0) { + z1[0]=1.0; + z1[1]=0.0; + } + else if (n==1) { + z1[0]=z[0]; + z1[1]=z[1]; + } + else { + temp1[0]=z[0]; temp1[1]=z[1]; + for (i=2; i<=n; i++) { + CMUL(temp1,z,temp); + temp1[0]=temp[0]; + temp1[1]=temp[1]; + } + z1[0]=temp[0]; + z1[1]=temp[1]; + } +} + +double Fact(int k) { + int i; double f; + f=ONE; + for (i=2; i<=k; i++) f *= ONE*i; + return f; +} + +/****************************************** + * FUNCTION GAMMA(X) * + * --------------------------------------- * + * Returns the value of Gamma(x) in double * + * precision as EXP(LN(GAMMA(X))) for X>0. * + ******************************************/ +double Gamma(double xx) { + double cof[6]; + double stp,x,tmp,ser; + int j; + cof[0]=76.18009173; + cof[1]=-86.50532033; + cof[2]=24.01409822; + cof[3]=-1.231739516; + cof[4]=0.120858003e-2; + cof[5]=-0.536382e-5; + stp=2.50662827465; + x=xx-ONE; + tmp=x+FPF; + tmp=(x+HALF)*log(tmp)-tmp; + ser=ONE; + for (j=0; j<6; j++) { + x += ONE; + ser += cof[j]/x; + } + return (exp(tmp+log(stp*ser))); +} + +// main subroutine +void CBESSJ(double *z, int nu, double *z1) { + /*-------------------------------------------------- + inf. (-z^2/4)^k + Jnu(z) = (z/2)^nu x Sum ------------------ + k=0 k! x Gamma(nu+k+1) + (nu must be >= 0). Here k=15. + ---------------------------------------------------*/ + int k; + double sum[2],tmp[2],tmp1[2],tmp2[2],tmp3[2]; + // calculate (z/2)^nu in tmp3 + tmp[0]=2.0; tmp[1]=0.0; + CDIV(z,tmp,tmp1); + IZPower(tmp1,nu,tmp3); + sum[0]=0.0; sum[1]=0.0; + //calculate Sum + for (k=0; k<=MAXK; k++) { + // calculate (-z^2/4)^k + IZPower(z,2,tmp); + tmp[0]=-tmp[0]; tmp[1]=-tmp[1]; + tmp1[0]=4.0; tmp1[1]=0.0; + CDIV(tmp,tmp1,tmp2); + + if (z[1]==0) { //case real number + tmp[0]=pow(tmp2[0],k); + tmp[1]=0.0; + } + else //case complex number + IZPower(tmp2,k,tmp); + + // divide by k! + tmp1[0]=Fact(k); tmp1[1]=0.0; + if (z[1]==0) { + tmp2[0]=tmp[0]/tmp1[0]; + tmp2[1]=0.0; + } + else + CDIV(tmp,tmp1,tmp2); + // divide by Gamma(nu+k+1) + tmp1[0]=Gamma(ONE*(nu+k+1)); tmp1[1]=0.0; + + if (z[1]==0) { + tmp[0]=tmp2[0]/tmp1[0]; + tmp[1]=0.0; + } + else + CDIV(tmp2,tmp1,tmp); + // actualize sum + sum[0] += tmp[0]; + sum[1] += tmp[1]; + + } + // multiply (z/2)^nu by sum + CMUL(tmp3,sum,z1); +} + + +// int main() { + +// double z1[2]; + +// printf("\n Complex Bessel Function of the 1st Kind of integer order\n\n"); +// printf(" Input complex argument (real imaginary): "); +// scanf("%lf %lf", &x, &y); +// z[0]=x; z[1]=y; +// printf(" Input integer order: "); scanf("%d", &nu); + +// CBESSJ(z,nu,z1); + +// printf("\n Function value: %f %f\n\n", z1[0], z1[1]); + +// } + +// end of file cbessj.cpp diff --git a/src/funcs/cbessj.h b/src/funcs/cbessj.h new file mode 100644 index 0000000..67ed0a7 --- /dev/null +++ b/src/funcs/cbessj.h @@ -0,0 +1,22 @@ +// cbessj.h +// +// Author: J.A. de Jong +// +// Description: +// Implementation of the besselJ function for complex numbers +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef CBESSJ_H +#define CBESSJ_H + +// z: input, z1: result; + +void CBESSJ(double *z, int nu, double *z1); + + +#endif // CBESSJ_H +////////////////////////////////////////////////////////////////////// + + + + diff --git a/src/funcs/rottfuncs.cpp b/src/funcs/rottfuncs.cpp new file mode 100644 index 0000000..ebd8697 --- /dev/null +++ b/src/funcs/rottfuncs.cpp @@ -0,0 +1,126 @@ +// rottfuncs.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// Wat doet dit bestand? +////////////////////////////////////////////////////////////////////// +#include "rottfuncs.h" +#include "bessel.h" + +namespace { + template + array_type f_inviscid(const array_type& rh_over_delta){ + TRACE(0,"f_inviscid"); + return 0.0*rh_over_delta; + } + + template + array_type f_vert(const array_type& rh_over_delta){ //Vertical plate f_nu/f_k + TRACE(0,"f_vert"); + c a=I+1.0; + return tanh(a*rh_over_delta)/(a*rh_over_delta); + } + + template + array_type f_blapprox(const array_type& rh_over_delta){ + TRACE(0,"f_blapprox()"); + return (1.0-I)/rh_over_delta/2.0; + } + + c f_circ(const c& rh_over_delta) { + TRACE(0,"f_circ(c rh_ov_delta)"); + // vc s=sq2*rh_over_delta; + c a(-1,1.0); + c R_over_delta=2.0*rh_over_delta; + c jvarg=a*R_over_delta; + c j0=besselj0(jvarg); + c j1_over_jxarg=besselj1_over_x(jvarg);; + return 2.0*j1_over_jxarg/j0; + } + vc f_circ(const vc& rh_over_delta) { + TRACE(0,"f_circ()"); + return vc(rh_over_delta).transform([] (c x) {return f_circ(x);}); + } + + c f_square(const c& rh_over_delta){ + TRACE(0,"f_square()"); + // Adjusted sligthly for improved speed. + // Needs checking if rh is defined right. + // For square pore with dimensions AxA : rh=S/PI = A^2 / 4A = A/4 + c F(0,0); + c C(0,0); + us n,m,msq,nsq; + d pisq=pow(number_pi,2); + for(n=1; n<11; n=n+2) + { + for(m = 1; m<11; m=m+2) + { + msq=m*m; + nsq=n*n; + C = 1.0-I*pisq/(32.0*pow(rh_over_delta,2))*((double)(msq+nsq)); + F+=(1.0/msq/nsq)/C; + } + } + return 1.0-64.0/(pisq*pisq)*F; + } + // To transform this one would introduce an unallowable loop overhead + vc f_square(const vc& rh_over_delta){ + TRACE(0,"f_square()"); + // Adjusted sligthly for improved speed. + // Needs checking if rh is defined right. + // For square pore with dimensions AxA : rh=S/PI = A^2 / 4A = A/4 + vc F(rh_over_delta.size(),fillwith::zeros); + vc C(rh_over_delta.size(),fillwith::zeros); + us n,m,msq,nsq; + d pisq=pow(number_pi,2); + for(n=1; n<11; n=n+2) + { + for(m = 1; m<11; m=m+2) + { + msq=m*m; + nsq=n*n; + C = 1.0-I*pisq/(32.0*pow(rh_over_delta,2))*(msq+nsq); + F+=(1.0/msq/nsq)/C; + } + } + return 1.0-64.0/(pisq*pisq)*F; + } +} +RottFuncs::RottFuncs(const cshape shape):_shape(shape){ + const string shapestr = cshapeToString(shape); + TRACE(8,"RottFuncs::RottFuncs(" << shapestr << ")"); + + switch (shape) { + case CIRC: + TRACE(10,"Fptr set to circ") + f_ptr=&f_circ; + f_ptrc=&f_circ; + break; + case VERT: + TRACE(10,"Fptr set to vert"); + f_ptr=&f_vert; + f_ptrc=&f_vert; + break; + case SQUARE: + TRACE(10,"Fptr set to square"); + f_ptr=&f_square; + f_ptrc=&f_square; + break; + case BLAPPROX: + TRACE(10,"Fptr set to blapprox"); + f_ptr=&f_blapprox; + f_ptrc=&f_blapprox; + break; + case INVISCID: + TRACE(10,"Fptr set to inviscid"); + f_ptr=&f_inviscid; + f_ptrc=&f_inviscid; + break; + default: + FATAL("Unhandled shape"); + break; + } +} + +////////////////////////////////////////////////////////////////////// diff --git a/src/funcs/rottfuncs.h b/src/funcs/rottfuncs.h new file mode 100644 index 0000000..78fd1b6 --- /dev/null +++ b/src/funcs/rottfuncs.h @@ -0,0 +1,43 @@ +// rottfuncs.h +// +// Author: J.A. de Jong +// +// Description: +// Implementation of thermoviscous (Rott's) functions +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef ROTTFUNCS_H +#define ROTTFUNCS_H +#include "tasmet_types.h" +#include "tasmet_tracer.h" +#include "tasmet_enum.h" + +DECLARE_ENUM(cshape,INVISCID,BLAPPROX,VERT,SQUARE,CIRC) + +class RottFuncs { // Thermoviscous Rott functions + vc (*f_ptr)(const vc& rh_over_delta); + c (*f_ptrc)(const c& rh_over_delta); + cshape _shape; +public: + RottFuncs(const cshape shape=cshape::INVISCID); + RottFuncs(const RottFuncs& other) =delete; + RottFuncs& operator=(const RottFuncs&) =delete; + ~RottFuncs(){} + + vc operator()(const vc& rh_over_delta) const { + TRACE(5,"fx vc"); + return f_ptr(rh_over_delta); + } + vc operator()(const vd& rh_over_delta) const { + TRACE(5,"fx vd"); + return f_ptr((1.0+0.0*I)*rh_over_delta); + } + + c operator()(const c& rh_over_delta) const {return f_ptrc(rh_over_delta);} + c operator()(const d& rh_over_delta) const {c n(rh_over_delta,0); return f_ptrc(n);} + cshape getCshape() const {return _shape;} +}; + + +#endif // ROTTFUNCS_H +////////////////////////////////////////////////////////////////////// diff --git a/src/funcs/skewsine.cpp b/src/funcs/skewsine.cpp new file mode 100644 index 0000000..ca934ba --- /dev/null +++ b/src/funcs/skewsine.cpp @@ -0,0 +1,27 @@ +// skewsine.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// Implementation of a skew sine +////////////////////////////////////////////////////////////////////// + +#include "skewsine.h" + +// Equation for a skew sine between x=0 and x=1 +d skewsine(d x){ + if (x<0.0) + return 0.0; + else if (x<1.0) + return x-(1.0/(2.0*number_pi))*sin(x*number_pi*2.0); + else + return 1.0; +} + +vd skewsine(const vd& x){ + vd result=x; + result.transform( [](d val) {return skewsine(val);}); + return result; +} + +////////////////////////////////////////////////////////////////////// diff --git a/src/funcs/skewsine.h b/src/funcs/skewsine.h new file mode 100644 index 0000000..97d806b --- /dev/null +++ b/src/funcs/skewsine.h @@ -0,0 +1,13 @@ +#pragma once + +#ifndef _SKEWSINE_H_ +#define _SKEWSINE_H_ +#include "tasmet_types.h" + + +d skewsine(d x); // Skew sine between x=0 and x=1 +vd skewsine(const vd& x); // Skew sine between x=0 and x=1 + + + +#endif /* _SKEWSINE_H_ */ diff --git a/src/funcs/sph_j.cpp b/src/funcs/sph_j.cpp new file mode 100644 index 0000000..ef1fd6d --- /dev/null +++ b/src/funcs/sph_j.cpp @@ -0,0 +1,32 @@ +// sph_j.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#include +#include "sph_j.h" +#include "tasmet_exception.h" + + +vd sph_j(int n,d x){ + if (n < 1) + throw TaSMETError("n needs to be > 1"); + vd result(n); + gsl_sf_bessel_jl_steed_array(n,x,result.memptr()); + return result; +} + +dmat sph_j(int n, const vd& x){ + if (n < 1) + throw TaSMETError("n needs to be > 1"); + + dmat result(x.size(),n); + for (us i = 0; i < x.size(); i++) + result.row(i) = sph_j(n, x(i)).t(); + return result; +} + + +////////////////////////////////////////////////////////////////////// diff --git a/src/funcs/sph_j.h b/src/funcs/sph_j.h new file mode 100644 index 0000000..95096d7 --- /dev/null +++ b/src/funcs/sph_j.h @@ -0,0 +1,32 @@ +// sph_j.h +// +// Author: J.A. de Jong +// +// Description: +// Implementation of spherical bessel functions +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef SPH_J_H +#define SPH_J_H +#include "tasmet_types.h" + + +// Returns the first n spherical Bessel functions for a given double value +vd sph_j(int n,d x); +// Returns a 2D array. The first axis is the value of x and the +// second axis (say the column number) corresponds to the n'th +// spherical Bessel function value. Hence the result vector looks like: + +// [ j_0(x_0) j_1(x_0) ... j_n(x_0) ] +// [ j_0(x_1) j_1(x_1_ ... j_n(x_1) ] +// . . +// . . +// . . +// [ j_0(x_n) j_1(x_n) ... j_n(x_n) ] + +dmat sph_j(int n, const vd& x); + + + +#endif // SPH_J_H +////////////////////////////////////////////////////////////////////// diff --git a/src/material/CMakeLists.txt b/src/material/CMakeLists.txt new file mode 100644 index 0000000..96e6462 --- /dev/null +++ b/src/material/CMakeLists.txt @@ -0,0 +1,7 @@ +add_library(material + gas.cpp + air.cpp + helium.cpp + nitrogen.cpp + solid.cpp + ) diff --git a/src/material/air.cpp b/src/material/air.cpp new file mode 100644 index 0000000..5a8f321 --- /dev/null +++ b/src/material/air.cpp @@ -0,0 +1,54 @@ +// air.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// + +#include "air.h" + +namespace { + // Air-specific data + d kappac[]={-0.00227583562,1.15480022E-4, \ + -7.90252856E-8,4.11702505E-11, \ + -7.43864331E-15}; + + d cpc[]={1047.63657,-0.372589265, \ + 9.45304214E-4,-6.02409443E-7, \ + 1.2858961E-10}; + d muc[]={-8.38278E-7,8.35717342E-8, \ + -7.69429583E-11,4.6437266E-14, \ + -1.06585607E-17}; +} // namespace + +d Air::h(d T,d p) const { + return cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5); +} +vd Air::h(const vd& T,const vd& p) const { + return cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5); +} +vd Air::cp(const vd& T,const vd& p) const { + return cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4); +} +d Air::cp(d T,d p) const { + return cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4); +} + +vd Air::kappa(const vd& T,const vd& p) const { + return kappac[0]+kappac[1]*T+kappac[2]*pow(T,2)+kappac[3]*pow(T,3)+kappac[4]*pow(T,4); +} +d Air::kappa(d T,d p) const { + return kappac[0]+kappac[1]*T+kappac[2]*pow(T,2)+kappac[3]*pow(T,3)+kappac[4]*pow(T,4); +} +vd Air::mu(const vd& T,const vd& p) const { + return muc[1]*T+muc[2]*pow(T,2)+muc[3]*pow(T,3)+muc[4]*pow(T,4)+muc[0]; +} +d Air::mu(d T,d p) const { + return muc[1]*T+muc[2]*pow(T,2)+muc[3]*pow(T,3)+muc[4]*pow(T,4)+muc[0]; +} + + + + +////////////////////////////////////////////////////////////////////// diff --git a/src/material/air.h b/src/material/air.h new file mode 100644 index 0000000..1fa221d --- /dev/null +++ b/src/material/air.h @@ -0,0 +1,34 @@ +// air.h +// +// Author: J.A. de Jong +// +// Description: +// Implementation of the gas air +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef AIR_H +#define AIR_H +#include "perfectgas.h" + + +class Air : public PerfectGas { +protected: + d Rs() const {return 287;} +public: + Air():PerfectGas(air){} + d cp(d T,d p) const; + vd cp(const vd& T,const vd& p) const; + d h(d T,d p) const; + vd h(const vd& T,const vd& p) const; + d mu(d T,d p) const; + vd mu(const vd& T,const vd& p) const; + d kappa(d T,d p) const; + vd kappa(const vd& T,const vd& p) const; + ~Air(){} +}; + + + + +#endif // AIR_H +////////////////////////////////////////////////////////////////////// diff --git a/src/material/gas.cpp b/src/material/gas.cpp new file mode 100644 index 0000000..f398d8a --- /dev/null +++ b/src/material/gas.cpp @@ -0,0 +1,34 @@ +// gas.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#include "tasmet_tracer.h" +#include "tasmet_exception.h" +#include "gas.h" +#include "helium.h" +#include "air.h" +#include "nitrogen.h" + +Gas* Gas::newGas(const GasType gastype) { + + switch (gastype) { + case helium: + return new Helium(); + break; + case air: + return new Air(); + case nitrogen: + return new Nitrogen(); + default: + FATAL("Gas type not known"); + break; + } + return nullptr; +} + + +////////////////////////////////////////////////////////////////////// + diff --git a/src/material/gas.h b/src/material/gas.h new file mode 100644 index 0000000..3a401bc --- /dev/null +++ b/src/material/gas.h @@ -0,0 +1,96 @@ +// gas.h +// +// Author: J.A. de Jong +// +// Description: +// Gas properties interface +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef GAS_H +#define GAS_H +#include +#include +#include "tasmet_enum.h" +#include "tasmet_types.h" + + + +class Gas{ +public: + DECLARE_ENUM(GasType,air,helium,nitrogen); +private: + GasType _gastype; +protected: + Gas(GasType gastype):_gastype(gastype){} +public: + + + + Gas(const Gas& other) =delete; + Gas& operator=(const Gas&)=delete; + virtual ~Gas(){} + + // Static method to generate a Gas + static Gas* newGas(const GasType gastype); + operator GasType() { return _gastype;} + + // Dimensionless numbers: + + // Specific heat ratio + d gamma(d T,d p) const {return cp(T,p)/cv(T,p);} + vd gamma(const vd& T,const vd& p) const {return cp(T,p)/cv(T,p);} + + // Prandtl number + vd pr(const vd& T,const vd& p) const {return mu(T,p)%cp(T,p)/kappa(T,p);} + d pr(d T,d p) const {return mu(T,p)*cp(T,p)/kappa(T,p);} + + // Virtuals that are part of the interface + // Density [kg/m^3] + virtual d rho(d T,d p) const=0; + virtual vd rho(const vd& T,const vd& p) const=0; + + // Internal energy [J/kg] + virtual d e(d T,d p) const=0; + virtual vd e(const vd& T,const vd& p) const=0; + + // Static enthalpy [J/kg] + virtual d h(d T,d p) const=0; + virtual vd h(const vd& T,const vd& p) const=0; + + // Specific heat at constant pressure + virtual d cp(d T,d p) const=0; + virtual vd cp(const vd& T,const vd& p) const=0; + + // Specific heat at constant density + virtual d cv(d T,d p) const=0; + virtual vd cv(const vd& T,const vd& p) const=0; + + virtual d beta(d T,d p) const=0; + virtual vd beta(const vd& T,const vd& p) const=0; + + // Adiabatic speed of sound [m/s] + virtual d cm(d T,d p) const=0; + virtual vd cm(const vd& T,const vd& p) const=0; + + // Dynamic viscosity [Pa*s] + virtual d mu(d T,d p) const=0; + virtual vd mu(const vd& T,const vd& p) const=0; + + // Thermal conductivity [W/mK] + virtual d kappa(d T,d p) const=0; + virtual vd kappa(const vd& T,const vd& p) const=0; + + + + + + + + + +}; + + + +#endif // GAS_H +////////////////////////////////////////////////////////////////////// diff --git a/src/material/helium.cpp b/src/material/helium.cpp new file mode 100644 index 0000000..19c84d1 --- /dev/null +++ b/src/material/helium.cpp @@ -0,0 +1,33 @@ +// helium.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// + +#include "helium.h" + + +vd Helium::cp(const vd& T,const vd& p) const { + return cp(0.0,0.0)*arma::ones(T.size()); +} +vd Helium::h(const vd& T,const vd& p) const { + return cp(T,p)%T; +} +d Helium::mu(d T,d p) const { + return 0.412e-6*pow(T,0.68014); +} +d Helium::kappa(d T,d p) const { + return 0.0025672*pow(T,0.716); +} +vd Helium::mu(const vd& T,const vd& p) const { + return 0.412e-6*pow(T,0.68014); +} +vd Helium::kappa(const vd& T,const vd& p) const { + return 0.0025672*pow(T,0.716); +} + + + +////////////////////////////////////////////////////////////////////// diff --git a/src/material/helium.h b/src/material/helium.h new file mode 100644 index 0000000..a422dc5 --- /dev/null +++ b/src/material/helium.h @@ -0,0 +1,39 @@ +// helium.h +// +// Author: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef HELIUM_H +#define HELIUM_H +#include "perfectgas.h" + +class Helium :public PerfectGas { +protected: + d Rs() const {return 2070;} +public: + Helium():PerfectGas(helium){} + + d cp(d T,d p) const { return 5195;} + vd cp(const vd& T,const vd& p) const; + + d h(d T,d p) const {return cp(0.0,0.0)*T;} + vd h(const vd& T,const vd& p) const; + + d mu(d T,d p) const; + vd mu(const vd& T,const vd& p) const; + + d kappa(d T,d p) const; + vd kappa(const vd& T,const vd& p) const; + + virtual ~Helium(){} +}; + + + + + +#endif // HELIUM_H +////////////////////////////////////////////////////////////////////// diff --git a/src/material/nitrogen.cpp b/src/material/nitrogen.cpp new file mode 100644 index 0000000..8ee8bb8 --- /dev/null +++ b/src/material/nitrogen.cpp @@ -0,0 +1,49 @@ +// nitrogen.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// + +#include "nitrogen.h" + +namespace { + // Nitrogen-specific data + d cpc[]={3.29868,0.00140824, \ + -3.96322e-6,5.64152e-9, \ + -2.44486e-12}; + + d kappavals[]={6.995e-3,0.0631e-3}; +} // namespace + +d Nitrogen::h(d T,d p) const { + return Rs()*(cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5)); +} +vd Nitrogen::h(const vd& T,const vd& p) const { + return Rs()*(cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5)); +} +vd Nitrogen::cp(const vd& T,const vd& p) const { + return Rs()*(cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4)); +} +d Nitrogen::cp(d T,d p) const { + return Rs()*(cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4)); +} +// Document Mina +vd Nitrogen::kappa(const vd& T,const vd& p) const { + return kappavals[1]*T+kappavals[0]; +} +d Nitrogen::kappa(d T,d p) const { + return kappavals[1]*T+kappavals[0]; +} +// http://www.lmnoeng.com/Flow/GasViscosity.php +vd Nitrogen::mu(const vd& T,const vd& p) const { + return (0.01781/1000)*(411.55/(T+111))%pow(T/300.55,1.5); +} +d Nitrogen::mu(d T,d p) const { + return (0.01781/1000)*(411.55/(T+111))*pow(T/300.55,1.5); +} + + + +////////////////////////////////////////////////////////////////////// diff --git a/src/material/nitrogen.h b/src/material/nitrogen.h new file mode 100644 index 0000000..c0fefb3 --- /dev/null +++ b/src/material/nitrogen.h @@ -0,0 +1,36 @@ +// nitrogen.h +// +// Author: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef NITROGEN_H +#define NITROGEN_H + +#include "perfectgas.h" + + +class Nitrogen : public PerfectGas { +protected: + d Rs() const {return 297;} +public: + Nitrogen():PerfectGas(nitrogen){} + + d cp(d T,d p) const; + vd cp(const vd& T,const vd& p) const; + d h(d T,d p) const; + vd h(const vd& T,const vd& p) const; + d mu(d T,d p) const; + vd mu(const vd& T,const vd& p) const; + d kappa(d T,d p) const; + vd kappa(const vd& T,const vd& p) const; + ~Nitrogen(){} +}; + + + +#endif // NITROGEN_H + +////////////////////////////////////////////////////////////////////// diff --git a/src/material/perfectgas.cpp b/src/material/perfectgas.cpp new file mode 100644 index 0000000..9483462 --- /dev/null +++ b/src/material/perfectgas.cpp @@ -0,0 +1,7 @@ +// perfectgas.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// diff --git a/src/material/perfectgas.h b/src/material/perfectgas.h new file mode 100644 index 0000000..5ec3258 --- /dev/null +++ b/src/material/perfectgas.h @@ -0,0 +1,78 @@ +// perfectgas.h +// +// Author: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef PERFECTGAS_H +#define PERFECTGAS_H +#include "gas.h" +#ifdef DEBUG_GAS +#define checkzero(x) \ + TRACE(-1,"Divide by zero testing entered."); \ + try {if(min(abs(x))<1e-13) throw 0; } \ + catch(int a){ \ + TRACE(0,"Divide by zero encountered."); \ + } +#else +#define checkzero(x) +#endif + + +class PerfectGas: public Gas { +protected: + PerfectGas(GasType gastype): Gas(gastype){} + + // Not implemented for a perfect gas: + // mu, kappa, h, cp, Rs + + virtual d Rs() const=0; + +public: + ~PerfectGas(){} + + vd rho(const vd&T,const vd&p) const { + checkzero(T); + return p/Rs()/T; + } + vd p(const vd& T,const vd& rho) const { + return rho%(Rs()*T); + } + vd cv(const vd& T,const vd& p) const { + return cp(T,p)-Rs(); + } + vd e(const vd& T,const vd& p) const { + return h(T,p)-Rs()*T; + } + d beta(d T,d p) const { + checkzero(T); + return 1/T; + } + vd beta(const vd& T,const vd& p) const { + return 1/T; + } + vd cm(const vd& T,const vd& p) const { + return sqrt(gamma(T,p)*Rs()%T); + } + d rho(d T,d p) const { + checkzero(T); + return p/Rs()/T; + } + d cv(d T,d p) const { + return cp(T,p)-Rs(); + } + d e(d T,d p) const { + return h(T,p)-Rs()*T; + } + d cm(d T,d p) const { + d csq=gamma(T,p)*Rs()*T; + return sqrt(csq); + } + +}; + +#endif // PERFECTGAS_H +////////////////////////////////////////////////////////////////////// + diff --git a/src/material/solid.cpp b/src/material/solid.cpp new file mode 100644 index 0000000..9eb2268 --- /dev/null +++ b/src/material/solid.cpp @@ -0,0 +1,66 @@ +#include "solid.h" +#include "tasmet_tracer.h" + +//Stainless steel +//Container class +Solid* Solid::newSolid(const SolidType solidtype){ + switch (solidtype) { + case stainless: + return new Stainless(); + break; + case kapton: + return new Kapton(); + break; + case copper: + return new Copper(); + break; + case stainless_hopkins: + return new StainlessHopkins(); + break; + default: + FATAL("Undefined solid type"); + break; + } + return nullptr; +} + +vd StainlessHopkins::c(const vd& T) const{return 490*pow(T,0);} +d StainlessHopkins::c(const d& T) const{return 490;} +vd StainlessHopkins::rho(const vd& T) const{return 7900*pow(T,0);} +d StainlessHopkins::rho(const d& T) const{return 7900;} +vd StainlessHopkins::kappa(const vd& T) const{return 14.9*pow(T,0);} +d StainlessHopkins::kappa(const d& T) const{return 14.9;} + + +vd Stainless::c(const vd& T) const { + vd arg=1.7054e-6*pow(T,-0.88962)+22324.0/pow(T,6); + return pow(arg,-1.0/3.0)+15/T; +} +d Stainless::c(const d& T) const { + d arg=1.7054e-6*pow(T,-0.88962)+22324.0/pow(T,6); + return pow(arg,-1.0/3.0)+15/T; +} +vd Stainless::rho(const vd& T) const { + return 8274.55-1055.23*exp(-1.0*pow((T-273.15-2171.05)/2058.08,2)); +} +d Stainless::rho(const d& T) const { + return 8274.55-1055.23*exp(-1.0*pow((T-273.15-2171.05)/2058.08,2)); +} +vd Stainless::kappa(const vd& T) const { + return pow(266800.0*pow(T,-5.2)+0.21416*pow(T,-1.6),-0.25); +} +d Stainless::kappa(const d& T) const { + return pow(266800.0*pow(T,-5.2)+0.21416*pow(T,-1.6),-0.25); +} + +//Copper +vd Copper::kappa(const vd& T) const { + return 398.0-0.567*(T-300.0); +} +d Copper::kappa(const d& T) const { return 398.0-0.567*(T-300.0); } +vd Copper::c(const vd& T) const {return 420.0*pow(T,0);} +d Copper::c(const d& T) const {return 420.0*pow(T,0);} +vd Copper::rho(const vd& T) const {return 9000.0*pow(T,0);} +d Copper::rho(const d& T) const {return 9000.0*pow(T,0);} + + diff --git a/src/material/solid.h b/src/material/solid.h new file mode 100644 index 0000000..aece917 --- /dev/null +++ b/src/material/solid.h @@ -0,0 +1,78 @@ +#pragma once +#ifndef SOLID_H +#define SOLID_H + +#include "tasmet_types.h" +#include "tasmet_enum.h" + + +class Solid{ +public: + DECLARE_ENUM(SolidType,stainless,stainless_hopkins,copper,kapton); +private: + SolidType _solidtype; +protected: + Solid(const SolidType); +public: + Solid(const Solid&) = delete; + Solid& operator=(const Solid&) = delete; + operator SolidType() const {return _solidtype;} + + Solid* newSolid(const SolidType); + + virtual vd kappa(const vd& T) const=0; + virtual d kappa(const d& T) const=0; + virtual vd c(const vd& T) const=0; + virtual d c(const d& T) const=0; + virtual vd rho(const vd& T) const=0; + virtual d rho(const d& T) const=0; + virtual ~Solid(); + +}; + +class StainlessHopkins:public Solid +{ +public: + StainlessHopkins(): Solid(stainless_hopkins){} + vd kappa(const vd& T) const; + d kappa(const d& T) const; + vd c(const vd& T) const; + d c(const d& T) const; + vd rho(const vd& T) const; + d rho(const d& T) const; +}; + +class Stainless: public Solid{ +public: + Stainless(): Solid(stainless){} + vd kappa(const vd& T) const; + d kappa(const d& T) const; + vd c(const vd& T) const; + d c(const d& T) const; + vd rho(const vd& T) const; + d rho(const d& T) const; + +}; +class Copper: public Solid{ +public: + Copper(): Solid(copper){} + vd kappa(const vd& T) const; + d kappa(const d& T) const; + vd c(const vd& T) const; + d c(const d& T) const; + vd rho(const vd& T) const; + d rho(const d& T) const; +}; + +class Kapton: public Solid { +public: + Kapton():Solid(kapton){} + vd kappa(const vd& T) const { return 0.2*(1.0-exp(-T/100.0));} + d kappa(const d& T) const { return 0.2*(1.0-exp(-T/100.0));} + vd c(const vd& T) const {return 3.64*T;} + d c(const d& T) const {return 3.64*T;} + vd rho(const vd& T) const {return 1445.0-0.085*T;} + d rho(const d& T) const {return 1445.0-0.085*T;} +}; + +#endif /* SOLID_H */ diff --git a/src/tasmet_consolecolors.h b/src/tasmet_consolecolors.h new file mode 100644 index 0000000..2a7c4c7 --- /dev/null +++ b/src/tasmet_consolecolors.h @@ -0,0 +1,27 @@ +// consolecolors.h +// +// Author: J.A. de Jong +// +// Description: +// Print text from C++ to stdout in color +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef CONSOLECOLORS_H +#define CONSOLECOLORS_H +#include + +#ifndef SWIG +#define red "\e[31m" +#define green "\e[32m" +#define def " \e[39m" + +#endif // SWIG + +// Command to clear the content of a console +inline void clearConsole(){ + std::cout << "\033c" << std::endl; +} + +#endif // CONSOLECOLORS_H +////////////////////////////////////////////////////////////////////// + diff --git a/src/tasmet_enum.h b/src/tasmet_enum.h new file mode 100644 index 0000000..ac8825d --- /dev/null +++ b/src/tasmet_enum.h @@ -0,0 +1,61 @@ +// tasmet_enum.h +// +// Author: J.A. de Jong +// +// Description: Copied and adapted code from +// http://stackoverflow.com/questions/5093460/how-to- +// convert-an-enum-type-variable-to-a-string + +// Used to automatically couple a enum to the describing string. +// Use it as follows: +// +// #include "tasmet_enum.h" +// DECLARE_ENUM(OS_type, Windows, Linux, Apple) + +// void main() { +// cout << OS_typeToString(Windows) << endl; +// if (StringToOS_type("Linux") == Linux) { cout << "yay!" << endl; } +// cout << "Max: " << MAX_NUMBER_OF_OS_type << endl; +// } + +#pragma once +#ifndef TASMET_ENUM_H +#define TASMET_ENUM_H +#include +#include + +// Search and remove whitespace from both ends of the string +inline std::string TrimEnumString(const std::string &s) +{ + std::string::const_iterator it = s.begin(); + while (it != s.end() && isspace(*it)) { it++; } + std::string::const_reverse_iterator rit = s.rbegin(); + while (rit.base() != it && isspace(*rit)) { rit++; } + return std::string(it, rit.base()); +} +inline void SplitEnumArgs(const char* szArgs, std::string Array[], int nMax) +{ + std::stringstream ss(szArgs); + std::string strSub; + int nIdx = 0; + while (ss.good() && (nIdx < nMax)) { + getline(ss, strSub, ','); + Array[nIdx] = TrimEnumString(strSub); + nIdx++; + } +}; +#define DECLARE_ENUM(ename, ...) \ + enum ename { __VA_ARGS__, MAX_NUMBER_OF_##ename }; \ + static std::string ename##Strings[MAX_NUMBER_OF_##ename]; \ + inline const char* ename##ToString(ename e) { \ + if (ename##Strings[0].empty()) { SplitEnumArgs(#__VA_ARGS__, ename##Strings, MAX_NUMBER_OF_##ename); } \ + return ename##Strings[e].c_str(); \ + } \ + inline ename StringTo##ename(const char* szEnum) { \ + for (int i = 0; i < MAX_NUMBER_OF_##ename; i++) { if (ename##Strings[i] == szEnum) { return (ename)i; } } \ + return MAX_NUMBER_OF_##ename; \ + } + + +#endif // TASMET_ENUM_H +////////////////////////////////////////////////////////////////////// diff --git a/src/tasmet_exception.cpp b/src/tasmet_exception.cpp new file mode 100644 index 0000000..45e9e1a --- /dev/null +++ b/src/tasmet_exception.cpp @@ -0,0 +1,18 @@ +// tasmet_exception.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// Implementation of exception TasMETBadAlloc +////////////////////////////////////////////////////////////////////// + +#include "tasmet_exception.h" + +const char* TasMETBadAlloc::what() const throw() { + return "Error: memory allocation failed. " + "Please make sure enough memory is available and restart the application"; +} + + +////////////////////////////////////////////////////////////////////// + diff --git a/src/tasmet_exception.h b/src/tasmet_exception.h new file mode 100644 index 0000000..44a455c --- /dev/null +++ b/src/tasmet_exception.h @@ -0,0 +1,26 @@ +// tasmet_exception.h +// +// Author: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef TASMET_EXCEPTION_H +#define TASMET_EXCEPTION_H + +#include +#include + +class TaSMETError : public std::runtime_error { +public: + TaSMETError(const std::string& msg = "") : std::runtime_error(msg) {} +}; + +class TasMETBadAlloc: public std::bad_alloc { + virtual const char* what() const throw(); +}; + +#endif // TASMET_EXCEPTION_H +////////////////////////////////////////////////////////////////////// + diff --git a/src/tasmet_io.h b/src/tasmet_io.h new file mode 100644 index 0000000..96181ad --- /dev/null +++ b/src/tasmet_io.h @@ -0,0 +1,25 @@ +// tasmet_io.h +// +// Author: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef TASMET_IO_H +#define TASMET_IO_H +#include +#include + +//Spoiling global namespace with often used functions and variables +using std::cout; /* Output to stdout */ +using std::endl; +using std::setprecision; +using std::setiosflags; +using std::cin; /* Input from stdin */ +using std::ostream; +using std::cerr; +using std::ios; + +#endif // TASMET_IO_H +////////////////////////////////////////////////////////////////////// diff --git a/src/tasmet_tracer.cpp b/src/tasmet_tracer.cpp new file mode 100644 index 0000000..d7f1d5f --- /dev/null +++ b/src/tasmet_tracer.cpp @@ -0,0 +1,13 @@ +// tasmet_tracer.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// + +#include "tasmet_tracer.h" + + + +////////////////////////////////////////////////////////////////////// diff --git a/src/tasmet_tracer.h b/src/tasmet_tracer.h new file mode 100644 index 0000000..ed16978 --- /dev/null +++ b/src/tasmet_tracer.h @@ -0,0 +1,89 @@ +// tasmet_tracer.h +// +// Author: J.A. de Jong +// Description: +// When this file is included, we are able to use the following macro's: +// TRACE(10, "Test"); +// +// double x=2; +// VARTRACE(2,x); +// +// WARN("Warning message"); +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef TASMET_TRACER_H +#define TASMET_TRACER_H +#include + +// **************************************** Tracer code +// If PP variable TRACER is not defined, we automatically set it on. +#ifndef TRACER +static_assert(false,"TRACER macro not defined"); +#endif + +// Not so interesting part +#define rawstr(x) #x +#define namestr(x) rawstr(x) +#define annestr(x) namestr(x) +#define FILEWITHOUTPATH ( strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : strrchr(__FILE__, '\\') ? strrchr(__FILE__, '\\') + 1 : __FILE__ ) +#define POS FILEWITHOUTPATH << ":" << __LINE__ << ": " +// End not so interesting part + +#define RAWWARNING(a) std::cout << red << a << def << "\n"; +#define WARN(a) RAWWARNING(POS << "WARNING: " << a) + +// SHOULD NOT BE USED IN A LIBRARY!! +#define FATAL(a) WARN(a); abort(); + +#if TRACER == 1 +#include "tasmet_consolecolors.h" +#include + +// This variable can be used to increase the trace level of one +// specific file at compile time +#ifndef TRACERPLUS +#define TRACERPLUS (0) +#endif + + +// Initialize MAXTRACELEVEL and BUILDINTRACELEVEL +#ifndef MAXTRACELEVEL +#define MAXTRACELEVEL (5000) // To which the TRACELEVEL initially is set +#endif + +// Define this preprocessor definition to overwrite +// Use -O flag for compiler to remove the dead functions! +// In that case all cout's for TRACE() are removed from code +#ifndef BUILDINTRACELEVEL +#define BUILDINTRACELEVEL (-10) +#endif + +extern int tasmet_tracer_level; +// Use this preprocessor command to introduce one tasmet_tracer_level integer per unit +/* Introduce one static logger */ +// We trust that the compiler will eliminate 'dead code', which means +// that if variable BUILDINTRACERLEVEL is set, the inner if statement +// will not be reached. +#define WRITETRACE(level,a) \ + if(level>=BUILDINTRACELEVEL && level >=tasmet_tracer_level) { \ + std::cout << a << "\n"; \ + } + + +#define TRACE(l,a) WRITETRACE(l+TRACERPLUS,annestr(tasmet_tracer_level) << "-" << (l)<< "-" << POS << a) +#define VARTRACE(l,a) TRACE(l,#a " = " << a) + +#define INITTRACE(ll) \ + std::cout << "INITTRACE with tracelevel " << ll << " for " << annestr(tasmet_tracer_level) << "\n"; \ + tasmet_tracer_level=ll; +#else // TRACER !=1 + +#define VARTRACE(l,a) +#define TRACE(l,a) +#define INITTRACE(a) +#define WARN(a) + + +#endif // ######################################## TRACER ==1 +#endif // TASMET_TRACER_H +////////////////////////////////////////////////////////////////////// diff --git a/src/tasmet_types.h b/src/tasmet_types.h new file mode 100644 index 0000000..94e1729 --- /dev/null +++ b/src/tasmet_types.h @@ -0,0 +1,81 @@ +// tasmet_types.h +// +// Author: J.A. de Jong +// +// Description: +// Typedefs and namespace pollution for stuff that is almost always +// needed. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef TASMET_TYPES_H +#define TASMET_TYPES_H +#include + +// We require C++ +#ifndef __cplusplus +#error The c++ compiler should be used. +#endif + +// Some header files I (almost) always need +#include // C++ std vector +#include // C++ string class +#include +#include // Our linear algebra package + +typedef size_t us; /* Size type I always use */ + +#define fillwith arma::fill // Fill vector or matrix with ... + +// To change the whole code to 32-bit floating points, change this to +// float. +#if TASMET_FLOAT == 32 +// typedef float d; /* Shortcut for double */ +#elif TASMET_FLOAT == 64 +typedef double d; /* Shortcut for double */ +#else +static_assert(false,"TASMET_FLOAT should be either 32 or 64"); +#endif + +typedef std::complex c; /* Shortcut for c++ complex number */ + +using arma::abs; +using std::string; /* C++ Strings */ +using std::vector; // C++ Vectors + +/* Armadillo functions and objects we often need */ +using arma::sin; +using arma::cos; +using arma::exp; +using arma::ones; +using arma::zeros; +using arma::eye; +using arma::inv; +using arma::det; +using arma::log_det; +using arma::norm; +using arma::span; +using arma::linspace; +using arma::diagmat; + +// Constants +// I need these numbers so often, that they can be in the global namespace + +const c I(0,1); //Complex unity +const c sqI=sqrt(I); +const d sq2=sqrt(2.0); +const c sqmI=sqrt(-1.0*I); +const c minI(0,-1); //Minus complex unity +const d number_pi=arma::datum::pi; // The number 3.14159265359.. + + +typedef arma::Col vd; /* Column vector of doubles */ +typedef arma::Col vc; /* Column vector of complex numbers */ +typedef arma::Mat dmat; /* (Dense) Matrix of doubles */ +typedef arma::Mat cmat; /* (Dense) matrix of complex numbers */ + +typedef arma::SpMat sdmat; /* Sparse matrix of doubles */ + + +#endif // TASMET_TYPES_H +////////////////////////////////////////////////////////////////////// + diff --git a/src/tasmet_utils.h b/src/tasmet_utils.h new file mode 100644 index 0000000..8db2a96 --- /dev/null +++ b/src/tasmet_utils.h @@ -0,0 +1,71 @@ +// utils.h +// +// Author: J.A. de Jong +// +// Description: +// Some generic utils. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef UTILS_H +#define UTILS_H +#include +#include +#include "tracer.h" +#include "vtypes.h" +#include +#include +#include + + +// Purge a vector of components +template +void purge(std::vector& vec){ + TRACE(10,"purge(vector)"); + for (T& it: vec){ + delete it; + it=nullptr; + } + vec.clear(); +} + +// Purge a vector of components +template +void purge(std::map& map){ + TRACE(10,"purge(map)"); + for (auto& it: map){ + delete it.second; + it.second=nullptr; + } + map.clear(); +} + +template +const T& min(const T& x,const T& y) { + return x<=y? x : y; +} + +template +const T& max(const T& x,const T& y) { + return x<=y? y : x; +} + +template +SegType* copySeg(const SegType& t,const Sys& sys) { + SegType* newt=new SegType(t); + if(!newt){ + WARN("Copying " << typeid(t).name() << "failed!"); + } + return newt; +} // copySeg +template +void makeNormal(T& c) { + for(auto& val: c) + if(!std::isnormal(val)) + val=0; +} + +} // namespace utils + + +////////////////////////////////////////////////////////////////////// +