diff --git a/CMakeLists.txt b/CMakeLists.txt index 04425c743..ed4cc92a7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -64,6 +64,7 @@ set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/Code/CMake" #----------------------------------------------------------------------------- #----------------------------------------------------------------------------- + # CMake Includes include(CheckLibraryExists) include(GetPrerequisites) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 01e247d12..a620a268e 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -611,18 +611,40 @@ class faceType double qmTRI3 = 2.0/3.0; }; +/// @brief Store options for output types. +// +struct OutputOptions { + bool boundary_integral = false; + bool spatial = false; + bool volume_integral = false; + + bool no_options_set() { + return !(boundary_integral | spatial | volume_integral); + } + + void set_option(consts::OutputType type, bool value) + { + if (type == consts::OutputType::boundary_integral) { + boundary_integral = value; + } else if (type == consts::OutputType::spatial) { + spatial = value; + } else if (type == consts::OutputType::volume_integral) { + volume_integral = value; + } + } +}; + /// @brief Declared type for outputed variables // class outputType { public: - // Is this output suppose to be written into VTK, boundary, vol - std::vector wtn{false, false, false}; + // Options to write various output types. + OutputOptions options; // The group that this belong to (one of outType_*) - consts::OutputType grp = consts::OutputType::outGrp_NA; - //int grp; + consts::OutputNameType grp = consts::OutputNameType::outGrp_NA; // Length of the outputed variable int l = 0; diff --git a/Code/Source/solver/all_fun.cpp b/Code/Source/solver/all_fun.cpp index 89362b3f1..68c86ceb4 100644 --- a/Code/Source/solver/all_fun.cpp +++ b/Code/Source/solver/all_fun.cpp @@ -346,7 +346,8 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array equation_name_to_type = { }; +// Map from output type string to OutputType. +// +const std::map output_type_name_to_type = { + {"B_INT", OutputType::boundary_integral}, + {"Boundary_integral", OutputType::boundary_integral}, + {"Spatial", OutputType::spatial}, + {"V_INT", OutputType::volume_integral}, + {"Volume_integral", OutputType::volume_integral} +}; + const std::map mesh_generator_name_to_type = { {"Tetgen", MeshGeneratorType::RMSH_TETGEN}, {"Meshsim", MeshGeneratorType::RMSH_MESHSIM} diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 754133d57..b4d2ff173 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -61,7 +61,6 @@ int enum_int(T value) return static_cast(value); } - /// Check if a value is set to infinity. template bool present(T value) @@ -332,7 +331,7 @@ enum class MeshGeneratorType /// Map for string to MeshGeneratorType. extern const std::map mesh_generator_name_to_type; -enum class OutputType +enum class OutputNameType { outGrp_NA = 500, outGrp_A = 501, @@ -391,6 +390,17 @@ enum class OutputType out_CGInv1 = 572 }; +/// @brief Simulation output file types. +// +enum class OutputType +{ + boundary_integral, + spatial, + volume_integral +}; + +extern const std::map output_type_name_to_type; + /// @brief Possible physical properties. Current maxNPror is 20. // enum class PhysicalProperyType diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 09bc061a5..08a92820e 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -942,7 +942,7 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: { using namespace consts; - #define n_dist_eq + #define ndist_eq #ifdef dist_eq DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); @@ -1122,11 +1122,14 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: for (int iOut = 0; iOut < lEq.nOutput; iOut++) { auto& output = lEq.output[iOut]; - cm.bcast(cm_mod, output.wtn); cm.bcast_enum(cm_mod, &output.grp); cm.bcast(cm_mod, &output.o); cm.bcast(cm_mod, &output.l); cm.bcast(cm_mod, output.name); + + cm.bcast(cm_mod, &output.options.boundary_integral); + cm.bcast(cm_mod, &output.options.spatial); + cm.bcast(cm_mod, &output.options.volume_integral); } // Distribute BC information diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 68f3f6431..3d80dcc87 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -47,7 +47,7 @@ namespace post { void all_post(Simulation* simulation, Array& res, const Array& lY, const Array& lD, - consts::OutputType outGrp, const int iEq) + consts::OutputNameType outGrp, const int iEq) { using namespace consts; @@ -66,14 +66,14 @@ void all_post(Simulation* simulation, Array& res, const Array& l auto& msh = com_mod.msh[iM]; Array tmpV(maxNSD,msh.nNo); - if (outGrp == OutputType::outGrp_WSS || outGrp == OutputType::outGrp_trac) { + if (outGrp == OutputNameType::outGrp_WSS || outGrp == OutputNameType::outGrp_trac) { bpost(simulation, msh, tmpV, lY, lD, outGrp); for (int a = 0; a < com_mod.msh[iM].nNo; a++) { int Ac = msh.gN(a); res.set_col(Ac, tmpV.col(a)); } - } else if (outGrp == OutputType::outGrp_J) { + } else if (outGrp == OutputNameType::outGrp_J) { Array tmpV(1,msh.nNo); Vector tmpVe(msh.nEl); tpost(simulation, msh, 1, tmpV, tmpVe, lD, lY, iEq, outGrp); @@ -83,7 +83,7 @@ void all_post(Simulation* simulation, Array& res, const Array& l res(0,Ac) = tmpV(0,a); } - } else if (outGrp == OutputType::outGrp_mises) { + } else if (outGrp == OutputNameType::outGrp_mises) { Array tmpV(1,msh.nNo); Vector tmpVe(msh.nEl); tpost(simulation, msh, 1, tmpV, tmpVe, lD, lY, iEq, outGrp); @@ -93,7 +93,7 @@ void all_post(Simulation* simulation, Array& res, const Array& l res(0,Ac) = tmpV(0,a); } - } else if (outGrp == OutputType::outGrp_divV) { + } else if (outGrp == OutputNameType::outGrp_divV) { Array tmpV(1,msh.nNo); div_post(simulation, msh, tmpV, lY, lD, iEq); res = 0.0; @@ -117,7 +117,7 @@ void all_post(Simulation* simulation, Array& res, const Array& l /// Here t is stress tensor: t = \mu (grad(u) + grad(u)^T) // void bpost(Simulation* simulation, const mshType& lM, Array& res, const Array& lY, const Array& lD, - consts::OutputType outGrp) + consts::OutputNameType outGrp) { using namespace consts; @@ -132,7 +132,7 @@ void bpost(Simulation* simulation, const mshType& lM, Array& res, const dmsg << "outGrp: " << outGrp; #endif - if ((outGrp != OutputType::outGrp_WSS) && (outGrp != OutputType::outGrp_trac)) { + if ((outGrp != OutputNameType::outGrp_WSS) && (outGrp != OutputNameType::outGrp_trac)) { throw std::runtime_error("Invalid output group. Correction is required in BPOST"); } @@ -315,12 +315,12 @@ void bpost(Simulation* simulation, const mshType& lM, Array& res, const taue = Tdn - ndTdn*nV; Vector lRes(maxNSD); - if (outGrp == OutputType::outGrp_WSS) { + if (outGrp == OutputNameType::outGrp_WSS) { for (int i = 0; i < nsd; i++) { lRes(i) = -taue(i); } - } else if (outGrp == OutputType::outGrp_trac) { + } else if (outGrp == OutputNameType::outGrp_trac) { for (int i = 0; i < nsd; i++) { lRes(i) = p*nV(i) - Tdn(i); } @@ -858,7 +858,7 @@ void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const } void post(Simulation* simulation, const mshType& lM, Array& res, const Array& lY, const Array& lD, - consts::OutputType outGrp, const int iEq) + consts::OutputNameType outGrp, const int iEq) { using namespace consts; auto& com_mod = simulation->com_mod; @@ -886,7 +886,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const A // int nsd = com_mod.nsd; - if ((outGrp == OutputType::outGrp_eFlx) && (com_mod.dmnId.size() == 0)) { + if ((outGrp == OutputNameType::outGrp_eFlx) && (com_mod.dmnId.size() == 0)) { double rho = eq.dmn[0].prop[PhysicalProperyType::fluid_density]; for (int a = 0; a < lM.nNo; a++) { int Ac = lM.gN(a); @@ -963,7 +963,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const A // Vorticity calculation // - if (outGrp == OutputType::outGrp_vort) { + if (outGrp == OutputNameType::outGrp_vort) { for (int a = 0; a < eNoN; a++) { if (nsd == 2) { lRes(2) = lRes(2) + Nx(0,a)*yl(1,a)- Nx(1,a)*yl(0,a); @@ -975,7 +975,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const A } // Vortex Identification Criterion (lamda_ci) - } else if (outGrp == OutputType::outGrp_vortex) { + } else if (outGrp == OutputNameType::outGrp_vortex) { Array ux(nsd,nsd); for (int a = 0; a < eNoN; a++) { for (int i = 0; i < nsd; i++) { @@ -993,7 +993,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const A // Energy flux calculation // - } else if (outGrp == OutputType::outGrp_eFlx) { + } else if (outGrp == OutputNameType::outGrp_eFlx) { double rho = eq.dmn[cDmn].prop[PhysicalProperyType::fluid_density]; double p = 0.0; Vector u(nsd); @@ -1013,7 +1013,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const A // Heat flux calculation // - } else if (outGrp == OutputType::outGrp_hFlx) { + } else if (outGrp == OutputNameType::outGrp_hFlx) { double kappa = eq.dmn[cDmn].prop[PhysicalProperyType::conductivity]; int i = eq.s; @@ -1047,7 +1047,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const A // Strain tensor invariants calculation // - } else if (outGrp == OutputType::outGrp_stInv) { + } else if (outGrp == OutputNameType::outGrp_stInv) { Array ksix(nsd,nsd); Vector lRes(maxNSD); @@ -1079,7 +1079,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const A // Viscosity // - } else if (outGrp == OutputType::outGrp_Visc) { + } else if (outGrp == OutputNameType::outGrp_Visc) { Array ux(nsd,nsd); Vector lRes(maxNSD); @@ -1210,7 +1210,7 @@ void ppbin2vtk(Simulation* simulation) // Reproduces Fortran SHLPOST. // void shl_post(Simulation* simulation, const mshType& lM, const int m, Array& res, - Vector& resE, const Array& lD, const int iEq, consts::OutputType outGrp) + Vector& resE, const Array& lD, const int iEq, consts::OutputNameType outGrp) { using namespace consts; using namespace mat_fun; @@ -1531,13 +1531,13 @@ void shl_post(Simulation* simulation, const mshType& lM, const int m, Array S(3,3); @@ -1643,7 +1643,7 @@ void shl_post(Simulation* simulation, const mshType& lM, const int m, Array& res, Vector& resE, const Array& lD, - const Array& lY, const int iEq, consts::OutputType outGrp) + const Array& lY, const int iEq, consts::OutputNameType outGrp) { using namespace consts; using namespace mat_fun; @@ -1849,13 +1849,13 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array switch (outGrp) { // Jacobian := determinant of deformation gradient tensor - case OutputType::outGrp_J: + case OutputNameType::outGrp_J: resl(0) = detF; sE(e) = sE(e) + w*detF; break; // Deformation gradient tensor (F) - case OutputType::outGrp_F: + case OutputNameType::outGrp_F: if (nsd == 3) { resl(0) = F(0,0); resl(1) = F(0,1); @@ -1875,7 +1875,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array break; // Green-Lagrange strain tensor - case OutputType::outGrp_strain: + case OutputNameType::outGrp_strain: if (cPhys == EquationType::phys_lElas) { resl = ed; } else { @@ -1898,9 +1898,9 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array } break; - case OutputType::outGrp_stress: - case OutputType::outGrp_cauchy: - case OutputType::outGrp_mises: + case OutputNameType::outGrp_stress: + case OutputNameType::outGrp_cauchy: + case OutputNameType::outGrp_mises: Array sigma(nsd,nsd); Array S(nsd,nsd); @@ -1961,7 +1961,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array } // 2nd Piola-Kirchhoff stress tensor - if (outGrp == OutputType::outGrp_stress) { + if (outGrp == OutputNameType::outGrp_stress) { if (nsd == 3) { resl(0) = S(0,0); resl(1) = S(1,1); @@ -1976,7 +1976,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array } //Cauchy stress tensor - } else if (outGrp == OutputType::outGrp_cauchy) { + } else if (outGrp == OutputNameType::outGrp_cauchy) { if (nsd == 3) { resl(0) = sigma(0,0); resl(1) = sigma(1,1); @@ -1991,7 +1991,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array } // Von Mises stress - } else if (outGrp == OutputType::outGrp_mises) { + } else if (outGrp == OutputNameType::outGrp_mises) { double trS = mat_trace(sigma, nsd) / static_cast(nsd); for (int l = 0; l < nsd; l++) { sigma(l,l) = sigma(l,l) - trS; diff --git a/Code/Source/solver/post.h b/Code/Source/solver/post.h index 552e1fa58..f5c5a4385 100644 --- a/Code/Source/solver/post.h +++ b/Code/Source/solver/post.h @@ -37,10 +37,10 @@ namespace post { void all_post(Simulation* simulation, Array& res, const Array& lY, const Array& lD, - consts::OutputType outGrp, const int iEq); + consts::OutputNameType outGrp, const int iEq); void bpost(Simulation* simulation, const mshType& lM, Array& res, const Array& lY, const Array& lD, - consts::OutputType outGrp); + consts::OutputNameType outGrp); void div_post(Simulation* simulation, const mshType& lM, Array& res, const Array& lY, const Array& lD, const int iEq); @@ -51,15 +51,15 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const Array& lD, Vector& res); void post(Simulation* simulation, const mshType& lM, Array& res, const Array& lY, const Array& lD, - consts::OutputType outGrp, const int iEq); + consts::OutputNameType outGrp, const int iEq); void ppbin2vtk(Simulation* simulation); void shl_post(Simulation* simulation, const mshType& lM, const int m, Array& res, - Vector& resE, const Array& lD, const int iEq, consts::OutputType outGrp); + Vector& resE, const Array& lD, const int iEq, consts::OutputNameType outGrp); void tpost(Simulation* simulation, const mshType& lM, const int m, Array& res, Vector& resE, const Array& lD, - const Array& lY, const int iEq, consts::OutputType outGrp); + const Array& lY, const int iEq, consts::OutputNameType outGrp); }; diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 5cfaeb08c..8ca08df4b 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -2239,39 +2239,31 @@ void read_outputs(Simulation* simulation, EquationParameters* eq_params, eqType& } } - // These are the default values, we use the first nDef/nBDef outputs - // - for (int j = 0; j < 3; j++) { - for (int k = 0; k < nDOP[j+1]; k++) { - lEq.output[k].wtn[j] = true; - } - } - // First reading the outputs for VTK files and then for boundaries and last for the volume // int nOut = eq_params->outputs.size(); - std::map output_types = { {"Spatial", 0}, {"B_INT", 1}, {"Boundary_integral",1}, - {"V_INT", 2}, {"Volume_integral", 2}}; - int j = 0; + OutputType output_type; + for (int iOut = 0; iOut < nOut; iOut++) { auto& output_params = eq_params->outputs[iOut]; - auto output_type = output_params->type.value(); - if (output_type == "Alias") { + auto output_type_str = output_params->type.value(); + if (output_type_str == "Alias") { continue; } + try { - j = output_types.at(output_type); + output_type = output_type_name_to_type.at(output_type_str); } catch (const std::out_of_range& exception) { - throw std::runtime_error("Unknown output type '" + output_type + "'."); + throw std::runtime_error("Unknown output type '" + output_type_str + "'."); } auto& output_list = output_params->output_list; for (int i = 0; i < lEq.nOutput; i++) { - lEq.output[i].wtn[j] = output_params->get_output_value(lEq.output[i].name); + lEq.output[i].options.set_option(output_type, output_params->get_output_value(lEq.output[i].name)); if (lEq.output[i].name == "Vortex") { if (nsd != maxNSD) { - lEq.output[i].wtn[j] = false; + lEq.output[i].options.set_option(output_type, false); } } } diff --git a/Code/Source/solver/read_files.h b/Code/Source/solver/read_files.h index 656b190cb..45c1a6463 100644 --- a/Code/Source/solver/read_files.h +++ b/Code/Source/solver/read_files.h @@ -44,7 +44,7 @@ namespace read_files_ns { const int maxOutput = 22; using EquationNdop = std::array; - using EquationOutputs = std::array; + using EquationOutputs = std::array; using EquationPhys = std::vector; using EquationProps = std::array, 20>; diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index d59c2cb88..2d16b9648 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -67,7 +67,7 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); nDOP = {1, 1, 0, 0}; - outPuts[0] = OutputType::out_voltage; + outPuts[0] = OutputNameType::out_voltage; // Set solver parameters. read_ls(simulation, eq_params, SolverType::lSolver_CG, lEq); @@ -156,18 +156,18 @@ SetEquationPropertiesMapType set_equation_props = { nDOP = {12, 4, 3, 0}; outPuts = { - OutputType::out_velocity, - OutputType::out_pressure, - OutputType::out_WSS, - OutputType::out_displacement, - OutputType::out_energyFlux, - OutputType::out_traction, - OutputType::out_vorticity, - OutputType::out_vortex, - OutputType::out_strainInv, - OutputType::out_viscosity, - OutputType::out_divergence, - OutputType::out_acceleration + OutputNameType::out_velocity, + OutputNameType::out_pressure, + OutputNameType::out_WSS, + OutputNameType::out_displacement, + OutputNameType::out_energyFlux, + OutputNameType::out_traction, + OutputNameType::out_vorticity, + OutputNameType::out_vortex, + OutputNameType::out_strainInv, + OutputNameType::out_viscosity, + OutputNameType::out_divergence, + OutputNameType::out_acceleration }; } else { @@ -185,10 +185,10 @@ SetEquationPropertiesMapType set_equation_props = { if (pstEq) { nDOP = {2, 2, 0, 0}; - outPuts = { OutputType::out_displacement, OutputType::out_stress }; + outPuts = { OutputNameType::out_displacement, OutputNameType::out_stress }; } else { nDOP = {1, 1, 0, 0}; - outPuts[0] = OutputType::out_displacement; + outPuts[0] = OutputNameType::out_displacement; } } @@ -232,17 +232,17 @@ SetEquationPropertiesMapType set_equation_props = { nDOP = {11, 2, 3, 0}; outPuts = { - OutputType::out_velocity, - OutputType::out_pressure, - OutputType::out_WSS, - OutputType::out_traction, - OutputType::out_vorticity, - OutputType::out_vortex, - OutputType::out_strainInv, - OutputType::out_energyFlux, - OutputType::out_viscosity, - OutputType::out_divergence, - OutputType::out_acceleration + OutputNameType::out_velocity, + OutputNameType::out_pressure, + OutputNameType::out_WSS, + OutputNameType::out_traction, + OutputNameType::out_vorticity, + OutputNameType::out_vortex, + OutputNameType::out_strainInv, + OutputNameType::out_energyFlux, + OutputNameType::out_viscosity, + OutputNameType::out_divergence, + OutputNameType::out_acceleration }; // Set solver parameters. @@ -267,9 +267,9 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); nDOP = {3,1,1,0}; - outPuts = {OutputType::out_temperature, - OutputType::out_heatFlux, - OutputType::out_velocity}; + outPuts = {OutputNameType::out_temperature, + OutputNameType::out_heatFlux, + OutputNameType::out_velocity}; // Set solver parameters. read_ls(simulation, eq_params, SolverType::lSolver_GMRES, lEq); @@ -294,7 +294,7 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); nDOP = {2,1,1,0}; - outPuts = {OutputType::out_temperature, OutputType::out_heatFlux}; + outPuts = {OutputNameType::out_temperature, OutputNameType::out_heatFlux}; // Set solver parameters. read_ls(simulation, eq_params, SolverType::lSolver_CG, lEq); @@ -367,30 +367,30 @@ SetEquationPropertiesMapType set_equation_props = { nDOP = {22, 4, 2, 0}; outPuts = { - OutputType::out_velocity, - OutputType::out_pressure, - OutputType::out_displacement, - OutputType::out_mises, - - OutputType::out_WSS, - OutputType::out_traction, - OutputType::out_vorticity, - OutputType::out_vortex, - OutputType::out_strainInv, - OutputType::out_energyFlux, - OutputType::out_viscosity, - OutputType::out_absVelocity, - OutputType::out_stress, - OutputType::out_cauchy, - OutputType::out_strain, - OutputType::out_jacobian, - OutputType::out_defGrad, - OutputType::out_integ, - OutputType::out_fibDir, - OutputType::out_fibAlign, - - OutputType::out_divergence, - OutputType::out_acceleration + OutputNameType::out_velocity, + OutputNameType::out_pressure, + OutputNameType::out_displacement, + OutputNameType::out_mises, + + OutputNameType::out_WSS, + OutputNameType::out_traction, + OutputNameType::out_vorticity, + OutputNameType::out_vortex, + OutputNameType::out_strainInv, + OutputNameType::out_energyFlux, + OutputNameType::out_viscosity, + OutputNameType::out_absVelocity, + OutputNameType::out_stress, + OutputNameType::out_cauchy, + OutputNameType::out_strain, + OutputNameType::out_jacobian, + OutputNameType::out_defGrad, + OutputNameType::out_integ, + OutputNameType::out_fibDir, + OutputNameType::out_fibAlign, + + OutputNameType::out_divergence, + OutputNameType::out_acceleration }; // Set solver parameters. @@ -426,13 +426,13 @@ SetEquationPropertiesMapType set_equation_props = { if (eq_params->prestress.defined() && eq_params->prestress.value()) { nDOP = {3,2,0,0}; - outPuts = {OutputType::out_displacement, OutputType::out_stress, OutputType::out_strain}; + outPuts = {OutputNameType::out_displacement, OutputNameType::out_stress, OutputNameType::out_strain}; } else { nDOP = {8,2,0,0}; outPuts = { - OutputType::out_displacement, OutputType::out_mises, OutputType::out_stress, - OutputType::out_strain, OutputType::out_velocity, OutputType::out_acceleration, - OutputType::out_integ, OutputType::out_jacobian + OutputNameType::out_displacement, OutputNameType::out_mises, OutputNameType::out_stress, + OutputNameType::out_strain, OutputNameType::out_velocity, OutputNameType::out_acceleration, + OutputNameType::out_integ, OutputNameType::out_jacobian }; } @@ -469,7 +469,7 @@ SetEquationPropertiesMapType set_equation_props = { } nDOP = {3, 1, 0, 0}; - outPuts = {OutputType::out_displacement, OutputType::out_velocity, OutputType::out_acceleration }; + outPuts = {OutputNameType::out_displacement, OutputNameType::out_velocity, OutputNameType::out_acceleration }; lEq.ls.relTol = 0.2; @@ -503,15 +503,15 @@ SetEquationPropertiesMapType set_equation_props = { nDOP = {9,1,0,0}; outPuts = { - OutputType::out_displacement, - OutputType::out_stress, - OutputType::out_strain, - OutputType::out_jacobian, - OutputType::out_defGrad, - OutputType::out_velocity, - OutputType::out_integ, - OutputType::out_CGstrain, - OutputType::out_CGInv1 + OutputNameType::out_displacement, + OutputNameType::out_stress, + OutputNameType::out_strain, + OutputNameType::out_jacobian, + OutputNameType::out_defGrad, + OutputNameType::out_velocity, + OutputNameType::out_integ, + OutputNameType::out_CGstrain, + OutputNameType::out_CGInv1 }; // Set solver parameters. @@ -540,14 +540,14 @@ SetEquationPropertiesMapType set_equation_props = { nDOP = {8, 2, 3, 0}; outPuts = { - OutputType::out_velocity, - OutputType::out_pressure, - OutputType::out_WSS, - OutputType::out_vorticity, - OutputType::out_traction, - OutputType::out_strainInv, - OutputType::out_viscosity, - OutputType::out_divergence + OutputNameType::out_velocity, + OutputNameType::out_pressure, + OutputNameType::out_WSS, + OutputNameType::out_vorticity, + OutputNameType::out_traction, + OutputNameType::out_strainInv, + OutputNameType::out_viscosity, + OutputNameType::out_divergence }; // Set solver parameters. @@ -580,15 +580,15 @@ SetEquationPropertiesMapType set_equation_props = { if (eq_params->prestress.defined() && eq_params->prestress.value()) { nDOP = {4,2,0,0}; - outPuts = {OutputType::out_displacement, OutputType::out_stress, OutputType::out_cauchy, OutputType::out_strain}; + outPuts = {OutputNameType::out_displacement, OutputNameType::out_stress, OutputNameType::out_cauchy, OutputNameType::out_strain}; //simulation->com_mod.pstEq = true; } else { nDOP = {12,2,0,0}; outPuts = { - OutputType::out_displacement, OutputType::out_mises, OutputType::out_stress, - OutputType::out_cauchy, OutputType::out_strain, OutputType::out_jacobian, - OutputType::out_defGrad, OutputType::out_integ, OutputType::out_fibDir, - OutputType::out_fibAlign, OutputType::out_velocity, OutputType::out_acceleration + OutputNameType::out_displacement, OutputNameType::out_mises, OutputNameType::out_stress, + OutputNameType::out_cauchy, OutputNameType::out_strain, OutputNameType::out_jacobian, + OutputNameType::out_defGrad, OutputNameType::out_integ, OutputNameType::out_fibDir, + OutputNameType::out_fibAlign, OutputNameType::out_velocity, OutputNameType::out_acceleration }; } @@ -625,20 +625,20 @@ SetEquationPropertiesMapType set_equation_props = { nDOP = {14, 2, 0, 0}; outPuts = { - OutputType::out_displacement, - OutputType::out_mises, - OutputType::out_stress, - OutputType::out_cauchy, - OutputType::out_strain, - OutputType::out_jacobian, - OutputType::out_defGrad, - OutputType::out_integ, - OutputType::out_fibDir, - OutputType::out_fibAlign, - OutputType::out_velocity, - OutputType::out_pressure, - OutputType::out_acceleration, - OutputType::out_divergence + OutputNameType::out_displacement, + OutputNameType::out_mises, + OutputNameType::out_stress, + OutputNameType::out_cauchy, + OutputNameType::out_strain, + OutputNameType::out_jacobian, + OutputNameType::out_defGrad, + OutputNameType::out_integ, + OutputNameType::out_fibDir, + OutputNameType::out_fibAlign, + OutputNameType::out_velocity, + OutputNameType::out_pressure, + OutputNameType::out_acceleration, + OutputNameType::out_divergence }; // Set solver parameters. diff --git a/Code/Source/solver/set_output_props.h b/Code/Source/solver/set_output_props.h index 69a9e9128..247428bc7 100644 --- a/Code/Source/solver/set_output_props.h +++ b/Code/Source/solver/set_output_props.h @@ -29,7 +29,7 @@ */ // The 'output_props_map' map defined here sets equation -// output properties from the given OutputType. +// output properties from the given OutputNameType. // #include @@ -41,46 +41,46 @@ /// output.l - Length of the outputed variable /// output.name - The name to be used for the output and also in input file // -using OutputProps = std::tuple; +using OutputProps = std::tuple; /// @brief Reproduces Fortran READOUTPUTS. // -std::map output_props_map = +std::map output_props_map = { // ----------------------------------------------------------------- // output group o l name // ----------------------------------------------------------------- - {OutputType::out_absVelocity, std::make_tuple(OutputType::outGrp_absV, 0, nsd, "Absolute_velocity") }, - {OutputType::out_acceleration, std::make_tuple(OutputType::outGrp_A, 0, nsd, "Acceleration") }, - {OutputType::out_cauchy, std::make_tuple(OutputType::outGrp_cauchy, 0, com_mod.nsymd, "Cauchy_stress") }, + {OutputNameType::out_absVelocity, std::make_tuple(OutputNameType::outGrp_absV, 0, nsd, "Absolute_velocity") }, + {OutputNameType::out_acceleration, std::make_tuple(OutputNameType::outGrp_A, 0, nsd, "Acceleration") }, + {OutputNameType::out_cauchy, std::make_tuple(OutputNameType::outGrp_cauchy, 0, com_mod.nsymd, "Cauchy_stress") }, - {OutputType::out_CGInv1, std::make_tuple(OutputType::out_CGInv1, 0, 1, "CG_Strain_Trace") }, - {OutputType::out_CGstrain, std::make_tuple(OutputType::outGrp_C, 0, com_mod.nsymd, "CG_Strain") }, + {OutputNameType::out_CGInv1, std::make_tuple(OutputNameType::out_CGInv1, 0, 1, "CG_Strain_Trace") }, + {OutputNameType::out_CGstrain, std::make_tuple(OutputNameType::outGrp_C, 0, com_mod.nsymd, "CG_Strain") }, - {OutputType::out_defGrad, std::make_tuple(OutputType::outGrp_F, 0, nsd*nsd, "Def_grad") }, - {OutputType::out_displacement, std::make_tuple(OutputType::outGrp_D, 0, nsd, "Displacement") }, - {OutputType::out_divergence, std::make_tuple(OutputType::outGrp_divV, 0, 1, "Divergence") }, - {OutputType::out_energyFlux, std::make_tuple(OutputType::outGrp_eFlx, 0, nsd, "Energy_flux") }, + {OutputNameType::out_defGrad, std::make_tuple(OutputNameType::outGrp_F, 0, nsd*nsd, "Def_grad") }, + {OutputNameType::out_displacement, std::make_tuple(OutputNameType::outGrp_D, 0, nsd, "Displacement") }, + {OutputNameType::out_divergence, std::make_tuple(OutputNameType::outGrp_divV, 0, 1, "Divergence") }, + {OutputNameType::out_energyFlux, std::make_tuple(OutputNameType::outGrp_eFlx, 0, nsd, "Energy_flux") }, - {OutputType::out_fibAlign, std::make_tuple(OutputType::outGrp_fA, 0, 1, "Fiber_alignment") }, - {OutputType::out_fibDir, std::make_tuple(OutputType::outGrp_fN, 0, nsd, "Fiber_direction") }, - {OutputType::out_fibStrn, std::make_tuple(OutputType::outGrp_fS, 0, 1, "Fiber_shortening") }, + {OutputNameType::out_fibAlign, std::make_tuple(OutputNameType::outGrp_fA, 0, 1, "Fiber_alignment") }, + {OutputNameType::out_fibDir, std::make_tuple(OutputNameType::outGrp_fN, 0, nsd, "Fiber_direction") }, + {OutputNameType::out_fibStrn, std::make_tuple(OutputNameType::outGrp_fS, 0, 1, "Fiber_shortening") }, - {OutputType::out_heatFlux, std::make_tuple(OutputType::outGrp_hFlx, 0, nsd, "Heat_flux") }, - {OutputType::out_integ, std::make_tuple(OutputType::outGrp_I, 0, 1, nsd == 2 ? "Area" : "Volume") }, - {OutputType::out_jacobian, std::make_tuple(OutputType::outGrp_J, 0, 1, "Jacobian") }, - {OutputType::out_mises, std::make_tuple(OutputType::outGrp_mises, 0, 1, "VonMises_stress") }, - {OutputType::out_pressure, std::make_tuple(OutputType::outGrp_Y, nsd, 1, "Pressure") }, - {OutputType::out_strain, std::make_tuple(OutputType::outGrp_strain, 0, com_mod.nsymd, "Strain") }, - {OutputType::out_strainInv, std::make_tuple(OutputType::outGrp_stInv, 0, nsd, "Strain_invariants") }, - {OutputType::out_stress, std::make_tuple(OutputType::outGrp_stress, 0, com_mod.nsymd, "Stress") }, - {OutputType::out_temperature, std::make_tuple(OutputType::outGrp_Y, 0, 1, "Temperature") }, - {OutputType::out_traction, std::make_tuple(OutputType::outGrp_trac, 0, nsd, "Traction") }, - {OutputType::out_velocity, std::make_tuple(OutputType::outGrp_Y, 0, nsd, "Velocity") }, - {OutputType::out_viscosity, std::make_tuple(OutputType::outGrp_Visc, 0, 1, "Viscosity") }, - {OutputType::out_voltage, std::make_tuple(OutputType::outGrp_Y, 0, 1, "Action_potential") }, - {OutputType::out_vortex, std::make_tuple(OutputType::outGrp_vortex, 0, 1, "Vortex") }, - {OutputType::out_vorticity, std::make_tuple(OutputType::outGrp_vort, 0, maxNSD, "Vorticity") }, - {OutputType::out_WSS, std::make_tuple(OutputType::outGrp_WSS, 0, maxNSD, "WSS") } + {OutputNameType::out_heatFlux, std::make_tuple(OutputNameType::outGrp_hFlx, 0, nsd, "Heat_flux") }, + {OutputNameType::out_integ, std::make_tuple(OutputNameType::outGrp_I, 0, 1, nsd == 2 ? "Area" : "Volume") }, + {OutputNameType::out_jacobian, std::make_tuple(OutputNameType::outGrp_J, 0, 1, "Jacobian") }, + {OutputNameType::out_mises, std::make_tuple(OutputNameType::outGrp_mises, 0, 1, "VonMises_stress") }, + {OutputNameType::out_pressure, std::make_tuple(OutputNameType::outGrp_Y, nsd, 1, "Pressure") }, + {OutputNameType::out_strain, std::make_tuple(OutputNameType::outGrp_strain, 0, com_mod.nsymd, "Strain") }, + {OutputNameType::out_strainInv, std::make_tuple(OutputNameType::outGrp_stInv, 0, nsd, "Strain_invariants") }, + {OutputNameType::out_stress, std::make_tuple(OutputNameType::outGrp_stress, 0, com_mod.nsymd, "Stress") }, + {OutputNameType::out_temperature, std::make_tuple(OutputNameType::outGrp_Y, 0, 1, "Temperature") }, + {OutputNameType::out_traction, std::make_tuple(OutputNameType::outGrp_trac, 0, nsd, "Traction") }, + {OutputNameType::out_velocity, std::make_tuple(OutputNameType::outGrp_Y, 0, nsd, "Velocity") }, + {OutputNameType::out_viscosity, std::make_tuple(OutputNameType::outGrp_Visc, 0, 1, "Viscosity") }, + {OutputNameType::out_voltage, std::make_tuple(OutputNameType::outGrp_Y, 0, 1, "Action_potential") }, + {OutputNameType::out_vortex, std::make_tuple(OutputNameType::outGrp_vortex, 0, 1, "Vortex") }, + {OutputNameType::out_vorticity, std::make_tuple(OutputNameType::outGrp_vort, 0, maxNSD, "Vorticity") }, + {OutputNameType::out_WSS, std::make_tuple(OutputNameType::outGrp_WSS, 0, maxNSD, "WSS") } }; diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 26338671b..948d48e10 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -43,82 +43,131 @@ namespace txt_ns { -/// @brief This is to check/create the txt file +/// @brief Create a text file for writing boundary integral quantities. // -void cctxt(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const std::array& fName, const std::vector& wtn) +void create_boundary_integral_file(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const std::string& file_name) { int fid = 1; - if (com_mod.cm.slv(cm_mod)) { return; } - // i=1 are the boundary values and i=2 are volume values + bool file_exists = false; + std::string geometry; + + if (FILE *fp = fopen(file_name.c_str(), "r")) { + file_exists = true; + fclose(fp); + } + + // [NOTE] What is this all about? + if (com_mod.cTS != 0 && file_exists) { + //CALL TRIMFILE(cTS+3,fName(i)) + return; + } + + auto fp = fopen(file_name.c_str(), "w"); + fprintf(fp, "# svMultiPhysics boundary integral results file. \n"); + fprintf(fp, "# Quantities represent averaged scalar or flux values over each mesh face.\n"); + fprintf(fp, "# \n"); + fprintf(fp, "# Format \n"); + fprintf(fp, "# ------ \n"); + fprintf(fp, "# face areas: [list of mesh face areas] \n"); + fprintf(fp, "# step time [list of mesh face names] \n"); + fprintf(fp, "# [time step] [time] [list of computed values for each mesh face]\n"); + + // Write face areas. // - for (int i = 0; i < 2; i++) { - if (!wtn[i]) { - continue; + fprintf(fp, "face areas: "); + for (int iM = 0; iM < com_mod.nMsh; iM++) { + auto& msh = com_mod.msh[iM]; + for (int iFa = 0; iFa < msh.nFa; iFa++) { + auto& fa = msh.fa[iFa]; + fprintf(fp, " %4.3e ", fa.area); } + } + fprintf(fp, "\n"); - bool flag = false; - auto file_name = fName[i]; - if (FILE *fp = fopen(file_name.c_str(), "r")) { - flag = true; - fclose(fp); + // Write face names. + // + fprintf(fp, "step time "); + + for (int iM = 0; iM < com_mod.nMsh; iM++) { + auto& msh = com_mod.msh[iM]; + for (int iFa = 0; iFa < msh.nFa; iFa++) { + auto& fa = msh.fa[iFa]; + auto stmp = fa.name; + fprintf(fp, " %s ", stmp.c_str()); } + } - // [NOTE] What is this all about? - if (com_mod.cTS != 0 && flag) { - //CALL TRIMFILE(cTS+3,fName(i)) - continue; - } + fprintf(fp, "\n"); + fclose(fp); +} - auto fp = fopen(file_name.c_str(), "w"); +/// @brief Create a text file for writing volume integral quantities. +// +void create_volume_integral_file(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const std::string& file_name) +{ + int fid = 1; + if (com_mod.cm.slv(cm_mod)) { + return; + } - if (i == 0) { - for (int iM = 0; iM < com_mod.nMsh; iM++) { - auto& msh = com_mod.msh[iM]; + bool file_exists = false; + std::string geometry; - for (int iFa = 0; iFa < msh.nFa; iFa++) { - auto& fa = msh.fa[iFa]; - auto stmp = fa.name; - fprintf(fp, " %s ", stmp.c_str()); - } - } + if (FILE *fp = fopen(file_name.c_str(), "r")) { + file_exists = true; + fclose(fp); + } - fprintf(fp, "\n"); + // [NOTE] What is this all about? + if (com_mod.cTS != 0 && file_exists) { + //CALL TRIMFILE(cTS+3,fName(i)) + return; + } - for (int iM = 0; iM < com_mod.nMsh; iM++) { - auto& msh = com_mod.msh[iM]; - for (int iFa = 0; iFa < msh.nFa; iFa++) { - auto& fa = msh.fa[iFa]; - fprintf(fp, " %4.3e ", fa.area); - } - } + auto fp = fopen(file_name.c_str(), "w"); + fprintf(fp, "# svMultiPhysics volume integral results file. \n"); + fprintf(fp, "# Quantities represent averaged scalar or flux values over each volume mesh domain.\n"); + fprintf(fp, "# \n"); + fprintf(fp, "# Format \n"); + fprintf(fp, "# ------ \n"); + fprintf(fp, "# domain volumes: [list of mesh domain volumes] \n"); + fprintf(fp, "# step time [list of mesh domain names] \n"); + fprintf(fp, "# [time step] [time] [list of computed values for each mesh domain]\n"); + + // Write domain volumes. + // + fprintf(fp, "domain volumes: "); + for (int iDmn = 0; iDmn < lEq.nDmn; iDmn++) { + fprintf(fp, " %4.3e ", lEq.dmn[iDmn].v); + } + fprintf(fp, "\n"); - } else { - for (int iDmn = 0; iDmn < lEq.nDmn; iDmn++) { - std::string stmp; - if (lEq.dmn[iDmn].Id == -1) { - stmp = "ENTIRE"; - } - } + // Write domain names. + // + fprintf(fp, "step time "); - for (int iDmn = 0; iDmn < lEq.nDmn; iDmn++) { - std::string stmp; - } + for (int iDmn = 0; iDmn < lEq.nDmn; iDmn++) { + std::string stmp = "domain-" + std::to_string(lEq.dmn[iDmn].Id); + if (lEq.dmn[iDmn].Id == -1) { + stmp = "domain"; } - - fprintf(fp, "\n"); - fprintf(fp, "\n"); - fclose(fp); + fprintf(fp, "%s ", stmp.c_str()); } + fprintf(fp, "\n"); + fclose(fp); } /// \todo [NOTE] This is not fully implemented. // -void txt(Simulation* simulation, const bool flag) +// init_write - if true then this is the start of the simulation and +// so create a new file to initialize output. +// +void txt(Simulation* simulation, const bool init_write) { using namespace consts; using namespace utils; @@ -146,7 +195,7 @@ void txt(Simulation* simulation, const bool flag) if (cplBC.coupled) { bool ltmp = false; - if (!flag) { + if (!init_write) { if (cplBC.useGenBC) { set_bc::genBC_Integ_X(com_mod, cm_mod, "L"); @@ -167,7 +216,7 @@ void txt(Simulation* simulation, const bool flag) } if (com_mod.cm.mas(cm_mod) && !cplBC.useGenBC) { - if (flag) { + if (init_write) { if (com_mod.cTS == 0 || !ltmp) { //OPEN(fid, FILE=cplBC.saveName) //CLOSE(fid, STATUS='DELETE') @@ -195,12 +244,9 @@ void txt(Simulation* simulation, const bool flag) // Process eq outputs // - //dmsg << "Process eq outputs ... "; - //dmsg << "com_mod.nEq: " << com_mod.nEq; - for (int iEq = 0; iEq < com_mod.nEq; iEq++) { #ifdef debug_txt - dmsg << "---------- iEq " << iEq+1 << " ----------"; + dmsg << "---------- iEq " << iEq+1; #endif Array tmpV(maxNSD,tnNo); auto& eq = com_mod.eq[iEq]; @@ -210,15 +256,13 @@ void txt(Simulation* simulation, const bool flag) // Don't write it when it doesn't suppose to be written #ifdef debug_txt dmsg; - dmsg << "----- iOut " << iOut+1 << " -----"; + dmsg << "----- iOut " << iOut+1; #endif - std::vector wtn; - for (int i = 1; i < 3; i++) { - wtn.push_back( eq.output[iOut].wtn[i] ); - } - - if (std::count(wtn.begin(), wtn.end(), false) == wtn.size()) { + // Get options for computing boundary and volume integral quantities. + // + auto output_options = eq.output[iOut].options; + if (!output_options.boundary_integral && !output_options.volume_integral) { continue; } @@ -245,11 +289,11 @@ void txt(Simulation* simulation, const bool flag) #endif switch (oGrp) { - case OutputType::outGrp_NA: + case OutputNameType::outGrp_NA: throw std::runtime_error("NA outGrp in TXT"); break; - case OutputType::outGrp_A: + case OutputNameType::outGrp_A: for (int j = 0; j < tnNo; j++) { for (int i = 0; i < l; i++) { tmpV(i,j) = com_mod.An(i+s,j); @@ -257,7 +301,7 @@ void txt(Simulation* simulation, const bool flag) } break; - case OutputType::outGrp_Y: + case OutputNameType::outGrp_Y: for (int j = 0; j < tnNo; j++) { for (int i = 0; i < l; i++) { tmpV(i,j) = com_mod.Yn(i+s,j); @@ -269,7 +313,7 @@ void txt(Simulation* simulation, const bool flag) } break; - case OutputType::outGrp_D: + case OutputNameType::outGrp_D: for (int j = 0; j < tnNo; j++) { for (int i = 0; i < l; i++) { tmpV(i,j) = com_mod.Dn(i+s,j); @@ -277,9 +321,9 @@ void txt(Simulation* simulation, const bool flag) } break; - case OutputType::outGrp_WSS: - case OutputType::outGrp_vort: - case OutputType::outGrp_trac: + case OutputNameType::outGrp_WSS: + case OutputNameType::outGrp_vort: + case OutputNameType::outGrp_trac: post::all_post(simulation, tmpV, com_mod.Yn, com_mod.Dn, oGrp, iEq); for (int a = 0; a < tnNo; a++) { auto vec = tmpV.col(a, {0,nsd-1}); @@ -288,15 +332,15 @@ void txt(Simulation* simulation, const bool flag) l = 1; break; - case OutputType::outGrp_eFlx: - case OutputType::outGrp_hFlx: - case OutputType::outGrp_divV: - case OutputType::outGrp_J: - case OutputType::outGrp_mises: + case OutputNameType::outGrp_eFlx: + case OutputNameType::outGrp_hFlx: + case OutputNameType::outGrp_divV: + case OutputNameType::outGrp_J: + case OutputNameType::outGrp_mises: post::all_post(simulation, tmpV, com_mod.Yn, com_mod.Dn, oGrp, iEq); break; - case OutputType::outGrp_absV: + case OutputNameType::outGrp_absV: for (int a = 0; a < tnNo; a++) { for (int i = 0; i < l; i++) { tmpV(i,a) = com_mod.Yn(i,a) - com_mod.Yn(i+nsd+1,a); @@ -304,7 +348,7 @@ void txt(Simulation* simulation, const bool flag) } break; - case OutputType::outGrp_I: + case OutputNameType::outGrp_I: div = false; for (int a = 0; a < tnNo; a++) { for (int i = 0; i < l; i++) { @@ -314,96 +358,49 @@ void txt(Simulation* simulation, const bool flag) break; default: - throw std::runtime_error("Undefined output"); + throw std::runtime_error("Undefined output '" + std::to_string(static_cast(oGrp)) + "'"); } - // Write average/flux of variables at the boundaries + // Write boundary and volume integral quantities. // #ifdef debug_txt dmsg << "Write average/flux of variables ..."; #endif const auto& appPath = chnl_mod.appPath; - std::array fName; - fName[0] = eq.sym + "_" + output.name; - fName[1] = fName[0]; + std::string file_name = eq.sym + "_" + output.name; + std::string boundary_file_name; if (l == nsd) { - fName[0] = appPath + "B_" + fName[0] + "_flux.txt"; + boundary_file_name = appPath + "B_" + file_name + "_flux.txt"; } else { - fName[0] = appPath + "B_" + fName[0] + "_average.txt"; + boundary_file_name = appPath + "B_" + file_name + "_average.txt"; } - fName[1] = appPath + "V_" + fName[1] + "_flux.txt"; + std::string volume_file_name = appPath + "V_" + file_name + "_flux.txt"; + #ifdef debug_txt - dmsg << "flag: " << flag; - dmsg << "fName[0]: " << fName[0]; - dmsg << "fName[1]: " << fName[1]; + dmsg << "init_write: " << init_write; + dmsg << "boundary_file_name: " << boundary_file_name; + dmsg << "volume_file_name: " << volume_file_name; #endif - if (flag) { - cctxt(com_mod, cm_mod, eq, fName, wtn); - } else { - wtxt(com_mod, cm_mod, eq, l, fName, tmpV, wtn, div, pflag); - } - } - - // IB outputs - - // TODO:DaveP] not implemented. - // -/* - if (.NOT.ibFlag) continue; - if (ALLOCATED(tmpV)) DEALLOCATE(tmpV) - ALLOCATE (tmpV(maxnsd,ib.tnNo)) - for (int iOut=1, eq(iEq).nOutIB - // Don't write it when it doesn't suppose to be written - wtn = eq(iEq).outIB(iOut).wtn(2:3) - if (ALL(.NOT.wtn)) continue; - l = eq(iEq).outIB(iOut).l - s = eq(iEq).s + eq(iEq).outIB(iOut).o - e = s + l - 1 - - oGrp = eq(iEq).outIB(iOut).grp - div = .TRUE. - - SELECT case (oGrp) - case (outGrp_NA) - err = "NA outGrp in TXT" - case (outGrp_Y) - for (int a=1, ib.tnNo - tmpV(1:l,a) = ib.Yb(s:e,a)/REAL(cm.np(), KIND=RKIND) - { - case (outGrp_D) - for (int a=1, ib.tnNo - tmpV(1:l,a) = ib.Ubo(s:e,a)/REAL(cm.np(), KIND=RKIND) - { - case (outGrp_I) - div = .FALSE. - for (int a=1, ib.tnNo - tmpV(1:l,a) = 1. - { - case DEFAULT - err = "Undefined output" - END SELECT - - fName = eq(iEq).sym//"_"//TRIM(eq(iEq).outIB(iOut).name) - if (l .EQ. nsd) { - fName(1) = TRIM(appPath)//"IB_B_"//TRIM(fName(1))// - 2 "_flux.txt" - } else { - fName(1) = TRIM(appPath)//"IB_B_"//TRIM(fName(1))// - 2 "_average.txt" + if (init_write) { + if (output_options.boundary_integral) { + create_boundary_integral_file(com_mod, cm_mod, eq, boundary_file_name); + } + if (output_options.volume_integral) { + create_volume_integral_file(com_mod, cm_mod, eq, volume_file_name); } - fName(2) = TRIM(appPath)//"IB_V_"//TRIM(fName(2))// - 2 "_average.txt" - if (flag) { - CALL IB_CCTXT(eq(iEq), fName, wtn) - } else { - CALL IB_WTXT(eq(iEq), l, fName, tmpV, wtn, div) + } else { + if (output_options.boundary_integral) { + write_boundary_integral_data(com_mod, cm_mod, eq, l, boundary_file_name, tmpV, div, pflag); } - } -*/ + if (output_options.volume_integral) { + write_volume_integral_data(com_mod, cm_mod, eq, l, volume_file_name, tmpV, div, pflag); + } + } + } // ECG leads output if (com_mod.cm.idcm() == cm_mod.master) { @@ -432,89 +429,72 @@ void txt(Simulation* simulation, const bool flag) #endif } -/// @brief This is to write to txt file +/// @brief Compute boundary integral for quantities at the current time step +/// and write them to a text file. /// -/// \todo TODO:DaveP] not fully implemented. -/// -/// WARNING: There is an indexing problem here because 'm' is a length and not an index. +/// NOTE: Be carefu of a potential indexing problem here because 'm' is a length and not an index. // -void wtxt(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, const std::array& fName, - const Array& tmpV, const std::vector& wtn, const bool div, const bool pFlag) +void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, + const std::string file_name, const Array& tmpV, const bool div, const bool pFlag) { - return; - - #ifdef debug_wtxt + #define n_debug_write_boundary_integral_data + #ifdef debug_write_boundary_integral_data DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); + dmsg << "file_name: " << file_name; dmsg << "m: " << m; dmsg << "div: " << div; dmsg << "pFlag: " << pFlag; #endif + double time = com_mod.time; + int time_step = com_mod.cTS; + #ifdef debug_wtxt + dmsg << "time_step: " << time_step; + dmsg << "time: " << time; + #endif + int fid = 1; int nsd = com_mod.nsd; FILE* fp = nullptr; - for (int i = 0; i < 2; i++) { - if (!wtn[i]) { - continue; - } - - if (com_mod.cm.mas(cm_mod)) { - fp = fopen(fName[i].c_str(), "a"); - } + if (com_mod.cm.mas(cm_mod)) { + fp = fopen(file_name.c_str(), "a"); + fprintf(fp, " %d %.10e ", time_step, time); + } - if (i == 0) { - for (int iM = 0; iM < com_mod.nMsh; iM++) { - auto& msh = com_mod.msh[iM]; - bool lTH = false; + for (int iM = 0; iM < com_mod.nMsh; iM++) { + auto& msh = com_mod.msh[iM]; + bool lTH = false; - if (msh.nFs == 2) { - lTH = true; - } + if (msh.nFs == 2) { + lTH = true; + } - for (int iFa = 0; iFa < msh.nFa; iFa++) { - auto& fa = msh.fa[iFa]; - double tmp = 0.0; - - if (m == 1) { - if (div) { - tmp = fa.area; - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0) / tmp; - } else { - if (pFlag && lTH) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, true); - } else { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0); - } - } + for (int iFa = 0; iFa < msh.nFa; iFa++) { + auto& fa = msh.fa[iFa]; + double tmp = 0.0; - } else if (m == nsd) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, m-1); + if (m == 1) { + if (div) { + tmp = fa.area; + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0) / tmp; + } else { + if (pFlag && lTH) { + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, true); } else { - throw std::runtime_error("WTXT only accepts 1 and nsd"); - } - - if (com_mod.cm.mas(cm_mod)) { - fprintf(fp, " %.10e ", tmp); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0); } } - } - } else { - for (int iDmn = 0; iDmn < lEq.nDmn; iDmn++) { - auto& dmn = lEq.dmn[iDmn]; - double tmp = 0.0; + } else if (m == nsd) { + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, m-1); + } else { + throw std::runtime_error("WTXT only accepts 1 and nsd"); + } - if (div) { - tmp = dmn.v; - tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1) / tmp; - } else { - tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1, pFlag); - } - if (com_mod.cm.mas(cm_mod)) { - fprintf(fp, " %.10e ", tmp); - } + if (com_mod.cm.mas(cm_mod)) { + fprintf(fp, " %.10e ", tmp); } } @@ -525,4 +505,60 @@ void wtxt(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, } } +/// @brief Compute volume integral for quantities at the current time step +/// and write them to a text file. +/// +/// NOTE: Be carefu of a potential indexing problem here because 'm' is a length and not an index. +// +void write_volume_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, + const std::string file_name, const Array& tmpV, const bool div, const bool pFlag) +{ + #define n_debug_write_volume_integral_data + #ifdef debug_write_volume_integral_data + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "file_name: " << file_name; + dmsg << "m: " << m; + dmsg << "div: " << div; + dmsg << "pFlag: " << pFlag; + #endif + + double time = com_mod.time; + int time_step = com_mod.cTS; + #ifdef debug_wtxt + dmsg << "time_step: " << time_step; + dmsg << "time: " << time; + #endif + + int fid = 1; + int nsd = com_mod.nsd; + FILE* fp = nullptr; + + if (com_mod.cm.mas(cm_mod)) { + fp = fopen(file_name.c_str(), "a"); + fprintf(fp, " %d %.10e ", time_step, time); + } + + for (int iDmn = 0; iDmn < lEq.nDmn; iDmn++) { + auto& dmn = lEq.dmn[iDmn]; + double tmp = 0.0; + + if (div) { + tmp = dmn.v; + tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1) / tmp; + } else { + tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1, pFlag); + } + + if (com_mod.cm.mas(cm_mod)) { + fprintf(fp, " %.10e ", tmp); + } + } + + if (com_mod.cm.mas(cm_mod)) { + fprintf(fp, "\n"); + fclose(fp); + } +} + }; diff --git a/Code/Source/solver/txt.h b/Code/Source/solver/txt.h index 6fee31287..c232d758f 100644 --- a/Code/Source/solver/txt.h +++ b/Code/Source/solver/txt.h @@ -39,12 +39,17 @@ namespace txt_ns { -void cctxt(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const std::array& fName, const std::vector& wtn); +void create_boundary_integral_file(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const std::string& file_name); + +void create_volume_integral_file(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const std::string& file_name); void txt(Simulation* simulation, const bool flag); -void wtxt(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, const std::array& fName, - const Array& tmpV, const std::vector& wtn, const bool div, const bool pflag); +void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, + const std::string file_name, const Array& tmpV, const bool div, const bool pFlag); + +void write_volume_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, + const std::string file_name, const Array& tmpV, const bool div, const bool pFlag); }; diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index e841b45e3..536127df6 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -739,7 +739,7 @@ void read_vtus(Simulation* simulation, Array& lA, Array& lY, Arr for (int iOut = 0; iOut < eq.nOutput; iOut++) { auto& output = eq.output[iOut]; - if (!output.wtn[0]) { + if (!output.options.spatial) { continue; } @@ -752,9 +752,9 @@ void read_vtus(Simulation* simulation, Array& lA, Array& lY, Arr Array tmpGS; switch (oGrp) { - case OutputType::outGrp_A: - case OutputType::outGrp_Y: - case OutputType::outGrp_D: + case OutputNameType::outGrp_A: + case OutputNameType::outGrp_Y: + case OutputNameType::outGrp_D: if (l > 1) { tmpGS.resize(maxNSD,nNo); } else { @@ -791,7 +791,7 @@ void read_vtus(Simulation* simulation, Array& lA, Array& lY, Arr } switch (oGrp) { - case OutputType::outGrp_A: + case OutputNameType::outGrp_A: for (int i = 0; i < l; i++) { for (int j = 0; j < gtnNo; j++) { lA(i+s,j) = tmpGS(i,j); @@ -799,7 +799,7 @@ void read_vtus(Simulation* simulation, Array& lA, Array& lY, Arr } break; - case OutputType::outGrp_Y: + case OutputNameType::outGrp_Y: for (int i = 0; i < l; i++) { for (int j = 0; j < gtnNo; j++) { lY(i+s,j) = tmpGS(i,j); @@ -807,7 +807,7 @@ void read_vtus(Simulation* simulation, Array& lA, Array& lY, Arr } break; - case OutputType::outGrp_D: + case OutputNameType::outGrp_D: for (int i = 0; i < l; i++) { for (int j = 0; j < gtnNo; j++) { lD(i+s,j) = tmpGS(i,j); @@ -959,7 +959,7 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array tmpVe; switch (oGrp) { - case OutputType::outGrp_NA: + case OutputNameType::outGrp_NA: throw std::runtime_error("Undefined output grp in VTK"); break; - case OutputType::outGrp_A: + case OutputNameType::outGrp_A: for (int a = 0; a < msh.nNo; a++) { int Ac = msh.gN(a); for (int i = 0; i < l; i++) { @@ -1078,7 +1078,7 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array true + + true + true + + + + true + + fsils