From 1328c7a2fe6bfbf7731f1fa8de29f022af607b50 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Sun, 24 Sep 2023 17:21:30 -0700 Subject: [PATCH 01/15] adding darcy equation --- Code/Source/svFSI/darcy.cpp | 3 +++ Code/Source/svFSI/darcy.h | 8 ++++++++ 2 files changed, 11 insertions(+) create mode 100644 Code/Source/svFSI/darcy.cpp create mode 100644 Code/Source/svFSI/darcy.h diff --git a/Code/Source/svFSI/darcy.cpp b/Code/Source/svFSI/darcy.cpp new file mode 100644 index 00000000..aeac8360 --- /dev/null +++ b/Code/Source/svFSI/darcy.cpp @@ -0,0 +1,3 @@ +// +// Created by Zack on 9/19/2023. +// diff --git a/Code/Source/svFSI/darcy.h b/Code/Source/svFSI/darcy.h new file mode 100644 index 00000000..cd4f9ba9 --- /dev/null +++ b/Code/Source/svFSI/darcy.h @@ -0,0 +1,8 @@ +// +// Created by Zack on 9/19/2023. +// + +#ifndef SV_TOP_DARCY_H +#define SV_TOP_DARCY_H + +#endif //SV_TOP_DARCY_H From 813262d655c8eca0416ca5ea72fcf258f27e7d98 Mon Sep 17 00:00:00 2001 From: zasexton Date: Mon, 25 Sep 2023 15:27:29 -0700 Subject: [PATCH 02/15] darcy equation can be built but generates a trivially zero solution --- Code/Source/svFSI/CMakeLists.txt | 1 + Code/Source/svFSI/ComMod.h | 7 + Code/Source/svFSI/Parameters.cpp | 1 + Code/Source/svFSI/Parameters.h | 1 + Code/Source/svFSI/all_fun.cpp | 42 ++++ Code/Source/svFSI/consts.cpp | 1 + Code/Source/svFSI/consts.h | 11 +- Code/Source/svFSI/darcy.cpp | 290 ++++++++++++++++++++++++- Code/Source/svFSI/darcy.h | 25 ++- Code/Source/svFSI/eq_assem.cpp | 9 + Code/Source/svFSI/read_files.cpp | 38 +++- Code/Source/svFSI/set_equation_dof.h | 1 + Code/Source/svFSI/set_equation_props.h | 26 ++- Code/Source/svFSI/set_output_props.h | 3 +- 14 files changed, 443 insertions(+), 13 deletions(-) diff --git a/Code/Source/svFSI/CMakeLists.txt b/Code/Source/svFSI/CMakeLists.txt index 4bfd2164..56459255 100644 --- a/Code/Source/svFSI/CMakeLists.txt +++ b/Code/Source/svFSI/CMakeLists.txt @@ -106,6 +106,7 @@ set(CSRCS fft.h fft.cpp heatf.h heatf.cpp heats.h heats.cpp + darcy.h darcy.cpp initialize.h initialize.cpp l_elas.h l_elas.cpp lhsa.h lhsa.cpp diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h index 03f18c91..e83f21e1 100644 --- a/Code/Source/svFSI/ComMod.h +++ b/Code/Source/svFSI/ComMod.h @@ -1557,6 +1557,13 @@ class ComMod { Array pSn; Vector pSa; + // Variables for perfusion simulations + Vector perfusion_pressure_source; + Vector perfusion_pressure_sink; + Vector perfusion_beta0; + Vector perfusion_beta1; + double permeability = 0.0; + // Temporary storage for initializing state variables Vector Pinit; Array Vinit; diff --git a/Code/Source/svFSI/Parameters.cpp b/Code/Source/svFSI/Parameters.cpp index 78454a9b..df611081 100644 --- a/Code/Source/svFSI/Parameters.cpp +++ b/Code/Source/svFSI/Parameters.cpp @@ -1373,6 +1373,7 @@ DomainParameters::DomainParameters() set_parameter("Solid_density", 0.5, !required, solid_density); set_parameter("Solid_viscosity", 0.9, !required, solid_viscosity); set_parameter("Source_term", 0.0, !required, source_term); + set_parameter("Permeability", 0.0, !required, permeability); set_parameter("Time_step_for_integration", 0.0, !required, time_step_for_integration); } diff --git a/Code/Source/svFSI/Parameters.h b/Code/Source/svFSI/Parameters.h index 737d9ab8..e53ce695 100644 --- a/Code/Source/svFSI/Parameters.h +++ b/Code/Source/svFSI/Parameters.h @@ -1074,6 +1074,7 @@ class DomainParameters : public ParameterLists Parameter solid_density; Parameter solid_viscosity; Parameter source_term; + Parameter permeability; Parameter time_step_for_integration; }; diff --git a/Code/Source/svFSI/all_fun.cpp b/Code/Source/svFSI/all_fun.cpp index df230deb..6abc171e 100644 --- a/Code/Source/svFSI/all_fun.cpp +++ b/Code/Source/svFSI/all_fun.cpp @@ -896,6 +896,48 @@ double jacobian(ComMod& com_mod, const int nDim, const int eNoN, const Array //kmenon_perfusion + local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Vector& U) + { + Vector local_vector; + + if (com_mod.ltg.size() == 0) { + throw std::runtime_error("ltg is not set yet"); + } + + if (cm.mas(cm_mod)) { + if (U.size() != com_mod.gtnNo) { + throw std::runtime_error("local_rs is only specified for vector with size gtnNo"); + } + } + + if (cm.seq()) { + local_vector.resize(com_mod.gtnNo); + local_vector = U; + return local_vector; + } + + local_vector.resize(com_mod.tnNo); + Vector tmpU(com_mod.gtnNo); + + if (cm.mas(cm_mod)) { + tmpU = U; + } + + cm.bcast(cm_mod, tmpU); + + for (int a = 0; a < com_mod.tnNo; a++) { + int Ac = com_mod.ltg[a]; + local_vector[a] = tmpU[Ac]; + } + + return local_vector; + } + //------- // local //------- diff --git a/Code/Source/svFSI/consts.cpp b/Code/Source/svFSI/consts.cpp index d720ab60..f1380960 100644 --- a/Code/Source/svFSI/consts.cpp +++ b/Code/Source/svFSI/consts.cpp @@ -172,6 +172,7 @@ const std::map equation_name_to_type = { {"heatS", EquationType::phys_heatS}, {"laplace", EquationType::phys_heatS}, {"poisson", EquationType::phys_heatS}, + {"perfusion", EquationType::phys_darcy}, {"stokes", EquationType::phys_stokes}, diff --git a/Code/Source/svFSI/consts.h b/Code/Source/svFSI/consts.h index 259eb193..17cce90e 100644 --- a/Code/Source/svFSI/consts.h +++ b/Code/Source/svFSI/consts.h @@ -296,7 +296,8 @@ enum class EquationType phys_CMM = 209, phys_CEP = 210, phys_ustruct = 211, // Nonlinear elastodynamics using mixed VMS-stabilized formulation - phys_stokes = 212 + phys_stokes = 212, + phys_darcy = 213 // Darcy flow for perfusion problems }; constexpr auto Equation_CMM = EquationType::phys_CMM; @@ -311,6 +312,7 @@ constexpr auto Equation_shell = EquationType::phys_shell; constexpr auto Equation_stokes = EquationType::phys_stokes; constexpr auto Equation_struct = EquationType::phys_struct; constexpr auto Equation_ustruct = EquationType::phys_ustruct; +constexpr auto Equation_darcy = EquationType::phys_darcy; extern const std::map equation_name_to_type; @@ -359,6 +361,7 @@ enum class OutputType outGrp_fS = 523, outGrp_C = 524, outGrp_I1 = 525, + outGrp_MBF = 526, out_velocity = 599, out_pressure = 598, @@ -387,7 +390,8 @@ enum class OutputType out_viscosity = 575, out_fibStrn = 574, out_CGstrain = 573, - out_CGInv1 = 572 + out_CGInv1 = 572, + out_MBF = 571 }; //--------------------- @@ -412,7 +416,8 @@ enum class PhysicalProperyType damping = 12, shell_thickness = 13, ctau_M = 14, // stabilization coeffs. for USTRUCT (momentum, continuity) - ctau_C = 15 + ctau_C = 15, + permeability = 16 }; //-------------------- diff --git a/Code/Source/svFSI/darcy.cpp b/Code/Source/svFSI/darcy.cpp index aeac8360..919e157a 100644 --- a/Code/Source/svFSI/darcy.cpp +++ b/Code/Source/svFSI/darcy.cpp @@ -1,3 +1,287 @@ -// -// Created by Zack on 9/19/2023. -// +/* + This code implements the Darcy Equation for 2D and 3D + problems in perfusion of porus media. + ------------------------------------------------------------- + Assumptions: + - Homogeneous Permeability + - Homogeneous Density + - Isotropic Permeability + - Assumptions of Stokes Flow + - Steady-State + ------------------------------------------------------------- + Strong form of the Single-Compartment Darcy equation: + u = -K∇(P) + ∇⋅u = β0(P_source - P) - β1(P - P_sink) + where: + u -> Volume flux vector + K -> Permeability tensor + P -> Pressure + ------------------------------------------------------------- + Weak form of the Single-Compartment Darcy equation: + -∫(∇q∇P)dΩ - λ∫qPdΩ = ∫qFdΩ - ∫q∇P⋅nvdΓ + where: + q -> Test function + λ -> (β0 + β1)/K + F -> -(β0(P_source) + β1(P_sink))/K + n -> Normal vector to the boundary + +*/ + +#include "darcy.h" + +#include "all_fun.h" +#include "lhsa.h" +#include "mat_fun.h" +#include "nn.h" +#include "utils.h" + +#ifdef WITH_TRILINOS +#include "trilinos_linear_solver.h" +#endif + +namespace darcy { + //--------- + // b_darcy + //--------- + // + void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR) + { + for (int a = 0; a < eNoN; a++) { + lR(0,a) = lR(0,a) + w * N(a) * h; + } + } + + void construct_darcy(ComMod& com_mod, const mshType& lM, const Array& Ag, const Array& Yg) { +#define n_debug_construct_darcy +#ifdef debug_construct_darcy + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); +#endif + + using namespace consts; + + const int nsd = com_mod.nsd; + const int tDof = com_mod.tDof; + const int dof = com_mod.dof; + const int cEq = com_mod.cEq; + const auto &eq = com_mod.eq[cEq]; + auto &cDmn = com_mod.cDmn; + + + int eNoN = lM.eNoN; +#ifdef debug_construct_darcy + dmsg << "cEq: " << cEq; + dmsg << "cDmn: " << cDmn; +#endif + + Vector ptr(eNoN); + Vector N(eNoN); + Array xl(nsd, eNoN), al(tDof, eNoN), yl(tDof, eNoN); + Array Nx(nsd, eNoN), lR(dof, eNoN); + Array3 lK(dof * dof, eNoN, eNoN); + Array ksix(nsd, nsd); + Vector local_source(eNoN), local_sink(eNoN), local_b0(eNoN), local_b1(eNoN); + + for (int e = 0; e < lM.nEl; e++) { + cDmn = all_fun::domain(com_mod, lM, cEq, e); + auto cPhys = eq.dmn[cDmn].phys; + if (cPhys != EquationType::phys_darcy) { + std::cout << "cPhys: " << cPhys << std::endl; + continue; + } + + // Update shape function for NURBS + if (lM.eType == ElementType::NRB) { + //CALL NRMNNX(lm, e) + } + + // Create local copies + for (int a = 0; a < eNoN; a++) { + int Ac = lM.IEN(a, e); + ptr(a) = Ac; + + for (int i = 0; i < nsd; i++) { + xl(i, a) = com_mod.x(i, Ac); + } + + for (int i = 0; i < tDof; i++) { + al(i, a) = Ag(i, Ac); + yl(i, a) = Yg(i, Ac); + } + + local_source(a) = com_mod.perfusion_pressure_source(Ac); + local_sink(a) = com_mod.perfusion_pressure_sink(Ac); + local_b0(a) = com_mod.perfusion_beta0(Ac); + local_b1(a) = com_mod.perfusion_beta1(Ac); + } + + // Gauss integration + + lR = 0.0; + lK = 0.0; + double Jac{0.0}; + + for (int g = 0; g < lM.nG; g++) { + if (g == 0 || !lM.lShpF) { + auto Nx_g = lM.Nx.slice(g); + nn::gnn(eNoN, nsd, nsd, Nx_g, xl, Nx, Jac, ksix); + if (utils::is_zero(Jac)) { + throw std::runtime_error( + "[construct_darcy] Jacobian for element " + std::to_string(e) + " is < 0."); + } + } + + double w = lM.w(g) * Jac; + N = lM.N.col(g); + + if (nsd == 3) { + darcy_3d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); + } else if (nsd == 2) { + darcy_2d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK, local_source, local_sink, local_b0, local_b1); + } else { + throw std::runtime_error("[construct_darcy] nsd must be 2 or 3."); + } + } + } + + // Assembly +#ifdef WITH_TRILINOS + if (eq.assmTLS) { + trilinos_doassem_(const_cast(eNoN), const_cast(ptr.data()), lK.data(), lR.data()); + } +#endif + lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR); + } + + //---------- + // darcy_2d + //---------- + // + void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK, + Vector& local_source, Vector& local_sink, Vector& b0, Vector& b1) { +#define n_debug_darcy_2d +#ifdef debug_darcy_2d + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "w: " << w; +#endif + using namespace consts; + using namespace mat_fun; + + const int nsd = com_mod.nsd; + const int cEq = com_mod.cEq; + auto &eq = com_mod.eq[cEq]; + const int cDmn = com_mod.cDmn; + auto &dmn = eq.dmn[cDmn]; + const double dt = com_mod.dt; + const int i = eq.s; + + double k = dmn.prop.at(PhysicalProperyType::permeability); // should be a tensor + double source = dmn.prop.at(PhysicalProperyType::source_term); //is this just a double? + double rho = dmn.prop.at(PhysicalProperyType::solid_density); // not needed + + double T1 = eq.af * eq.gam * dt; + double amd = eq.am * rho/ T1; + double wl = w * T1; + +#ifdef debug_darcy_2d + dmsg; + dmsg << "nu: " << nu; + dmsg << "source: " << source; + dmsg << "rho: " << rho ; + dmsg << "T1: " << T1; + dmsg << "i: " << i; + dmsg << "wl: " << wl; +#endif + double Td = -source; + Vector Tx(nsd); + + + for (int a = 0; a < eNoN; a++) { + Td = Td + N(a)*al(i,a); + Tx(0) = Tx(0) + Nx(0,a)*yl(i,a); + Tx(1) = Tx(1) + Nx(1,a)*yl(i,a); + } + + + Td = Td * rho; + +#ifdef debug_darcy_2d + dmsg; + dmsg << "Td: " << Td; + dmsg << "Tx: " << Tx; +#endif + + for (int a = 0; a < eNoN; a++) { + lR(0,a) = lR(0,a) + w*(N(a)*Td + (Nx(0,a)*Tx(0) + Nx(1,a)*Tx(1))*k); + //std::cout << "lR(0, a): " << lR(0, a) << std::endl; + for (int b = 0; b < eNoN; b++) { + lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*amd + k*(Nx(0,a)*Nx(0,b) + Nx(1,a)*Nx(1,b))); + } + } + } + + //---------- + // darcy_3d + //---------- + // + void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK) { +#define n_debug_darcy_3d +#ifdef debug_darcy_3d + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "w: " << w; +#endif + using namespace consts; + using namespace mat_fun; + + const int nsd = com_mod.nsd; + const int cEq = com_mod.cEq; + auto &eq = com_mod.eq[cEq]; + const int cDmn = com_mod.cDmn; + auto &dmn = eq.dmn[cDmn]; + const double dt = com_mod.dt; + const int i = eq.s; + + double k = dmn.prop.at(PhysicalProperyType::permeability); // should be a tensor + double source = dmn.prop.at(PhysicalProperyType::source_term); //is this just a double? + double rho = dmn.prop.at(PhysicalProperyType::solid_density); + + double T1 = eq.af * eq.gam * dt; + double amd = eq.am * rho / T1; + double wl = w * T1; + +#ifdef debug_darcy_3d + dmsg; + dmsg << "k: " << k; + dmsg << "source: " << source; + dmsg << "rho: " << rho ; + dmsg << "T1: " << T1; + dmsg << "i: " << i; + dmsg << "wl: " << wl; +#endif + double Td = -source; + Vector Tx(nsd); + + for (int a = 0; a < eNoN; a++){ + Td = Td + N(a) * al(i,a); + Tx(0) = Tx(0) + Nx(0,a) * yl(i,a); + Tx(1) = Tx(1) + Nx(1,a) * yl(i,a); + Tx(2) = Tx(2) + Nx(2,a) * yl(i,a); + } + + Td = Td * rho; + + for (int a = 0; a < eNoN; a++){ + lR(0,a) = lR(0, a) + w*(N(a)*Td + + k*(Nx(0,a)*Tx(0) + Nx(1,a)*Tx(1) + Nx(2,a)*Tx(2))); + for (int b = 0; b < eNoN; b++){ + lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*amd + + k*(Nx(0,a)*Nx(0,b) + Nx(1,a)*Nx(1,b) + + Nx(2,a)*Nx(2,b))); + } + } + } +} \ No newline at end of file diff --git a/Code/Source/svFSI/darcy.h b/Code/Source/svFSI/darcy.h index cd4f9ba9..914a3193 100644 --- a/Code/Source/svFSI/darcy.h +++ b/Code/Source/svFSI/darcy.h @@ -1,8 +1,25 @@ // -// Created by Zack on 9/19/2023. +// This code implements the Darcy Equation for 2D and 3D +// problems in perfusion of porus media. // -#ifndef SV_TOP_DARCY_H -#define SV_TOP_DARCY_H +#ifndef DARCY_H +#define DARCY_H -#endif //SV_TOP_DARCY_H +#include "ComMod.h" + +namespace darcy { + + void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR); + + void construct_darcy(ComMod& com_mod, const mshType& lM, const Array& Ag, const Array& Dg); + + void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK, Vector& local_source, + Vector& local_sink, Vector& local_b0, Vector& local_b1); + + void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK); +} + +#endif //DARCY_H diff --git a/Code/Source/svFSI/eq_assem.cpp b/Code/Source/svFSI/eq_assem.cpp index 47b6fb7d..11ed3ae3 100644 --- a/Code/Source/svFSI/eq_assem.cpp +++ b/Code/Source/svFSI/eq_assem.cpp @@ -20,6 +20,7 @@ #include "stokes.h" #include "sv_struct.h" #include "ustruct.h" +#include "darcy.h" #include @@ -108,6 +109,10 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& heats::b_heats(com_mod, eNoN, w, N, h, lR); break; + case EquationType::phys_darcy: + darcy::b_darcy(com_mod, eNoN, w, N, h, lR); + break; + case EquationType::phys_heatF: heatf::b_heatf(com_mod, eNoN, w, N, y, h, nV, lR, lK); break; @@ -336,6 +341,10 @@ void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const heats::construct_heats(com_mod, lM, Ag, Yg); break; + case EquationType::phys_darcy: + darcy::construct_darcy(com_mod, lM, Ag, Yg); + break; + case EquationType::phys_lElas: l_elas::construct_l_elas(com_mod, lM, Ag, Dg); break; diff --git a/Code/Source/svFSI/read_files.cpp b/Code/Source/svFSI/read_files.cpp index fc8d9487..44076e9d 100644 --- a/Code/Source/svFSI/read_files.cpp +++ b/Code/Source/svFSI/read_files.cpp @@ -1156,6 +1156,9 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& case PhysicalProperyType::source_term: rtmp = domain_params->source_term.value(); break; + + case PhysicalProperyType::permeability: + rtmp = domain_params->permeability.value(); } lEq.dmn[iDmn].prop[prop] = rtmp; @@ -1591,7 +1594,40 @@ void read_files(Simulation* simulation, const std::string& file_name) if (!com_mod.mvMsh) { throw std::runtime_error("mesh equation can only be specified after FSI equation"); } - } + } + + if (eq.phys == EquationType::phys_darcy) { + com_mod.perfusion_pressure_source.resize(com_mod.gtnNo); + com_mod.perfusion_pressure_sink.resize(com_mod.gtnNo); + com_mod.perfusion_beta0.resize(com_mod.gtnNo); + com_mod.perfusion_beta1.resize(com_mod.gtnNo); + std::ifstream perfusion_pressure_source_file; + std::ifstream perfusion_pressure_sink_file; + std::ifstream perfusion_pressure_b0_file; + std::ifstream perfusion_pressure_b1_file; + perfusion_pressure_source_file.open("perfusionSource"); + perfusion_pressure_sink_file.open("perfusionSink"); + perfusion_pressure_b0_file.open("perfusionB0"); + perfusion_pressure_b1_file.open("perfusionB1"); + if (!perfusion_pressure_source_file.is_open()) { + throw std::runtime_error("Failed to open the perfusion pressure source file 'perfusionSource'."); + } + if (!perfusion_pressure_sink_file.is_open()) { + throw std::runtime_error("Failed to open the perfusion pressure source file 'perfusionSource'."); + } + if (!perfusion_pressure_b0_file.is_open()) { + throw std::runtime_error("Failed to open the perfusion pressure source file 'perfusionSource'."); + } + if (!perfusion_pressure_b1_file.is_open()) { + throw std::runtime_error("Failed to open the perfusion pressure source file 'perfusionSource'."); + } + for (int i = 0; i < com_mod.gtnNo; i++) { + perfusion_pressure_source_file >> com_mod.perfusion_pressure_source[i]; + perfusion_pressure_sink_file >> com_mod.perfusion_pressure_sink[i]; + perfusion_pressure_b0_file >> com_mod.perfusion_beta0[i]; + perfusion_pressure_b1_file >> com_mod.perfusion_beta1[i]; + } + } } #ifdef debug_read_files dmsg << "Done Read equations " << " "; diff --git a/Code/Source/svFSI/set_equation_dof.h b/Code/Source/svFSI/set_equation_dof.h index 71ec1ae0..43dc4e4c 100644 --- a/Code/Source/svFSI/set_equation_dof.h +++ b/Code/Source/svFSI/set_equation_dof.h @@ -11,6 +11,7 @@ std::map equation_dof_map = {EquationType::phys_fluid, std::make_tuple(nsd+1, "NS") }, {EquationType::phys_heatF, std::make_tuple(1, "HF") }, {EquationType::phys_heatS, std::make_tuple(1, "HS") }, + {EquationType::phys_darcy, std::make_tuple(1, "HS") }, {EquationType::phys_lElas, std::make_tuple(nsd, "LE") }, {EquationType::phys_struct, std::make_tuple(nsd, "ST") }, {EquationType::phys_ustruct, std::make_tuple(nsd+1, "ST") }, diff --git a/Code/Source/svFSI/set_equation_props.h b/Code/Source/svFSI/set_equation_props.h index 17f8dbfb..a7437406 100644 --- a/Code/Source/svFSI/set_equation_props.h +++ b/Code/Source/svFSI/set_equation_props.h @@ -267,7 +267,31 @@ SetEquationPropertiesMapType set_equation_props = { } }, //---------------------------// -// phys_FSI // +// phys_darcy // +//---------------------------// +{consts::EquationType::phys_darcy, [](Simulation* simulation, EquationParameters* eq_params, eqType& lEq, EquationProps& propL, + EquationOutputs& outPuts, EquationNdop& nDOP) -> void +{ + using namespace consts; + auto& com_mod = simulation->get_com_mod(); + lEq.phys = consts::EquationType::phys_darcy; + + propL[0][0] = PhysicalProperyType::permeability; + propL[1][0] = PhysicalProperyType::source_term; + propL[2][0] = PhysicalProperyType::solid_density; + + read_domain(simulation, eq_params, lEq, propL); + + nDOP = {2,1,1,0}; + outPuts = {OutputType::out_temperature, OutputType::out_heatFlux}; + + // Set solver parameters. + read_ls(simulation, eq_params, SolverType::lSolver_CG, lEq); +} }, + + +//---------------------------// +// phys_FSI //const //---------------------------// // {consts::EquationType::phys_FSI, [](Simulation* simulation, EquationParameters* eq_params, eqType& lEq, EquationProps& propL, diff --git a/Code/Source/svFSI/set_output_props.h b/Code/Source/svFSI/set_output_props.h index 1272c47f..19426de6 100644 --- a/Code/Source/svFSI/set_output_props.h +++ b/Code/Source/svFSI/set_output_props.h @@ -54,6 +54,7 @@ std::map output_props_map = {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") } + {OutputType::out_WSS, std::make_tuple(OutputType::outGrp_WSS, 0, maxNSD, "WSS")}, + {OutputType::out_MBF, std::make_tuple(OutputType::outGrp_MBF, 0, 1, "MBF")} }; From 4e7d7486eb7d578b6869c284509f6e4a72a1215a Mon Sep 17 00:00:00 2001 From: zasexton Date: Thu, 5 Oct 2023 13:32:26 -0700 Subject: [PATCH 03/15] commenting files to find mesh indexing error --- Code/Source/svFSI/all_fun.cpp | 26 ++++++++++- Code/Source/svFSI/baf_ini.cpp | 3 +- Code/Source/svFSI/darcy.cpp | 20 ++++---- Code/Source/svFSI/heats.cpp | 9 +++- Code/Source/svFSI/initialize.cpp | 2 +- Code/Source/svFSI/load_msh.cpp | 3 +- Code/Source/svFSI/nn.cpp | 58 ++++++++++++++++++++++-- Code/Source/svFSI/read_msh.cpp | 57 ++++++++++++++--------- Code/Source/svFSI/vtk_xml.cpp | 8 ++-- Code/Source/svFSI/vtk_xml_parser.cpp | 68 ++++++++++++++++++++++++---- 10 files changed, 199 insertions(+), 55 deletions(-) diff --git a/Code/Source/svFSI/all_fun.cpp b/Code/Source/svFSI/all_fun.cpp index 6abc171e..3aab10e9 100644 --- a/Code/Source/svFSI/all_fun.cpp +++ b/Code/Source/svFSI/all_fun.cpp @@ -536,7 +536,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co dmsg << "lFa.name: " << lFa.name; dmsg << "lFa.eType: " << lFa.eType; #endif - + std::cout << "Integrate with Flag" << std::endl; bool flag = pFlag; int nsd = com_mod.nsd; int insd = nsd - 1; @@ -621,6 +621,29 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co #endif double result = 0.0; + std::cout << "lFa.gN" << std::endl; + std::cout << lFa.gN << std::endl; + std::cout << "lFa.gE" << std::endl; + std::cout << lFa.gE << std::endl; + std::cout << "lFa.IEN" << std::endl; + for (int i = 0; i < lFa.nEl; i++) { + std::cout << i << " - "; + for (int j = 0; j < lFa.eNoN; j++) { + std::cout << lFa.IEN(j, i) << " "; + } + std::cout << std::endl; + } + std::cout << std::endl; + std::cout << "msh.IEN" << std::endl; + for (int i = 0; i < com_mod.msh[lFa.iM].nEl; i++) { + std::cout << i << " - "; + for (int j = 0; j < com_mod.msh[lFa.iM].eNoN; j++) { + std::cout << com_mod.msh[lFa.iM].IEN(j,i) << " "; + } + std::cout << std::endl; + } + + for (int e = 0; e < lFa.nEl; e++) { // [TODO:DaveP] not implemented. if (lFa.eType == ElementType::NRB) { @@ -713,7 +736,6 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co } double result = 0.0; - for (int e = 0; e < lFa.nEl; e++) { //dmsg << "----- e " << e+1 << " -----"; // Updating the shape functions, if this is a NURB diff --git a/Code/Source/svFSI/baf_ini.cpp b/Code/Source/svFSI/baf_ini.cpp index da7eae15..aae3490c 100644 --- a/Code/Source/svFSI/baf_ini.cpp +++ b/Code/Source/svFSI/baf_ini.cpp @@ -425,11 +425,12 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa) // Vector sA(com_mod.tnNo); sA = 1.0; + std::cout << "Start Area Integration" << std::endl; double area = all_fun::integ(com_mod, cm_mod, lFa, sA); + std::cout << "Face '" << lFa.name << "' area: " << area << std::endl; #ifdef debug_face_ini dmsg << "Face '" << lFa.name << "' area: " << area; #endif - if (utils::is_zero(area)) { if (cm.mas(cm_mod)) { throw std::runtime_error("Face '" + lFa.name + "' has zero area."); diff --git a/Code/Source/svFSI/darcy.cpp b/Code/Source/svFSI/darcy.cpp index 919e157a..13fcbc4d 100644 --- a/Code/Source/svFSI/darcy.cpp +++ b/Code/Source/svFSI/darcy.cpp @@ -142,17 +142,16 @@ namespace darcy { throw std::runtime_error("[construct_darcy] nsd must be 2 or 3."); } } - } - // Assembly + // Assembly #ifdef WITH_TRILINOS - if (eq.assmTLS) { - trilinos_doassem_(const_cast(eNoN), const_cast(ptr.data()), lK.data(), lR.data()); - } + if (eq.assmTLS) { + trilinos_doassem_(const_cast(eNoN), const_cast(ptr.data()), lK.data(), lR.data()); + } #endif - lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR); + lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR); + } } - //---------- // darcy_2d //---------- @@ -214,10 +213,11 @@ namespace darcy { #endif for (int a = 0; a < eNoN; a++) { - lR(0,a) = lR(0,a) + w*(N(a)*Td + (Nx(0,a)*Tx(0) + Nx(1,a)*Tx(1))*k); - //std::cout << "lR(0, a): " << lR(0, a) << std::endl; + lR(0,a) = lR(0,a) + w*(N(a)*Td + + (Nx(0,a)*Tx(0) + Nx(1,a)*Tx(1))*k); for (int b = 0; b < eNoN; b++) { - lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*amd + k*(Nx(0,a)*Nx(0,b) + Nx(1,a)*Nx(1,b))); + lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*amd + + k*(Nx(0,a)*Nx(0,b) + Nx(1,a)*Nx(1,b))); } } } diff --git a/Code/Source/svFSI/heats.cpp b/Code/Source/svFSI/heats.cpp index 2cb24853..d3b57bc4 100644 --- a/Code/Source/svFSI/heats.cpp +++ b/Code/Source/svFSI/heats.cpp @@ -16,7 +16,12 @@ namespace heats { //--------- // b_heats //--------- -// + +// eNoN - number of nodes in element +// lR - local residual +// w - gauss point weight +// N - shape functions +// h - heat flux at boundary void b_heats(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR) { for (int a = 0; a < eNoN; a++) { @@ -184,7 +189,7 @@ void heats_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& timeP) // std::tie(eq.dof, eq.sym) = equation_dof_map.at(eq.phys); - if (std::set{Equation_fluid, Equation_heatF, Equation_heatS, Equation_CEP, Equation_stokes}.count(eq.phys) == 0) { + if (std::set{Equation_fluid, Equation_heatF, Equation_heatS, Equation_CEP, Equation_stokes, Equation_darcy}.count(eq.phys) == 0) { dFlag = true; } diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index 3e1b0de8..6001848a 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -161,7 +161,7 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p // If node IDs were not read then create them. if (face.gN.size() == 0) { read_msh_ns::calc_nbc(mesh, face); - + std::cout << "creating node ids!" << std::endl; // Reset the connecttivity with the new node IDs? for (int e = 0; e < face.nEl; e++) { for (int a = 0; a < face.eNoN; a++) { @@ -174,6 +174,7 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p } nn::select_eleb(simulation, mesh, face); + std::cout << "Done face.nEl: " << face.nEl << std::endl; } } diff --git a/Code/Source/svFSI/nn.cpp b/Code/Source/svFSI/nn.cpp index 8b1fd56b..ffc7772e 100644 --- a/Code/Source/svFSI/nn.cpp +++ b/Code/Source/svFSI/nn.cpp @@ -561,6 +561,7 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, auto& msh = com_mod.msh[iM]; int eNoN = msh.eNoN; + #ifdef debug_gnnb dmsg << "iM: " << iM+1; dmsg << "Ec: " << Ec+1; @@ -582,18 +583,32 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, for (int a = 0; a < eNoNb; a++) { int Ac = lFa.IEN(a,e); int b = 0; - + bool match = false; for (int ib = 0; ib < eNoN; ib++) { b = ib; if (setIt[ib]) { int Bc = msh.IEN(ib,Ec); if (Bc == Ac) { + //std::cout << "Match! " << Ac << std::endl; + match = true; break; + } else { + //std::cout << "No match! " << "Ac: " << Ac << " Bc: " << Bc << std::endl; + match = false; } } } - if (b > eNoN) { + if (!match) { + std::cout << "Face element: " << e << std::endl; + std::cout << "Global Node : " << Ac << std::endl; + std::cout << "Mesh element: " << Ec << std::endl; + + //std::cout << "Global Nodes in Mesh Element: IEN(" << ib << "," << Ec << ")" << std::endl; + for (int ib = 0; ib < eNoN; ib++) { + std::cout << msh.IEN(ib,Ec) << " "; + } + std::cout << std::endl; throw std::runtime_error("could not find matching face nodes"); } @@ -609,14 +624,21 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, a = a + 1; } } - + //std::cout << "ptr("<< Ec << ")" << std::endl; + for (int i = 0; i < eNoN; i++) { + //std::cout << ptr(i) << " "; + } + //std::cout << std::endl; // Correcting the position vector if mesh is moving // + //std::cout << "lX("<< Ec << ")" << std::endl; for (int a = 0; a < eNoN; a++) { int Ac = msh.IEN(a,Ec); for (int i = 0; i < lX.nrows(); i++) { lX(i,a) = com_mod.x(i,Ac); + //std::cout << lX(i,a) << " "; } + //std::cout << std::endl; if (com_mod.mvMsh) { for (int i = 0; i < lX.nrows(); i++) { @@ -647,6 +669,7 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, } } + auto v = utils::cross(xXi); for (int i = 0; i < nsd; i++) { v(i) = v(i) / sqrt(utils::norm(v)); @@ -690,16 +713,43 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, Array xXi(nsd,insd); + //std::cout << "Nx: " << std::endl; + for (int i = 0; i v1, v2, v3; int num_elems = mesh.gnEl; int num_reorder = 0; @@ -1149,6 +1149,7 @@ void read_msh(Simulation* simulation) // // READMSH - CALL READSV(lPM, msh(iM)) // + std::cout << "Number of Mesh Elements: " << com_mod.nMsh << std::endl; for (int iM = 0; iM < com_mod.nMsh; iM++) { auto param = simulation->parameters.mesh_parameters[iM]; auto& mesh = com_mod.msh[iM]; @@ -1185,9 +1186,11 @@ void read_msh(Simulation* simulation) #ifdef debug_read_msh dmsg << "Check for unique face names " << " ...."; #endif + std::cout << "Number of Faces: " << mesh.nFa << std::endl; for (int iFa = 0; iFa < mesh.nFa; iFa++) { mesh.fa[iFa].iM = iM; auto ctmp = mesh.fa[iFa].name; + std::cout << "Face Name: " << ctmp << std::endl; for (int i = 0; i < iM; i++) { for (int j = 0; j < com_mod.msh[i].nFa; j++) { if ((ctmp == com_mod.msh[i].fa[j].name) && ((i != iM || j != iFa))) { @@ -1196,7 +1199,7 @@ void read_msh(Simulation* simulation) } } } - + std::cout << "Done Face Name Checking." << std::endl; // Global total number of nodes. com_mod.gtnNo += mesh.gnNo; #ifdef debug_read_msh @@ -1226,7 +1229,7 @@ void read_msh(Simulation* simulation) } mesh.x.clear(); } - + std::cout << "Done creating global nodal coordinate array." << std::endl; com_mod.x = gX; // Check for shell elements. @@ -1273,6 +1276,7 @@ void read_msh(Simulation* simulation) dmsg << "Allocate new mesh nodes or something " << " ..."; dmsg << "com_mod.gtnNo: " << com_mod.gtnNo; #endif + std::cout << "gX size: " << com_mod.nsd << ", " << com_mod.gtnNo << std::endl; gX.resize(com_mod.nsd,com_mod.gtnNo); gX = com_mod.x; @@ -1289,7 +1293,8 @@ void read_msh(Simulation* simulation) msh.gN.resize(msh.gnNo); msh.gN = -1; } - } + } + std::cout << "Done allocating new mesh nodes or something." << std::endl; // Examining the existance of projection faces and setting %gN. // Reseting gtnNo and recounting nodes that are not duplicated @@ -1316,6 +1321,7 @@ void read_msh(Simulation* simulation) } } } + std::cout << "Done examining the existance of projection faces and setting %gN." << std::endl; #ifdef debug_read_msh dmsg << "Allocate com_mod.x " << " ..."; @@ -1345,7 +1351,7 @@ void read_msh(Simulation* simulation) com_mod.msh[iM].lN[Ac] = a; } } - + std::cout << "Done temporarily allocating msh%lN array." << std::endl; // Re-arranging x and finding the size of the entire domain // First rearrange 2D/3D mesh and then, 1D fiber mesh // @@ -1388,13 +1394,12 @@ void read_msh(Simulation* simulation) } b = b + 1; - } - + } #ifdef debug_read_msh dmsg << "b: " << b+1; #endif - } - + } + std::cout << "Done rearranging x." << std::endl; b = 0; for (int iM = 0; iM < com_mod.nMsh; iM++) { @@ -1418,24 +1423,23 @@ void read_msh(Simulation* simulation) } for (int i = 0; i < com_mod.nsd; i++) { - com_mod.x(i,Ac) = gX(i,b); + //com_mod.x(i,Ac) = gX(i,b); if (maxX[i] < com_mod.x(i,Ac)) { - maxX[i] = com_mod.x(i,Ac); + //maxX[i] = com_mod.x(i,Ac); } if (minX[i] > com_mod.x(i,Ac)) { - minX[i] = com_mod.x(i,Ac); + //minX[i] = com_mod.x(i,Ac); } } b = b + 1; } - #ifdef debug_read_msh dmsg << "b: " << b+1; #endif } - + std::cout << "Done rearranging x." << std::endl; #ifdef debug_read_msh dmsg << "minX: " << minX; dmsg << "maxX: " << maxX; @@ -1447,24 +1451,32 @@ void read_msh(Simulation* simulation) dmsg << "Rearrange " << " fa ... "; #endif + + std::cout << "Number of Meshes: " << com_mod.nMsh << std::endl; for (int iM = 0; iM < com_mod.nMsh; iM++) { for (int iFa = 0; iFa < com_mod.msh[iM].nFa; iFa++) { for (int a = 0; a < com_mod.msh[iM].fa[iFa].nNo; a++) { - int Ac = com_mod.msh[iM].fa[iFa].gN[a]; - Ac = com_mod.msh[iM].gN[Ac]; - com_mod.msh[iM].fa[iFa].gN[a] = Ac; + //int Ac = com_mod.msh[iM].fa[iFa].gN[a]; + //std::cout << "Ac: " << Ac << std::endl; + //Ac = com_mod.msh[iM].gN[Ac]; + //com_mod.msh[iM].fa[iFa].gN[a] = Ac; } - + std::cout << "Done rearranging fa inner most loop 1." << std::endl; for (int e = 0; e < com_mod.msh[iM].fa[iFa].nEl; e++) { for (int a = 0; a < com_mod.msh[iM].fa[iFa].eNoN; a++) { - int Ac = com_mod.msh[iM].fa[iFa].IEN(a,e); - Ac = com_mod.msh[iM].gN[Ac]; - com_mod.msh[iM].fa[iFa].IEN(a,e) = Ac; + //int Ac = com_mod.msh[iM].fa[iFa].IEN(a,e); + //Ac = com_mod.msh[iM].gN[Ac]; + //std::cout << "IEN(" << a << "," << e << "): " << Ac << std::endl; + //com_mod.msh[iM].fa[iFa].IEN(a,e) = Ac; + //std::cout << "IEN(" << a << "," << e << "): " << Ac << " Done." << std::endl; } - } + } + std::cout << "Done rearranging fa inner most loop 2." << std::endl; } } + std::cout << "Done rearranging fa." << std::endl; + if (com_mod.resetSim) { auto& rmsh = com_mod.rmsh; int gtnNo = com_mod.gtnNo; @@ -1502,6 +1514,7 @@ void read_msh(Simulation* simulation) } } + std::cout << "Done resetSim." << std::endl; // Setting dmnId parameter here, if there is at least one mesh that // has defined eId. // diff --git a/Code/Source/svFSI/vtk_xml.cpp b/Code/Source/svFSI/vtk_xml.cpp index 265d6e54..dc583278 100644 --- a/Code/Source/svFSI/vtk_xml.cpp +++ b/Code/Source/svFSI/vtk_xml.cpp @@ -426,9 +426,10 @@ void read_vtp(const std::string& file_name, faceType& face) } else { for (int e = 0; e < face.nEl; e++) { for (int a = 0; a < face.eNoN; a++) { - int Ac = face.IEN(a,e); - Ac = face.gN(Ac); - face.IEN(a,e) = Ac; + //int Ac = face.IEN(a,e); + //Ac = face.gN(Ac); + //face.IEN(a,e) = Ac; + //std::cout << "[read_vtp] IEN(" << a << "," << e << "): " << face.IEN(a,e) << std::endl; } } } @@ -451,6 +452,7 @@ void read_vtp(const std::string& file_name, faceType& face) //std::cout << std::endl; } } + } //---------------- diff --git a/Code/Source/svFSI/vtk_xml_parser.cpp b/Code/Source/svFSI/vtk_xml_parser.cpp index 91947d59..b431b000 100644 --- a/Code/Source/svFSI/vtk_xml_parser.cpp +++ b/Code/Source/svFSI/vtk_xml_parser.cpp @@ -83,14 +83,22 @@ void store_element_conn(vtkSmartPointer vtk_polydata, faceType& fac face.eNoN = np_elem; face.IEN = Array(np_elem, num_elems); + auto global_node_ar = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); for (int i = 0; i < num_elems; i++) { vtk_polydata->GetCell(i, cell); + // Map the Extracted Local Cell Point Ids to the Global Volume Mesh Node Ids + //std::cout << "Begin Safe Down Cast from the Data Storage" << std::endl; + //auto global_node_ar = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); auto num_cell_pts = cell->GetNumberOfPoints(); + //std::cout << "Success on DownCast Operation" << std::endl; for (int j = 0; j < num_cell_pts; j++) { - auto id = cell->PointIds->GetId(j); + auto id = global_node_ar->GetValue(cell->PointIds->GetId(j)); + std::cout << "Id type: " << static_cast(id) <<" id: " << id << std::endl; face.IEN(j,i) = id; + //std::cout << "face_" < vtk_polydata, mshType& mesh) @@ -115,12 +123,15 @@ void store_element_conn(vtkSmartPointer vtk_polydata, mshType& mesh std::cout << "[store_element_conn(polydata,mesh)] np_elem: " << np_elem << std::endl; #endif + auto global_node_ar = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); for (int i = 0; i < num_elems; i++) { - //int elem_id = elem_ids->GetValue(i); + // int elem_id = elem_ids->GetValue(i); + // Node id's within the total mesh are the same since interior elements are not neglected vtk_polydata->GetCell(i, cell); auto num_cell_pts = cell->GetNumberOfPoints(); for (int j = 0; j < num_cell_pts; j++) { - auto id = cell->PointIds->GetId(j); + //auto id = cell->PointIds->GetId(j); + auto id = global_node_ar->GetValue(cell->PointIds->GetId(j)); mesh.gIEN(j,i) = id; } } @@ -144,6 +155,7 @@ void store_element_conn(vtkSmartPointer vtk_ugrid, mshType& std::cout << "[store_element_conn(ugrid)] " << std::endl; std::cout << "[store_element_conn(ugrid)] ========== store_element_conn(ugrid) =========" << std::endl; #endif + std::cout << "[store_element_conn(polydata,mesh)] connecting mesh" << std::endl; vtkSmartPointer cell_types = vtk_ugrid->GetCellTypesArray(); auto num_elems = vtk_ugrid->GetNumberOfCells(); #ifdef debug_store_element_conn @@ -231,7 +243,7 @@ void store_element_conn(vtkSmartPointer vtk_ugrid, mshType& std::cout << "[store_element_conn] Number of elements: " << num_elems << std::endl; std::cout << "[store_element_conn] Number of nodes per element: " << np_elem << std::endl; #endif - + auto global_node_ar = vtkIntArray::SafeDownCast(vtk_ugrid->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); auto cell = vtkGenericCell::New(); for (int i = 0; i < num_elems; i++) { vtk_ugrid->GetCell(i, cell); @@ -242,6 +254,9 @@ void store_element_conn(vtkSmartPointer vtk_ugrid, mshType& } for (int j = 0; j < num_cell_pts; j++) { auto id = cell->PointIds->GetId(j); + //auto id2 = global_node_ar->GetValue(id); + //std::cout << "mesh.gIEN("<< j << "," << i << ") = " << id << std::endl; + // std::cout << "mesh.gIEN("<< j << "," << i << ") = " << id2 << std::endl; //mesh.gIEN(i,j) = id; mesh.gIEN(j,i) = id; } @@ -288,7 +303,8 @@ void store_element_ids(vtkSmartPointer vtk_polydata, faceType& face // [NOTE] It is not clear how these IDs are used but if they // index into arrays or vectors then they need to be offset by -1. for (int i = 0; i < num_elem_ids; i++) { - face.gE(i) = elem_ids->GetValue(i) - 1; + face.gE(i) = elem_ids->GetValue(i); // - 1; assume that the face meshes are zero indexed + //std::cout << "[store_element_ids] face.gE(" << i << "): " << face.gE(i) << std::endl; } } @@ -376,14 +392,24 @@ void store_nodal_ids(vtkSmartPointer vtk_ugrid, mshType& me void store_nodal_ids(vtkSmartPointer vtk_polydata, faceType& face) { vtkIdType num_nodes = vtk_polydata->GetNumberOfPoints(); + vtkIdType num_pt_arrays = vtk_polydata->GetPointData()->GetNumberOfArrays(); + vtkIdType num_cell_arrays = vtk_polydata->GetCellData()->GetNumberOfArrays(); auto node_ids = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); + //for (vtkIdType i = 0; i < num_pt_arrays; i++) { + // if (vtk_polydata->GetPointData()->GetArrayName(i) == NODE_IDS_NAME){ + // //std::cout << "GlobalNodeID Array Number: " << i << std::endl; + // node_ids = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(i)); + // //std::cout << node_ids << std::endl; + // } + //} if (node_ids == nullptr) { return; } face.gN = Vector(num_nodes); for (int i = 0; i < num_nodes; i++) { + //std::cout << "[store_nodal_ids] node_ids->GetValue(" << i << "): " << std::endl; // [NOTE] It seems that face node IDs are 1 based. - face.gN(i) = node_ids->GetValue(i) - 1; + face.gN(i) = node_ids->GetValue(i); //- 1; //std::cout << "[store_nodal_ids] face.gN(" << i << "): " << face.gN(i) << std::endl; } } @@ -398,7 +424,7 @@ void store_nodal_ids(vtkSmartPointer vtk_polydata, mshType& mesh) mesh.gN = Vector(num_nodes); for (int i = 0; i < num_nodes; i++) { // [NOTE] It seems that face node IDs are 1 based. - mesh.gN(i) = node_ids->GetValue(i) - 1; + mesh.gN(i) = node_ids->GetValue(i); // - 1; //std::cout << "[store_nodal_ids] face.gN(" << i << "): " << face.gN(i) << std::endl; } } @@ -491,8 +517,10 @@ void load_vtp(const std::string& file_name, faceType& face) std::cout << "[load_vtp] ===== vtk_xml_parser.cpp::load_vtp ===== " << std::endl; #endif auto reader = vtkSmartPointer::New(); + std::cout << "[load_vtp] Reading VTK file '" << file_name << "' ... " << std::endl; reader->SetFileName(file_name.c_str()); reader->Update(); + std::cout << "[load_vtp] Done. " << std::endl; vtkSmartPointer vtk_polydata = reader->GetOutput(); vtkIdType num_nodes = vtk_polydata->GetNumberOfPoints(); @@ -501,6 +529,8 @@ void load_vtp(const std::string& file_name, faceType& face) } vtkIdType num_elems = vtk_polydata->GetNumberOfCells(); + std::cout << "[load_vtp] Number of nodes: " << num_nodes << std::endl; + std::cout << "[load_vtp] Number of elements: " << num_elems << std::endl; #ifdef debug_load_vtp std::cout << "[load_vtp] Number of nodes: " << num_nodes << std::endl; std::cout << "[load_vtp] Number of elements: " << num_elems << std::endl; @@ -515,10 +545,27 @@ void load_vtp(const std::string& file_name, faceType& face) // Store element connectivity. store_element_conn(vtk_polydata, face); - + //for (int e = 0; e < face.nEl; e++) { + // for (int a = 0; a < face.eNoN; a++) { + //int Ac = face.IEN(a,e); + //Ac = face.gN(Ac); + //face.IEN(a,e) = Ac; + // std::cout << "[read_vtp] IEN(" << a << "," << e << "): " << face.IEN(a,e) << std::endl; + // } + //} // Store element IDs. store_element_ids(vtk_polydata, face); - + /* + for (int e = 0; e < face.nEl; e++) { + for (int a = 0; a < face.eNoN; a++) { + int Ac = face.IEN(a,e); + Ac = face.gN(Ac); + face.IEN(a,e) = Ac; + std::cout << "[read_vtp] IEN(" << a << "," << e << "): " << face.IEN(a,e) << std::endl; + } + } + */ + std::cout << "[store vtp] Complete. " << std::endl; #ifdef debug_load_vtp std::cout << "[load_vtp] Done. " << std::endl; #endif @@ -537,7 +584,9 @@ void load_vtp(const std::string& file_name, mshType& mesh) #endif auto reader = vtkSmartPointer::New(); reader->SetFileName(file_name.c_str()); + std::cout << "[load_vtp] Reading VTK file '" << file_name << "' ... " << std::endl; reader->Update(); + std::cout << "[load_vtp] Done. " << std::endl; vtkSmartPointer vtk_polydata = reader->GetOutput(); vtkIdType num_nodes = vtk_polydata->GetNumberOfPoints(); @@ -623,6 +672,7 @@ void load_vtu(const std::string& file_name, mshType& mesh) // Store element connectivity. store_element_conn(vtk_ugrid, mesh); + } } // namespace vtk_utils From e5ff25c4568031b423708a36d4c3f2e66ef215a6 Mon Sep 17 00:00:00 2001 From: zasexton Date: Mon, 9 Oct 2023 13:52:36 -0700 Subject: [PATCH 04/15] correcting gN, gE, and IEN after reading meshes --- Code/Source/svFSI/load_msh.cpp | 50 ++++++++++++++++++++++++++++++++-- Code/Source/svFSI/vtk_xml.cpp | 6 ++-- 2 files changed, 51 insertions(+), 5 deletions(-) diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index 92402afa..beaa46c6 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -154,8 +154,7 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p // If node IDs were not read then create them. if (face.gN.size() == 0) { read_msh_ns::calc_nbc(mesh, face); - - // Reset the connecttivity with the new node IDs? + // Reset the connectivity with the new node IDs? for (int e = 0; e < face.nEl; e++) { for (int a = 0; a < face.eNoN; a++) { int Ac = face.IEN(a,e); @@ -166,6 +165,53 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p } } + // Set face global node IDs + for (int a = 0; a Date: Tue, 10 Oct 2023 10:20:41 -0700 Subject: [PATCH 05/15] removing IEN resetting in vtk_xml.cpp since this functionality is moved into load_msh.cpp --- Code/Source/svFSI/vtk_xml.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Code/Source/svFSI/vtk_xml.cpp b/Code/Source/svFSI/vtk_xml.cpp index da67b998..340d7c4f 100644 --- a/Code/Source/svFSI/vtk_xml.cpp +++ b/Code/Source/svFSI/vtk_xml.cpp @@ -420,15 +420,16 @@ void read_vtp(const std::string& file_name, faceType& face) if (face.gN.size() == 0) { std::cout << "[WARNING] No node IDs found in the '" << file_name << "' face file."; - } else { - for (int e = 0; e < face.nEl; e++) { - for (int a = 0; a < face.eNoN; a++) { + } + //else { + // for (int e = 0; e < face.nEl; e++) { + // for (int a = 0; a < face.eNoN; a++) { //int Ac = face.IEN(a,e); //Ac = face.gN(Ac); //face.IEN(a,e) = Ac; - } - } - } + // } + // } + //} // Create essential BC array. // From 8e11b754919b7655c770598f6c19aaac5c828868 Mon Sep 17 00:00:00 2001 From: zasexton Date: Thu, 12 Oct 2023 09:37:28 -0700 Subject: [PATCH 06/15] adding general parameters for unsteady darcy problem --- Code/Source/svFSI/darcy.cpp | 36 +++++++++++++++++++++++------------- 1 file changed, 23 insertions(+), 13 deletions(-) diff --git a/Code/Source/svFSI/darcy.cpp b/Code/Source/svFSI/darcy.cpp index 13fcbc4d..5891fee0 100644 --- a/Code/Source/svFSI/darcy.cpp +++ b/Code/Source/svFSI/darcy.cpp @@ -176,10 +176,18 @@ namespace darcy { const double dt = com_mod.dt; const int i = eq.s; - double k = dmn.prop.at(PhysicalProperyType::permeability); // should be a tensor + Array K = dmn.prop.at(PhysicalProperyType::permeability); // should be a tensor double source = dmn.prop.at(PhysicalProperyType::source_term); //is this just a double? - double rho = dmn.prop.at(PhysicalProperyType::solid_density); // not needed - + double rho = dmn.prop.at(PhysicalProperyType::solid_density); // need fluid density instead + //double phi = dmn.prop.at(PhysicalProperyType::porosity); + //double p0 = dmn.prop.at(PhysicalProperyType::porosity_pressure); + //double beta = dmn.prop.at(PhysicalProperyType::compressibility); + //double p1 = dmn.prop.at(PhysicalProperyType::density_pressure); + //double gamma = dmn.prop.at(PhysicalProperyType::gamma); + //double B = dmn.prop.at(PhysicalProperyType::B); + //double mu = dmn.prop.at(PhysicalProperyType::viscosity); + + // need to figure out how these are edited for the darcy case double T1 = eq.af * eq.gam * dt; double amd = eq.am * rho/ T1; double wl = w * T1; @@ -193,18 +201,18 @@ namespace darcy { dmsg << "i: " << i; dmsg << "wl: " << wl; #endif - double Td = -source; - Vector Tx(nsd); - + double Pd = -source; + Vector Px(nsd); + // need to formulate this for the compressible fluid/porous media with respect to pressure only for (int a = 0; a < eNoN; a++) { - Td = Td + N(a)*al(i,a); - Tx(0) = Tx(0) + Nx(0,a)*yl(i,a); - Tx(1) = Tx(1) + Nx(1,a)*yl(i,a); + Pd = Pd + N(a)*al(i,a); + Px(0) = Px(0) + Nx(0,a)*yl(i,a); + Px(1) = Px(1) + Nx(1,a)*yl(i,a); } - Td = Td * rho; + Pd = Pd * rho; #ifdef debug_darcy_2d dmsg; @@ -213,11 +221,13 @@ namespace darcy { #endif for (int a = 0; a < eNoN; a++) { - lR(0,a) = lR(0,a) + w*(N(a)*Td + - (Nx(0,a)*Tx(0) + Nx(1,a)*Tx(1))*k); + lR(0,a) = lR(0,a) + w*(N(a)*Pd + + (Nx(0,a)*(K(0,0)*Px(0)+K(0,1)*Px(1)) + + Nx(1,a)*(K(1,0)*Px(0)+K(1,1)*Px(1)))); for (int b = 0; b < eNoN; b++) { lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*amd + - k*(Nx(0,a)*Nx(0,b) + Nx(1,a)*Nx(1,b))); + (Nx(0,a)*(K(0,0)*Nx(0,b) + K(0,1)*Nx(1,b)) + + Nx(1,a)*(K(1,0)*Nx(0,b) + K(1,1)*Nx(1,b)))); } } } From 6621f6f4be956428db26f89d4a06966ca0f00567 Mon Sep 17 00:00:00 2001 From: zasexton Date: Sat, 14 Oct 2023 21:32:14 -0700 Subject: [PATCH 07/15] added final version of the hash map for gN, gE, and IEN setting for face meshes --- Code/Source/svFSI/load_msh.cpp | 99 +++++++++++++++++++++------------- 1 file changed, 62 insertions(+), 37 deletions(-) diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index beaa46c6..1f629023 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -11,6 +11,12 @@ #include #include #include +#include +#include +#include +#include +#include + namespace load_msh { @@ -109,7 +115,43 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p if (com_mod.ichckIEN) { read_msh_ns::check_ien(simulation, mesh); } - + std::cout << "Hash Map Construction: Nodes " << std::endl; + auto start_time = std::chrono::high_resolution_clock::now(); + // Create a hash map for nodes and elements. + std::unordered_map mesh_node_map; + for (int i = 0; i(end_time - start_time).count() << " ms" << std::endl; + std::cout << "Hash Map Construction: Nodes. Done." << std::endl; + std::cout << "Hash Map Construction: Elements " << std::endl; + start_time = std::chrono::high_resolution_clock::now(); + // Create a hash map for elements IEN + std::unordered_map mesh_element_set; + for (int i = 0; i element_nodes; + for (int k = 0; k < mesh.eNoN; k++) { + if (k != j) { + element_nodes.push_back(mesh.gIEN(k, i)); + } + } + std::sort(element_nodes.begin(), element_nodes.end()); + std::string key = ""; + for (int node: element_nodes) { + key += std::to_string(node) + ","; + } + mesh_element_set[key] = i; + } + } + end_time = std::chrono::high_resolution_clock::now(); + std::cout << "Time to set element hash map: " << std::chrono::duration_cast(end_time - start_time).count() << " ms" << std::endl; + std::cout << "Hash Map Construction: Elements. Done." << std::endl; // Read face meshes. // // Creates nodal coordinates and element connectivity data. @@ -164,54 +206,37 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p } } } - + std::cout << "face.nNo: " << face.nNo << std::endl; + start_time = std::chrono::high_resolution_clock::now(); // Set face global node IDs for (int a = 0; a(end_time - start_time).count() << " ms" << std::endl; + std::cout << "face node set. Done" << std::endl; // Set face element IDs for (int e = 0; e < face.nEl; e++) { - for (int b = 0; b < mesh.gnEl; b++) { - int foundN = 0; - for (int c = 0; c element_nodes; for (int a = 0; a < face.eNoN; a++) { int Ac = face.IEN(a,e); Ac = face.gN(Ac); face.IEN(a,e) = Ac; + element_nodes.push_back(Ac); + } + std::string key = ""; + std::sort(element_nodes.begin(), element_nodes.end()); + for (int a = 0; a(end_time - start_time).count() << " ms" << std::endl; nn::select_eleb(simulation, mesh, face); } From ca53ab35d94e86e66667864dfda4e463d6fc49ed Mon Sep 17 00:00:00 2001 From: zasexton Date: Mon, 16 Oct 2023 21:18:17 -0700 Subject: [PATCH 08/15] adding binary mask permutation for non-tet elements (#114) --- Code/Source/svFSI/load_msh.cpp | 304 +++++++++++++++++---------------- 1 file changed, 158 insertions(+), 146 deletions(-) diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index 1f629023..e414a834 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -16,7 +16,7 @@ #include #include #include - +#include namespace load_msh { @@ -89,158 +89,170 @@ void read_ndnlff(const std::string& file_name, faceType& face) /// /// SUBROUTINE READSV(list, lM) // -void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_param) -{ - auto mesh_path = mesh_param->mesh_file_path(); - auto mesh_name = mesh_param->get_name(); - #define n_dbg_read_sv - #ifdef dbg_read_sv - DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); - dmsg.banner(); - dmsg << "Mesh name: " << mesh_name; - dmsg << "Mesh path: " << mesh_path; - #endif - - // Read in volume mesh. - vtk_xml::read_vtu(mesh_path, mesh); - - // Set mesh element properites for the input element type. - nn::select_ele(simulation->com_mod, mesh); - - // Check the mesh element node ordering. - // - // Note: This may change element node ordering. - // - auto& com_mod = simulation->get_com_mod(); - if (com_mod.ichckIEN) { - read_msh_ns::check_ien(simulation, mesh); - } - std::cout << "Hash Map Construction: Nodes " << std::endl; - auto start_time = std::chrono::high_resolution_clock::now(); - // Create a hash map for nodes and elements. - std::unordered_map mesh_node_map; - for (int i = 0; i(end_time - start_time).count() << " ms" << std::endl; - std::cout << "Hash Map Construction: Nodes. Done." << std::endl; - std::cout << "Hash Map Construction: Elements " << std::endl; - start_time = std::chrono::high_resolution_clock::now(); - // Create a hash map for elements IEN - std::unordered_map mesh_element_set; - for (int i = 0; i element_nodes; - for (int k = 0; k < mesh.eNoN; k++) { - if (k != j) { - element_nodes.push_back(mesh.gIEN(k, i)); - } - } - std::sort(element_nodes.begin(), element_nodes.end()); - std::string key = ""; - for (int node: element_nodes) { - key += std::to_string(node) + ","; - } - mesh_element_set[key] = i; - } - } - end_time = std::chrono::high_resolution_clock::now(); - std::cout << "Time to set element hash map: " << std::chrono::duration_cast(end_time - start_time).count() << " ms" << std::endl; - std::cout << "Hash Map Construction: Elements. Done." << std::endl; - // Read face meshes. - // - // Creates nodal coordinates and element connectivity data. - // - mesh.nFa = mesh_param->face_parameters.size(); - mesh.fa.resize(mesh.nFa); +void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_param) { + auto mesh_path = mesh_param->mesh_file_path(); + auto mesh_name = mesh_param->get_name(); +#define n_dbg_read_sv +#ifdef dbg_read_sv + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "Mesh name: " << mesh_name; + dmsg << "Mesh path: " << mesh_path; +#endif - if (mesh.lFib && (mesh.nFa > 1)) { - throw std::runtime_error("There are " + std::to_string(mesh.nFa) + " faces defined for the '" + - mesh.name + "' mesh. Only one face is allowed for a 1D fiber-based mesh."); - } + // Read in volume mesh. + vtk_xml::read_vtu(mesh_path, mesh); - for (int i = 0; i < mesh.nFa; i++) { - auto face_param = mesh_param->face_parameters[i]; - auto& face = mesh.fa[i]; - face.name = face_param->name(); - #ifdef dbg_read_sv - dmsg << "face.name: " << face.name; - #endif - - face.qmTRI3 = face_param->quadrature_modifier_TRI3(); - if (face.qmTRI3 < (1.0/3.0) || face.qmTRI3 > 1.0) { - throw std::runtime_error("Quadrature_modifier_TRI3 must be in the range [1/3, 1]."); - } + // Set mesh element properites for the input element type. + nn::select_ele(simulation->com_mod, mesh); - if (mesh.lFib) { - auto face_path = face_param->end_nodes_face_file_path(); - #ifdef dbg_read_sv - dmsg << "Read end nodes face file ... " << " "; - dmsg << "face_path: " << face_path; - #endif - if (face_path == "") { - throw std::runtime_error("No end nodes face file path provided."); - } - - read_ndnlff(face_path, face); - - } else { - auto face_path = face_param->face_file_path(); - vtk_xml::read_vtp(face_path, face); - - // If node IDs were not read then create them. - if (face.gN.size() == 0) { - read_msh_ns::calc_nbc(mesh, face); - // Reset the connectivity with the new node IDs? - for (int e = 0; e < face.nEl; e++) { - for (int a = 0; a < face.eNoN; a++) { - int Ac = face.IEN(a,e); - Ac = face.gN(Ac); - face.IEN(a,e) = Ac; - } - } - } - } - std::cout << "face.nNo: " << face.nNo << std::endl; - start_time = std::chrono::high_resolution_clock::now(); - // Set face global node IDs - for (int a = 0; aget_com_mod(); + if (com_mod.ichckIEN) { + read_msh_ns::check_ien(simulation, mesh); } - face.gN(a) = mesh_node_map[key]; - } - end_time = std::chrono::high_resolution_clock::now(); - std::cout << "Time to set face global IDs: " << std::chrono::duration_cast(end_time - start_time).count() << " ms" << std::endl; - std::cout << "face node set. Done" << std::endl; - // Set face element IDs - for (int e = 0; e < face.nEl; e++) { - std::vector element_nodes; - for (int a = 0; a < face.eNoN; a++) { - int Ac = face.IEN(a,e); - Ac = face.gN(Ac); - face.IEN(a,e) = Ac; - element_nodes.push_back(Ac); + // Read face meshes. + // + // Creates nodal coordinates and element connectivity data. + // + mesh.nFa = mesh_param->face_parameters.size(); + mesh.fa.resize(mesh.nFa); + + if (mesh.lFib && (mesh.nFa > 1)) { + throw std::runtime_error("There are " + std::to_string(mesh.nFa) + " faces defined for the '" + + mesh.name + "' mesh. Only one face is allowed for a 1D fiber-based mesh."); } - std::string key = ""; - std::sort(element_nodes.begin(), element_nodes.end()); - for (int a = 0; aface_parameters[i]; + auto &face = mesh.fa[i]; + face.name = face_param->name(); +#ifdef dbg_read_sv + dmsg << "face.name: " << face.name; +#endif + + face.qmTRI3 = face_param->quadrature_modifier_TRI3(); + if (face.qmTRI3 < (1.0 / 3.0) || face.qmTRI3 > 1.0) { + throw std::runtime_error("Quadrature_modifier_TRI3 must be in the range [1/3, 1]."); + } + + if (mesh.lFib) { + auto face_path = face_param->end_nodes_face_file_path(); +#ifdef dbg_read_sv + dmsg << "Read end nodes face file ... " << " "; + dmsg << "face_path: " << face_path; +#endif + if (face_path == "") { + throw std::runtime_error("No end nodes face file path provided."); + } + + read_ndnlff(face_path, face); + + } else { + auto face_path = face_param->face_file_path(); + vtk_xml::read_vtp(face_path, face); + + // If node IDs were not read then create them. + if (face.gN.size() == 0) { + read_msh_ns::calc_nbc(mesh, face); + // Reset the connectivity with the new node IDs? + for (int e = 0; e < face.nEl; e++) { + for (int a = 0; a < face.eNoN; a++) { + int Ac = face.IEN(a, e); + Ac = face.gN(Ac); + face.IEN(a, e) = Ac; + } + } + } + } } - face.gE(e) = mesh_element_set[key]; - } - end_time = std::chrono::high_resolution_clock::now(); - std::cout << "Time to set face global IDs: " << std::chrono::duration_cast(end_time - start_time).count() << " ms" << std::endl; - nn::select_eleb(simulation, mesh, face); - } + if (!mesh.lFib) { + // Create a hash map for nodes and elements. + std::unordered_map mesh_node_map; + for (int i = 0; i < mesh.gnNo; i++) { + std::string key = ""; + for (int j = 0; j < com_mod.nsd; j++) { + key += std::to_string(mesh.x(j, i)) + ","; + } + mesh_node_map[key] = i; + } + // Create a hash map for elements IEN + std::unordered_map mesh_element_set; + // Generate all possible combinations of element nodes for faces + std::vector mask(mesh.eNoN); + Vector face_nodes(mesh.eNoN); + for (int i = 0; i < mesh.eNoN; i++) { + face_nodes(i) = i; + mask[i] = 0; + } + // This line calculates the number of combinations of nodes composing unique faces + int n_comb = (std::tgamma(mesh.eNoN + 1) / std::tgamma(mesh.fa[0].eNoN + 1)) * std::tgamma(mesh.eNoN - mesh.fa[0].eNoN + 1); + Array face_hash(mesh.fa[0].eNoN, n_comb); + for (int i = 0; i < mesh.fa[0].eNoN; i++) { + mask[i] = 1; + } + std::sort(mask.begin(), mask.end()); + int count = 0; + do { + int j = 0; + for (int i = 0; i < mesh.eNoN; i++) { + if (mask[i] == 1) { + face_hash(j, count) = face_nodes(i); + j++; + } + } + count++; + } while (std::next_permutation(mask.begin(), mask.end())); -} + for (int i = 0; i < mesh.gnEl; i++) { + for (int j = 0; j < n_comb; j++) { + std::vector element_nodes; + for (int k = 0; k < mesh.fa[0].eNoN; k++) { + element_nodes.push_back(mesh.gIEN(face_hash(k, j), i)); + } + std::sort(element_nodes.begin(), element_nodes.end()); + std::string key = ""; + for (int node: element_nodes) { + key += std::to_string(node) + ","; + } + mesh_element_set[key] = i; + } + } + for (int i = 0; i < mesh.nFa; i++) { + auto &face = mesh.fa[i]; + if (!mesh.lFib) { + // Set face global node IDs + for (int a = 0; a < face.nNo; a++) { + std::string key = ""; + for (int j = 0; j < com_mod.nsd; j++) { + key += std::to_string(face.x(j, a)) + ","; + } + face.gN(a) = mesh_node_map[key]; + } + // Set face element IDs + for (int e = 0; e < face.nEl; e++) { + std::vector element_nodes; + for (int a = 0; a < face.eNoN; a++) { + int Ac = face.IEN(a, e); + Ac = face.gN(Ac); + face.IEN(a, e) = Ac; + element_nodes.push_back(Ac); + } + std::string key = ""; + std::sort(element_nodes.begin(), element_nodes.end()); + for (int a = 0; a < face.eNoN; a++) { + key += std::to_string(element_nodes[a]) + ","; + } + face.gE(e) = mesh_element_set[key]; + } + nn::select_eleb(simulation, mesh, face); + } + } + } + } }; From fd04be0c476969788e54980ece00185a2c8d5ffe Mon Sep 17 00:00:00 2001 From: zasexton Date: Wed, 18 Oct 2023 18:09:44 -0700 Subject: [PATCH 09/15] adding boundary node counter-clockwise ordering to element matching (#114) --- Code/Source/svFSI/ComMod.h | 3 +++ Code/Source/svFSI/load_msh.cpp | 29 ++++++----------------- Code/Source/svFSI/vtk_xml_parser.cpp | 35 ++++++++++++++++++++++++++++ 3 files changed, 45 insertions(+), 22 deletions(-) diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h index 6482ff1b..cb23139a 100644 --- a/Code/Source/svFSI/ComMod.h +++ b/Code/Source/svFSI/ComMod.h @@ -833,6 +833,9 @@ class mshType /// @brief IB: Mesh size parameter double dx = 0.0; + /// @breif ordering: node ordering for boundaries + std::vector> ordering; + /// @brief Element distribution between processors Vector eDist; diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index e414a834..ad9bc4af 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -188,30 +188,12 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p face_nodes(i) = i; mask[i] = 0; } - // This line calculates the number of combinations of nodes composing unique faces - int n_comb = (std::tgamma(mesh.eNoN + 1) / std::tgamma(mesh.fa[0].eNoN + 1)) * std::tgamma(mesh.eNoN - mesh.fa[0].eNoN + 1); - Array face_hash(mesh.fa[0].eNoN, n_comb); - for (int i = 0; i < mesh.fa[0].eNoN; i++) { - mask[i] = 1; - } - std::sort(mask.begin(), mask.end()); - int count = 0; - do { - int j = 0; - for (int i = 0; i < mesh.eNoN; i++) { - if (mask[i] == 1) { - face_hash(j, count) = face_nodes(i); - j++; - } - } - count++; - } while (std::next_permutation(mask.begin(), mask.end())); for (int i = 0; i < mesh.gnEl; i++) { - for (int j = 0; j < n_comb; j++) { + for (unsigned int j = 0; j < mesh.ordering.size(); j++) { std::vector element_nodes; - for (int k = 0; k < mesh.fa[0].eNoN; k++) { - element_nodes.push_back(mesh.gIEN(face_hash(k, j), i)); + for (unsigned int k = 0; k < mesh.ordering[j].size(); k++) { + element_nodes.push_back(mesh.gIEN(mesh.ordering[j][k], i)); } std::sort(element_nodes.begin(), element_nodes.end()); std::string key = ""; @@ -249,10 +231,13 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p } face.gE(e) = mesh_element_set[key]; } - nn::select_eleb(simulation, mesh, face); } } } + for (int i = 0; i #include +#include namespace vtk_xml_parser { @@ -38,6 +39,32 @@ std::map vtk_cell_to_elem { {VTK_WEDGE, 6} }; +std::map>> vtk_cell_ordering { + {VTK_HEXAHEDRON, {{0,3,2,1}, + {4,5,6,7}, + {0,1,5,4}, + {1,2,6,5}, + {2,3,7,6}, + {3,0,4,7}}}, + {VTK_LINE, {{0}, + {1}}}, + {VTK_QUAD, {{0,1}, + {1,2}, + {2,3}, + {3,0}}}, + {VTK_TETRA, {{0,1,2}, + {0,1,3}, + {1,2,3}, + {2,0,3}}}, + {VTK_TRIANGLE, {{0,1}, + {1,2}, + {2,0}}}, + {VTK_WEDGE, {{0,1,2}, + {3,4,5}, + {0,1,4,3}, + {1,2,5,4}, + {2,0,3,5}}} +}; /// Names of data arrays store in VTK mesh files. const std::string NODE_IDS_NAME("GlobalNodeID"); const std::string ELEMENT_IDS_NAME("GlobalElementID"); @@ -194,18 +221,25 @@ void store_element_conn(vtkSmartPointer vtk_ugrid, mshType& #endif int np_elem = 0; + std::vector> ordering; if (num_line != 0) { np_elem = vtk_cell_to_elem[VTK_LINE]; + ordering = vtk_cell_ordering[VTK_LINE]; } if (num_hex != 0) { np_elem = vtk_cell_to_elem[VTK_HEXAHEDRON]; + ordering = vtk_cell_ordering[VTK_HEXAHEDRON]; } if (num_quad != 0) { np_elem = vtk_cell_to_elem[VTK_QUAD]; + ordering = vtk_cell_ordering[VTK_QUAD]; } if (num_tet != 0) { np_elem = vtk_cell_to_elem[VTK_TETRA]; + ordering = vtk_cell_ordering[VTK_TETRA]; } if (num_tri != 0) { np_elem = vtk_cell_to_elem[VTK_TRIANGLE]; + ordering = vtk_cell_ordering[VTK_TRIANGLE]; } if (num_wedge != 0) { np_elem = vtk_cell_to_elem[VTK_WEDGE]; + ordering = vtk_cell_ordering[VTK_WEDGE]; } // For higher-order elements with mid-side nodes. @@ -219,6 +253,7 @@ void store_element_conn(vtkSmartPointer vtk_ugrid, mshType& mesh.gnEl = num_elems; mesh.eNoN = np_elem; mesh.gIEN = Array(np_elem, num_elems); + mesh.ordering = ordering; #ifdef debug_store_element_conn std::cout << "[store_element_conn] np_elem: " << np_elem << std::endl; From 0858ce09f1337615f0c065671bcab49fcd5fec26 Mon Sep 17 00:00:00 2001 From: zasexton Date: Thu, 26 Oct 2023 18:03:40 -0700 Subject: [PATCH 10/15] completed 2D darcy general model formulation for (#71) --- Code/Source/svFSI/ComMod.h | 4 +++ Code/Source/svFSI/Parameters.cpp | 6 ++++ Code/Source/svFSI/Parameters.h | 6 ++++ Code/Source/svFSI/consts.h | 8 ++++- Code/Source/svFSI/darcy.cpp | 49 ++++++++++++++++---------- Code/Source/svFSI/read_files.cpp | 25 +++++++++++++ Code/Source/svFSI/set_equation_props.h | 7 ++++ 7 files changed, 86 insertions(+), 19 deletions(-) diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h index e83f21e1..e18835fd 100644 --- a/Code/Source/svFSI/ComMod.h +++ b/Code/Source/svFSI/ComMod.h @@ -1563,6 +1563,10 @@ class ComMod { Vector perfusion_beta0; Vector perfusion_beta1; double permeability = 0.0; + double porosity = 0.0; + double media_compressibility = 0.0; + double fluid_compressibility = 0.0; + double fluid_viscosity = 1.0; // Temporary storage for initializing state variables Vector Pinit; diff --git a/Code/Source/svFSI/Parameters.cpp b/Code/Source/svFSI/Parameters.cpp index df611081..4af3903d 100644 --- a/Code/Source/svFSI/Parameters.cpp +++ b/Code/Source/svFSI/Parameters.cpp @@ -1374,6 +1374,12 @@ DomainParameters::DomainParameters() set_parameter("Solid_viscosity", 0.9, !required, solid_viscosity); set_parameter("Source_term", 0.0, !required, source_term); set_parameter("Permeability", 0.0, !required, permeability); + set_parameter("Porosity", 0.0, !required, porosity); + set_parameter("Porosity_pressure", 0.0, !required, porosity_pressure); + set_parameter("Media_compressibility", 0.0, !required, media_compressibility); + set_parameter("Fluid_compressibility", 0.0, !required, fluid_compressibility); + set_parameter("Fluid_viscosity", 1.0, !required, fluid_viscosity); + set_parameter("Density_pressure", 0.0, !required, density_pressure); set_parameter("Time_step_for_integration", 0.0, !required, time_step_for_integration); } diff --git a/Code/Source/svFSI/Parameters.h b/Code/Source/svFSI/Parameters.h index e53ce695..d55394d6 100644 --- a/Code/Source/svFSI/Parameters.h +++ b/Code/Source/svFSI/Parameters.h @@ -1075,7 +1075,13 @@ class DomainParameters : public ParameterLists Parameter solid_viscosity; Parameter source_term; Parameter permeability; + Parameter porosity; + Parameter porosity_pressure; Parameter time_step_for_integration; + Parameter media_compressibility; + Parameter fluid_compressibility; + Parameter fluid_viscosity; + Parameter density_pressure; }; //-------------------- diff --git a/Code/Source/svFSI/consts.h b/Code/Source/svFSI/consts.h index 17cce90e..4be62fee 100644 --- a/Code/Source/svFSI/consts.h +++ b/Code/Source/svFSI/consts.h @@ -417,7 +417,13 @@ enum class PhysicalProperyType shell_thickness = 13, ctau_M = 14, // stabilization coeffs. for USTRUCT (momentum, continuity) ctau_C = 15, - permeability = 16 + permeability = 16, + porosity = 17, + porosity_pressure = 18, + media_compressibility = 19, + fluid_compressibility = 20, + fluid_viscosity = 21, + density_pressure = 22 }; //-------------------- diff --git a/Code/Source/svFSI/darcy.cpp b/Code/Source/svFSI/darcy.cpp index 5891fee0..68419cc3 100644 --- a/Code/Source/svFSI/darcy.cpp +++ b/Code/Source/svFSI/darcy.cpp @@ -175,21 +175,26 @@ namespace darcy { auto &dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; const int i = eq.s; + Array K(nsd, nsd); - Array K = dmn.prop.at(PhysicalProperyType::permeability); // should be a tensor - double source = dmn.prop.at(PhysicalProperyType::source_term); //is this just a double? + double k = dmn.prop.at(PhysicalProperyType::permeability); // should be a tensor + K(0,0) = k; + K(0,1) = 0.0; + K(1,0) = 0.0; + K(1,1) = k; + double source = dmn.prop.at(PhysicalProperyType::source_term); // is this just a double? double rho = dmn.prop.at(PhysicalProperyType::solid_density); // need fluid density instead - //double phi = dmn.prop.at(PhysicalProperyType::porosity); - //double p0 = dmn.prop.at(PhysicalProperyType::porosity_pressure); - //double beta = dmn.prop.at(PhysicalProperyType::compressibility); - //double p1 = dmn.prop.at(PhysicalProperyType::density_pressure); - //double gamma = dmn.prop.at(PhysicalProperyType::gamma); - //double B = dmn.prop.at(PhysicalProperyType::B); - //double mu = dmn.prop.at(PhysicalProperyType::viscosity); + double phi = dmn.prop.at(PhysicalProperyType::porosity); // + double p0 = dmn.prop.at(PhysicalProperyType::porosity_pressure); + double beta_0 = dmn.prop.at(PhysicalProperyType::media_compressibility); + double beta_1 = dmn.prop.at(PhysicalProperyType::fluid_compressibility); + double p1 = dmn.prop.at(PhysicalProperyType::density_pressure); + double rho_0 = dmn.prop.at(PhysicalProperyType::fluid_density); + double mu = dmn.prop.at(PhysicalProperyType::fluid_viscosity); // need to figure out how these are edited for the darcy case double T1 = eq.af * eq.gam * dt; - double amd = eq.am * rho/ T1; + double amd = eq.am / T1; double wl = w * T1; #ifdef debug_darcy_2d @@ -212,8 +217,6 @@ namespace darcy { } - Pd = Pd * rho; - #ifdef debug_darcy_2d dmsg; dmsg << "Td: " << Td; @@ -221,13 +224,23 @@ namespace darcy { #endif for (int a = 0; a < eNoN; a++) { - lR(0,a) = lR(0,a) + w*(N(a)*Pd + - (Nx(0,a)*(K(0,0)*Px(0)+K(0,1)*Px(1)) + - Nx(1,a)*(K(1,0)*Px(0)+K(1,1)*Px(1)))); + lR(0,a) = lR(0,a) + w*(N(a)*Pd*(rho_0*beta_0-2*beta_0*beta_1*yl(i,a)+beta_0*beta_1*(p0+p1) + phi*beta_1)+ + ((beta_1*yl(i,a))/mu)*(Nx(0,a)*(K(0,0)*Px(0)+K(0,1)*Px(1)) + + Nx(1,a)*(K(1,0)*Px(0)+K(1,1)*Px(1))) - + -rho_0/mu*(Nx(0,a)*(K(0,0)*Px(0)+K(0,1)*Px(1)) + + Nx(1,a)*(K(1,0)*Px(0)+K(1,1)*Px(1)))); for (int b = 0; b < eNoN; b++) { - lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*amd + - (Nx(0,a)*(K(0,0)*Nx(0,b) + K(0,1)*Nx(1,b)) + - Nx(1,a)*(K(1,0)*Nx(0,b) + K(1,1)*Nx(1,b)))); + lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*(-rho_0*beta_0-phi*beta_1+ + 2*beta_0*beta_1*yl(i,a)- + beta_0*beta_1*(p0+p1)+ + 2*beta_0*beta_1*N(b))*amd+ + Pd*N(a)*N(b)*2*beta_0*beta_1 - + ((rho_0-beta_1*(N(b)-p1))/-mu)*(Nx(0,a)*(K(0,0)*Nx(0,b) + K(0,1)*Nx(1,b)) + + Nx(1,a)*(K(1,0)*Nx(0,b) + K(1,1)*Nx(1,b))) - + ((rho_0-beta_1*(N(b)-p1))/-mu)*(Nx(0,a)*(K(0,0)*Px(0) + K(0,1)*Px(1)) + + Nx(1,a)*(K(1,0)*Px(0) + K(1,1)*Px(1))) - + ((beta_1*yl(i,a))/mu)*(Nx(0,a)*(K(0,0)*Nx(0,b) + K(0,1)*Nx(1,b)) + + Nx(1,a)*(K(1,0)*Nx(0,b) + K(1,1)*Nx(1,b)))); } } } diff --git a/Code/Source/svFSI/read_files.cpp b/Code/Source/svFSI/read_files.cpp index 44076e9d..c20fa967 100644 --- a/Code/Source/svFSI/read_files.cpp +++ b/Code/Source/svFSI/read_files.cpp @@ -1159,6 +1159,31 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& case PhysicalProperyType::permeability: rtmp = domain_params->permeability.value(); + break; + + case PhysicalProperyType::porosity: + rtmp = domain_params->porosity.value(); + break; + + case PhysicalProperyType::porosity_pressure: + rtmp = domain_params->porosity_pressure.value(); + break; + + case PhysicalProperyType::media_compressibility: + rtmp = domain_params->media_compressibility.value(); + break; + + case PhysicalProperyType::fluid_compressibility: + rtmp = domain_params->fluid_compressibility.value(); + break; + + case PhysicalProperyType::fluid_viscosity: + rtmp = domain_params->fluid_viscosity.value(); + break; + + case PhysicalProperyType::density_pressure: + rtmp = domain_params->density_pressure.value(); + break; } lEq.dmn[iDmn].prop[prop] = rtmp; diff --git a/Code/Source/svFSI/set_equation_props.h b/Code/Source/svFSI/set_equation_props.h index a7437406..c0116070 100644 --- a/Code/Source/svFSI/set_equation_props.h +++ b/Code/Source/svFSI/set_equation_props.h @@ -279,6 +279,13 @@ SetEquationPropertiesMapType set_equation_props = { propL[0][0] = PhysicalProperyType::permeability; propL[1][0] = PhysicalProperyType::source_term; propL[2][0] = PhysicalProperyType::solid_density; + propL[3][0] = PhysicalProperyType::porosity; + propL[4][0] = PhysicalProperyType::fluid_density; + propL[5][0] = PhysicalProperyType::porosity_pressure; + propL[6][0] = PhysicalProperyType::media_compressibility; + propL[7][0] = PhysicalProperyType::fluid_compressibility; + propL[8][0] = PhysicalProperyType::fluid_viscosity; + propL[9][0] = PhysicalProperyType::density_pressure; read_domain(simulation, eq_params, lEq, propL); From 7a97ef62e0b7b2e1d380230b081f62e508dce471 Mon Sep 17 00:00:00 2001 From: zasexton Date: Thu, 26 Oct 2023 18:52:16 -0700 Subject: [PATCH 11/15] adding scientific notation and 16 precision to global node matching (#114) --- Code/Source/svFSI/load_msh.cpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index 0795eb57..be95781e 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -43,8 +43,8 @@ #include #include #include -#include #include +#include #include #include @@ -203,11 +203,12 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p // Create a hash map for nodes and elements. std::unordered_map mesh_node_map; for (int i = 0; i < mesh.gnNo; i++) { - std::string key = ""; + std::ostringstream key; + key << std::scientific << std::setprecision(16); for (int j = 0; j < com_mod.nsd; j++) { - key += std::to_string(mesh.x(j, i)) + ","; + key << mesh.x(j,i) <<","; } - mesh_node_map[key] = i; + mesh_node_map[key.str()] = i; } // Create a hash map for elements IEN std::unordered_map mesh_element_set; @@ -239,11 +240,12 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p if (!mesh.lFib) { // Set face global node IDs for (int a = 0; a < face.nNo; a++) { - std::string key = ""; + std::ostringstream key; + key << std::scientific << std::setprecision(16); for (int j = 0; j < com_mod.nsd; j++) { - key += std::to_string(face.x(j, a)) + ","; + key << face.x(j, a) << ","; } - face.gN(a) = mesh_node_map[key]; + face.gN(a) = mesh_node_map[key.str()]; } // Set face element IDs for (int e = 0; e < face.nEl; e++) { From dbfd582f1a8f3ee6c5ec1d53b18d82f7eee60f90 Mon Sep 17 00:00:00 2001 From: Zachary Sexton <47196674+zasexton@users.noreply.github.com> Date: Fri, 27 Oct 2023 01:28:53 -0700 Subject: [PATCH 12/15] Revert "Merge branch 'Tissue-perfusion-model-#71' into 0-indexed-meshes(#114)" This reverts commit 7cbf43a0b21cdb5706b5292f8c5f62cd8b0d4070, reversing changes made to 7a97ef62e0b7b2e1d380230b081f62e508dce471. --- Code/Source/svFSI/CMakeLists.txt | 1 - Code/Source/svFSI/ComMod.h | 12 - Code/Source/svFSI/Parameters.cpp | 7 - Code/Source/svFSI/Parameters.h | 7 - Code/Source/svFSI/all_fun.cpp | 69 +----- Code/Source/svFSI/baf_ini.cpp | 3 +- Code/Source/svFSI/consts.cpp | 1 - Code/Source/svFSI/consts.h | 17 +- Code/Source/svFSI/darcy.cpp | 310 ------------------------- Code/Source/svFSI/darcy.h | 25 -- Code/Source/svFSI/eq_assem.cpp | 9 - Code/Source/svFSI/heats.cpp | 3 +- Code/Source/svFSI/initialize.cpp | 2 +- Code/Source/svFSI/load_msh.cpp | 2 - Code/Source/svFSI/nn.cpp | 58 +---- Code/Source/svFSI/read_files.cpp | 63 +---- Code/Source/svFSI/read_msh.cpp | 57 ++--- Code/Source/svFSI/set_equation_dof.h | 1 - Code/Source/svFSI/set_equation_props.h | 33 +-- Code/Source/svFSI/set_output_props.h | 3 +- Code/Source/svFSI/vtk_xml.cpp | 1 - Code/Source/svFSI/vtk_xml_parser.cpp | 68 +----- 22 files changed, 46 insertions(+), 706 deletions(-) delete mode 100644 Code/Source/svFSI/darcy.cpp delete mode 100644 Code/Source/svFSI/darcy.h diff --git a/Code/Source/svFSI/CMakeLists.txt b/Code/Source/svFSI/CMakeLists.txt index 56459255..4bfd2164 100644 --- a/Code/Source/svFSI/CMakeLists.txt +++ b/Code/Source/svFSI/CMakeLists.txt @@ -106,7 +106,6 @@ set(CSRCS fft.h fft.cpp heatf.h heatf.cpp heats.h heats.cpp - darcy.h darcy.cpp initialize.h initialize.cpp l_elas.h l_elas.cpp lhsa.h lhsa.cpp diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h index f3ca9408..e538751b 100644 --- a/Code/Source/svFSI/ComMod.h +++ b/Code/Source/svFSI/ComMod.h @@ -1541,18 +1541,6 @@ class ComMod { Array pSn; Vector pSa; - - /// @breif Variables for perfusion simulations - Vector perfusion_pressure_source; - Vector perfusion_pressure_sink; - Vector perfusion_beta0; - Vector perfusion_beta1; - double permeability = 0.0; - double porosity = 0.0; - double media_compressibility = 0.0; - double fluid_compressibility = 0.0; - double fluid_viscosity = 1.0; - /// @brief Temporary storage for initializing state variables Vector Pinit; Array Vinit; diff --git a/Code/Source/svFSI/Parameters.cpp b/Code/Source/svFSI/Parameters.cpp index d851c7d3..8e4a0933 100644 --- a/Code/Source/svFSI/Parameters.cpp +++ b/Code/Source/svFSI/Parameters.cpp @@ -1244,13 +1244,6 @@ DomainParameters::DomainParameters() set_parameter("Solid_density", 0.5, !required, solid_density); set_parameter("Solid_viscosity", 0.9, !required, solid_viscosity); set_parameter("Source_term", 0.0, !required, source_term); - set_parameter("Permeability", 0.0, !required, permeability); - set_parameter("Porosity", 0.0, !required, porosity); - set_parameter("Porosity_pressure", 0.0, !required, porosity_pressure); - set_parameter("Media_compressibility", 0.0, !required, media_compressibility); - set_parameter("Fluid_compressibility", 0.0, !required, fluid_compressibility); - set_parameter("Fluid_viscosity", 1.0, !required, fluid_viscosity); - set_parameter("Density_pressure", 0.0, !required, density_pressure); set_parameter("Time_step_for_integration", 0.0, !required, time_step_for_integration); } diff --git a/Code/Source/svFSI/Parameters.h b/Code/Source/svFSI/Parameters.h index 5221bb6a..a4233000 100644 --- a/Code/Source/svFSI/Parameters.h +++ b/Code/Source/svFSI/Parameters.h @@ -1018,14 +1018,7 @@ class DomainParameters : public ParameterLists Parameter solid_density; Parameter solid_viscosity; Parameter source_term; - Parameter permeability; - Parameter porosity; - Parameter porosity_pressure; Parameter time_step_for_integration; - Parameter media_compressibility; - Parameter fluid_compressibility; - Parameter fluid_viscosity; - Parameter density_pressure; }; /// @brief The RemesherParameters class stores parameters for the diff --git a/Code/Source/svFSI/all_fun.cpp b/Code/Source/svFSI/all_fun.cpp index 0b6f749f..cd62e395 100644 --- a/Code/Source/svFSI/all_fun.cpp +++ b/Code/Source/svFSI/all_fun.cpp @@ -548,7 +548,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co dmsg << "lFa.name: " << lFa.name; dmsg << "lFa.eType: " << lFa.eType; #endif - std::cout << "Integrate with Flag" << std::endl; + bool flag = pFlag; int nsd = com_mod.nsd; int insd = nsd - 1; @@ -633,29 +633,6 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co #endif double result = 0.0; - std::cout << "lFa.gN" << std::endl; - std::cout << lFa.gN << std::endl; - std::cout << "lFa.gE" << std::endl; - std::cout << lFa.gE << std::endl; - std::cout << "lFa.IEN" << std::endl; - for (int i = 0; i < lFa.nEl; i++) { - std::cout << i << " - "; - for (int j = 0; j < lFa.eNoN; j++) { - std::cout << lFa.IEN(j, i) << " "; - } - std::cout << std::endl; - } - std::cout << std::endl; - std::cout << "msh.IEN" << std::endl; - for (int i = 0; i < com_mod.msh[lFa.iM].nEl; i++) { - std::cout << i << " - "; - for (int j = 0; j < com_mod.msh[lFa.iM].eNoN; j++) { - std::cout << com_mod.msh[lFa.iM].IEN(j,i) << " "; - } - std::cout << std::endl; - } - - for (int e = 0; e < lFa.nEl; e++) { // [TODO:DaveP] not implemented. if (lFa.eType == ElementType::NRB) { @@ -745,6 +722,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co } double result = 0.0; + for (int e = 0; e < lFa.nEl; e++) { //dmsg << "----- e " << e+1 << " -----"; // Updating the shape functions, if this is a NURB @@ -920,49 +898,6 @@ double jacobian(ComMod& com_mod, const int nDim, const int eNoN, const Array //kmenon_perfusion - local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Vector& U) - { - Vector local_vector; - - if (com_mod.ltg.size() == 0) { - throw std::runtime_error("ltg is not set yet"); - } - - if (cm.mas(cm_mod)) { - if (U.size() != com_mod.gtnNo) { - throw std::runtime_error("local_rs is only specified for vector with size gtnNo"); - } - } - - if (cm.seq()) { - local_vector.resize(com_mod.gtnNo); - local_vector = U; - return local_vector; - } - - local_vector.resize(com_mod.tnNo); - Vector tmpU(com_mod.gtnNo); - - if (cm.mas(cm_mod)) { - tmpU = U; - } - - cm.bcast(cm_mod, tmpU); - - for (int a = 0; a < com_mod.tnNo; a++) { - int Ac = com_mod.ltg[a]; - local_vector[a] = tmpU[Ac]; - } - - return local_vector; - } - - Vector local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Vector& U) { diff --git a/Code/Source/svFSI/baf_ini.cpp b/Code/Source/svFSI/baf_ini.cpp index 03ea6331..a960d666 100644 --- a/Code/Source/svFSI/baf_ini.cpp +++ b/Code/Source/svFSI/baf_ini.cpp @@ -454,12 +454,11 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa) // Vector sA(com_mod.tnNo); sA = 1.0; - std::cout << "Start Area Integration" << std::endl; double area = all_fun::integ(com_mod, cm_mod, lFa, sA); - std::cout << "Face '" << lFa.name << "' area: " << area << std::endl; #ifdef debug_face_ini dmsg << "Face '" << lFa.name << "' area: " << area; #endif + if (utils::is_zero(area)) { if (cm.mas(cm_mod)) { throw std::runtime_error("Face '" + lFa.name + "' has zero area."); diff --git a/Code/Source/svFSI/consts.cpp b/Code/Source/svFSI/consts.cpp index 602315e1..795a40cb 100644 --- a/Code/Source/svFSI/consts.cpp +++ b/Code/Source/svFSI/consts.cpp @@ -188,7 +188,6 @@ const std::map equation_name_to_type = { {"heatS", EquationType::phys_heatS}, {"laplace", EquationType::phys_heatS}, {"poisson", EquationType::phys_heatS}, - {"perfusion", EquationType::phys_darcy}, {"stokes", EquationType::phys_stokes}, diff --git a/Code/Source/svFSI/consts.h b/Code/Source/svFSI/consts.h index 60bc4c91..89cabbfd 100644 --- a/Code/Source/svFSI/consts.h +++ b/Code/Source/svFSI/consts.h @@ -304,8 +304,7 @@ enum class EquationType phys_CMM = 209, phys_CEP = 210, phys_ustruct = 211, // Nonlinear elastodynamics using mixed VMS-stabilized formulation - phys_stokes = 212, - phys_darcy = 213 // Darcy flow for perfusion problems + phys_stokes = 212 }; constexpr auto Equation_CMM = EquationType::phys_CMM; @@ -320,7 +319,6 @@ constexpr auto Equation_shell = EquationType::phys_shell; constexpr auto Equation_stokes = EquationType::phys_stokes; constexpr auto Equation_struct = EquationType::phys_struct; constexpr auto Equation_ustruct = EquationType::phys_ustruct; -constexpr auto Equation_darcy = EquationType::phys_darcy; extern const std::map equation_name_to_type; @@ -361,7 +359,6 @@ enum class OutputType outGrp_fS = 523, outGrp_C = 524, outGrp_I1 = 525, - outGrp_MBF = 526, out_velocity = 599, out_pressure = 598, @@ -390,8 +387,7 @@ enum class OutputType out_viscosity = 575, out_fibStrn = 574, out_CGstrain = 573, - out_CGInv1 = 572, - out_MBF = 571 + out_CGInv1 = 572 }; /// @brief Possible physical properties. Current maxNPror is 20. @@ -413,14 +409,7 @@ enum class PhysicalProperyType damping = 12, shell_thickness = 13, ctau_M = 14, // stabilization coeffs. for USTRUCT (momentum, continuity) - ctau_C = 15, - permeability = 16, - porosity = 17, - porosity_pressure = 18, - media_compressibility = 19, - fluid_compressibility = 20, - fluid_viscosity = 21, - density_pressure = 22 + ctau_C = 15 }; enum class PreconditionerType diff --git a/Code/Source/svFSI/darcy.cpp b/Code/Source/svFSI/darcy.cpp deleted file mode 100644 index 68419cc3..00000000 --- a/Code/Source/svFSI/darcy.cpp +++ /dev/null @@ -1,310 +0,0 @@ -/* - This code implements the Darcy Equation for 2D and 3D - problems in perfusion of porus media. - ------------------------------------------------------------- - Assumptions: - - Homogeneous Permeability - - Homogeneous Density - - Isotropic Permeability - - Assumptions of Stokes Flow - - Steady-State - ------------------------------------------------------------- - Strong form of the Single-Compartment Darcy equation: - u = -K∇(P) - ∇⋅u = β0(P_source - P) - β1(P - P_sink) - where: - u -> Volume flux vector - K -> Permeability tensor - P -> Pressure - ------------------------------------------------------------- - Weak form of the Single-Compartment Darcy equation: - -∫(∇q∇P)dΩ - λ∫qPdΩ = ∫qFdΩ - ∫q∇P⋅nvdΓ - where: - q -> Test function - λ -> (β0 + β1)/K - F -> -(β0(P_source) + β1(P_sink))/K - n -> Normal vector to the boundary - -*/ - -#include "darcy.h" - -#include "all_fun.h" -#include "lhsa.h" -#include "mat_fun.h" -#include "nn.h" -#include "utils.h" - -#ifdef WITH_TRILINOS -#include "trilinos_linear_solver.h" -#endif - -namespace darcy { - //--------- - // b_darcy - //--------- - // - void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR) - { - for (int a = 0; a < eNoN; a++) { - lR(0,a) = lR(0,a) + w * N(a) * h; - } - } - - void construct_darcy(ComMod& com_mod, const mshType& lM, const Array& Ag, const Array& Yg) { -#define n_debug_construct_darcy -#ifdef debug_construct_darcy - DebugMsg dmsg(__func__, com_mod.cm.idcm()); - dmsg.banner(); -#endif - - using namespace consts; - - const int nsd = com_mod.nsd; - const int tDof = com_mod.tDof; - const int dof = com_mod.dof; - const int cEq = com_mod.cEq; - const auto &eq = com_mod.eq[cEq]; - auto &cDmn = com_mod.cDmn; - - - int eNoN = lM.eNoN; -#ifdef debug_construct_darcy - dmsg << "cEq: " << cEq; - dmsg << "cDmn: " << cDmn; -#endif - - Vector ptr(eNoN); - Vector N(eNoN); - Array xl(nsd, eNoN), al(tDof, eNoN), yl(tDof, eNoN); - Array Nx(nsd, eNoN), lR(dof, eNoN); - Array3 lK(dof * dof, eNoN, eNoN); - Array ksix(nsd, nsd); - Vector local_source(eNoN), local_sink(eNoN), local_b0(eNoN), local_b1(eNoN); - - for (int e = 0; e < lM.nEl; e++) { - cDmn = all_fun::domain(com_mod, lM, cEq, e); - auto cPhys = eq.dmn[cDmn].phys; - if (cPhys != EquationType::phys_darcy) { - std::cout << "cPhys: " << cPhys << std::endl; - continue; - } - - // Update shape function for NURBS - if (lM.eType == ElementType::NRB) { - //CALL NRMNNX(lm, e) - } - - // Create local copies - for (int a = 0; a < eNoN; a++) { - int Ac = lM.IEN(a, e); - ptr(a) = Ac; - - for (int i = 0; i < nsd; i++) { - xl(i, a) = com_mod.x(i, Ac); - } - - for (int i = 0; i < tDof; i++) { - al(i, a) = Ag(i, Ac); - yl(i, a) = Yg(i, Ac); - } - - local_source(a) = com_mod.perfusion_pressure_source(Ac); - local_sink(a) = com_mod.perfusion_pressure_sink(Ac); - local_b0(a) = com_mod.perfusion_beta0(Ac); - local_b1(a) = com_mod.perfusion_beta1(Ac); - } - - // Gauss integration - - lR = 0.0; - lK = 0.0; - double Jac{0.0}; - - for (int g = 0; g < lM.nG; g++) { - if (g == 0 || !lM.lShpF) { - auto Nx_g = lM.Nx.slice(g); - nn::gnn(eNoN, nsd, nsd, Nx_g, xl, Nx, Jac, ksix); - if (utils::is_zero(Jac)) { - throw std::runtime_error( - "[construct_darcy] Jacobian for element " + std::to_string(e) + " is < 0."); - } - } - - double w = lM.w(g) * Jac; - N = lM.N.col(g); - - if (nsd == 3) { - darcy_3d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); - } else if (nsd == 2) { - darcy_2d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK, local_source, local_sink, local_b0, local_b1); - } else { - throw std::runtime_error("[construct_darcy] nsd must be 2 or 3."); - } - } - - // Assembly -#ifdef WITH_TRILINOS - if (eq.assmTLS) { - trilinos_doassem_(const_cast(eNoN), const_cast(ptr.data()), lK.data(), lR.data()); - } -#endif - lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR); - } - } - //---------- - // darcy_2d - //---------- - // - void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, - const Array& al, const Array& yl, Array& lR, Array3& lK, - Vector& local_source, Vector& local_sink, Vector& b0, Vector& b1) { -#define n_debug_darcy_2d -#ifdef debug_darcy_2d - DebugMsg dmsg(__func__, com_mod.cm.idcm()); - dmsg.banner(); - dmsg << "w: " << w; -#endif - using namespace consts; - using namespace mat_fun; - - const int nsd = com_mod.nsd; - const int cEq = com_mod.cEq; - auto &eq = com_mod.eq[cEq]; - const int cDmn = com_mod.cDmn; - auto &dmn = eq.dmn[cDmn]; - const double dt = com_mod.dt; - const int i = eq.s; - Array K(nsd, nsd); - - double k = dmn.prop.at(PhysicalProperyType::permeability); // should be a tensor - K(0,0) = k; - K(0,1) = 0.0; - K(1,0) = 0.0; - K(1,1) = k; - double source = dmn.prop.at(PhysicalProperyType::source_term); // is this just a double? - double rho = dmn.prop.at(PhysicalProperyType::solid_density); // need fluid density instead - double phi = dmn.prop.at(PhysicalProperyType::porosity); // - double p0 = dmn.prop.at(PhysicalProperyType::porosity_pressure); - double beta_0 = dmn.prop.at(PhysicalProperyType::media_compressibility); - double beta_1 = dmn.prop.at(PhysicalProperyType::fluid_compressibility); - double p1 = dmn.prop.at(PhysicalProperyType::density_pressure); - double rho_0 = dmn.prop.at(PhysicalProperyType::fluid_density); - double mu = dmn.prop.at(PhysicalProperyType::fluid_viscosity); - - // need to figure out how these are edited for the darcy case - double T1 = eq.af * eq.gam * dt; - double amd = eq.am / T1; - double wl = w * T1; - -#ifdef debug_darcy_2d - dmsg; - dmsg << "nu: " << nu; - dmsg << "source: " << source; - dmsg << "rho: " << rho ; - dmsg << "T1: " << T1; - dmsg << "i: " << i; - dmsg << "wl: " << wl; -#endif - double Pd = -source; - Vector Px(nsd); - - // need to formulate this for the compressible fluid/porous media with respect to pressure only - for (int a = 0; a < eNoN; a++) { - Pd = Pd + N(a)*al(i,a); - Px(0) = Px(0) + Nx(0,a)*yl(i,a); - Px(1) = Px(1) + Nx(1,a)*yl(i,a); - } - - -#ifdef debug_darcy_2d - dmsg; - dmsg << "Td: " << Td; - dmsg << "Tx: " << Tx; -#endif - - for (int a = 0; a < eNoN; a++) { - lR(0,a) = lR(0,a) + w*(N(a)*Pd*(rho_0*beta_0-2*beta_0*beta_1*yl(i,a)+beta_0*beta_1*(p0+p1) + phi*beta_1)+ - ((beta_1*yl(i,a))/mu)*(Nx(0,a)*(K(0,0)*Px(0)+K(0,1)*Px(1)) + - Nx(1,a)*(K(1,0)*Px(0)+K(1,1)*Px(1))) - - -rho_0/mu*(Nx(0,a)*(K(0,0)*Px(0)+K(0,1)*Px(1)) + - Nx(1,a)*(K(1,0)*Px(0)+K(1,1)*Px(1)))); - for (int b = 0; b < eNoN; b++) { - lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*(-rho_0*beta_0-phi*beta_1+ - 2*beta_0*beta_1*yl(i,a)- - beta_0*beta_1*(p0+p1)+ - 2*beta_0*beta_1*N(b))*amd+ - Pd*N(a)*N(b)*2*beta_0*beta_1 - - ((rho_0-beta_1*(N(b)-p1))/-mu)*(Nx(0,a)*(K(0,0)*Nx(0,b) + K(0,1)*Nx(1,b)) + - Nx(1,a)*(K(1,0)*Nx(0,b) + K(1,1)*Nx(1,b))) - - ((rho_0-beta_1*(N(b)-p1))/-mu)*(Nx(0,a)*(K(0,0)*Px(0) + K(0,1)*Px(1)) + - Nx(1,a)*(K(1,0)*Px(0) + K(1,1)*Px(1))) - - ((beta_1*yl(i,a))/mu)*(Nx(0,a)*(K(0,0)*Nx(0,b) + K(0,1)*Nx(1,b)) + - Nx(1,a)*(K(1,0)*Nx(0,b) + K(1,1)*Nx(1,b)))); - } - } - } - - //---------- - // darcy_3d - //---------- - // - void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, - const Array& al, const Array& yl, Array& lR, Array3& lK) { -#define n_debug_darcy_3d -#ifdef debug_darcy_3d - DebugMsg dmsg(__func__, com_mod.cm.idcm()); - dmsg.banner(); - dmsg << "w: " << w; -#endif - using namespace consts; - using namespace mat_fun; - - const int nsd = com_mod.nsd; - const int cEq = com_mod.cEq; - auto &eq = com_mod.eq[cEq]; - const int cDmn = com_mod.cDmn; - auto &dmn = eq.dmn[cDmn]; - const double dt = com_mod.dt; - const int i = eq.s; - - double k = dmn.prop.at(PhysicalProperyType::permeability); // should be a tensor - double source = dmn.prop.at(PhysicalProperyType::source_term); //is this just a double? - double rho = dmn.prop.at(PhysicalProperyType::solid_density); - - double T1 = eq.af * eq.gam * dt; - double amd = eq.am * rho / T1; - double wl = w * T1; - -#ifdef debug_darcy_3d - dmsg; - dmsg << "k: " << k; - dmsg << "source: " << source; - dmsg << "rho: " << rho ; - dmsg << "T1: " << T1; - dmsg << "i: " << i; - dmsg << "wl: " << wl; -#endif - double Td = -source; - Vector Tx(nsd); - - for (int a = 0; a < eNoN; a++){ - Td = Td + N(a) * al(i,a); - Tx(0) = Tx(0) + Nx(0,a) * yl(i,a); - Tx(1) = Tx(1) + Nx(1,a) * yl(i,a); - Tx(2) = Tx(2) + Nx(2,a) * yl(i,a); - } - - Td = Td * rho; - - for (int a = 0; a < eNoN; a++){ - lR(0,a) = lR(0, a) + w*(N(a)*Td + - k*(Nx(0,a)*Tx(0) + Nx(1,a)*Tx(1) + Nx(2,a)*Tx(2))); - for (int b = 0; b < eNoN; b++){ - lK(0,a,b) = lK(0,a,b) + wl*(N(a)*N(b)*amd + - k*(Nx(0,a)*Nx(0,b) + Nx(1,a)*Nx(1,b) + - Nx(2,a)*Nx(2,b))); - } - } - } -} \ No newline at end of file diff --git a/Code/Source/svFSI/darcy.h b/Code/Source/svFSI/darcy.h deleted file mode 100644 index 914a3193..00000000 --- a/Code/Source/svFSI/darcy.h +++ /dev/null @@ -1,25 +0,0 @@ -// -// This code implements the Darcy Equation for 2D and 3D -// problems in perfusion of porus media. -// - -#ifndef DARCY_H -#define DARCY_H - -#include "ComMod.h" - -namespace darcy { - - void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR); - - void construct_darcy(ComMod& com_mod, const mshType& lM, const Array& Ag, const Array& Dg); - - void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, - const Array& al, const Array& yl, Array& lR, Array3& lK, Vector& local_source, - Vector& local_sink, Vector& local_b0, Vector& local_b1); - - void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, - const Array& al, const Array& yl, Array& lR, Array3& lK); -} - -#endif //DARCY_H diff --git a/Code/Source/svFSI/eq_assem.cpp b/Code/Source/svFSI/eq_assem.cpp index 02ec5688..c7d45f1e 100644 --- a/Code/Source/svFSI/eq_assem.cpp +++ b/Code/Source/svFSI/eq_assem.cpp @@ -49,7 +49,6 @@ #include "stokes.h" #include "sv_struct.h" #include "ustruct.h" -#include "darcy.h" #include @@ -134,10 +133,6 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& heats::b_heats(com_mod, eNoN, w, N, h, lR); break; - case EquationType::phys_darcy: - darcy::b_darcy(com_mod, eNoN, w, N, h, lR); - break; - case EquationType::phys_heatF: heatf::b_heatf(com_mod, eNoN, w, N, y, h, nV, lR, lK); break; @@ -363,10 +358,6 @@ void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const heats::construct_heats(com_mod, lM, Ag, Yg); break; - case EquationType::phys_darcy: - darcy::construct_darcy(com_mod, lM, Ag, Yg); - break; - case EquationType::phys_lElas: l_elas::construct_l_elas(com_mod, lM, Ag, Dg); break; diff --git a/Code/Source/svFSI/heats.cpp b/Code/Source/svFSI/heats.cpp index e4aa8a9b..a7a3b438 100644 --- a/Code/Source/svFSI/heats.cpp +++ b/Code/Source/svFSI/heats.cpp @@ -43,7 +43,6 @@ namespace heats { - void b_heats(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR) { for (int a = 0; a < eNoN; a++) { @@ -205,7 +204,7 @@ void heats_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& timeP) // std::tie(eq.dof, eq.sym) = equation_dof_map.at(eq.phys); - if (std::set{Equation_fluid, Equation_heatF, Equation_heatS, Equation_CEP, Equation_stokes, Equation_darcy}.count(eq.phys) == 0) { + if (std::set{Equation_fluid, Equation_heatF, Equation_heatS, Equation_CEP, Equation_stokes}.count(eq.phys) == 0) { dFlag = true; } diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index 35eca08b..be95781e 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -151,13 +151,11 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p mesh.nFa = mesh_param->face_parameters.size(); mesh.fa.resize(mesh.nFa); - if (mesh.lFib && (mesh.nFa > 1)) { throw std::runtime_error("There are " + std::to_string(mesh.nFa) + " faces defined for the '" + mesh.name + "' mesh. Only one face is allowed for a 1D fiber-based mesh."); } - for (int i = 0; i < mesh.nFa; i++) { auto face_param = mesh_param->face_parameters[i]; auto &face = mesh.fa[i]; diff --git a/Code/Source/svFSI/nn.cpp b/Code/Source/svFSI/nn.cpp index 455633e6..90261e9b 100644 --- a/Code/Source/svFSI/nn.cpp +++ b/Code/Source/svFSI/nn.cpp @@ -555,7 +555,6 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, auto& msh = com_mod.msh[iM]; int eNoN = msh.eNoN; - #ifdef debug_gnnb dmsg << "iM: " << iM+1; dmsg << "Ec: " << Ec+1; @@ -577,32 +576,18 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, for (int a = 0; a < eNoNb; a++) { int Ac = lFa.IEN(a,e); int b = 0; - bool match = false; + for (int ib = 0; ib < eNoN; ib++) { b = ib; if (setIt[ib]) { int Bc = msh.IEN(ib,Ec); if (Bc == Ac) { - //std::cout << "Match! " << Ac << std::endl; - match = true; break; - } else { - //std::cout << "No match! " << "Ac: " << Ac << " Bc: " << Bc << std::endl; - match = false; } } } - if (!match) { - std::cout << "Face element: " << e << std::endl; - std::cout << "Global Node : " << Ac << std::endl; - std::cout << "Mesh element: " << Ec << std::endl; - - //std::cout << "Global Nodes in Mesh Element: IEN(" << ib << "," << Ec << ")" << std::endl; - for (int ib = 0; ib < eNoN; ib++) { - std::cout << msh.IEN(ib,Ec) << " "; - } - std::cout << std::endl; + if (b > eNoN) { throw std::runtime_error("could not find matching face nodes"); } @@ -618,21 +603,14 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, a = a + 1; } } - //std::cout << "ptr("<< Ec << ")" << std::endl; - for (int i = 0; i < eNoN; i++) { - //std::cout << ptr(i) << " "; - } - //std::cout << std::endl; + // Correcting the position vector if mesh is moving // - //std::cout << "lX("<< Ec << ")" << std::endl; for (int a = 0; a < eNoN; a++) { int Ac = msh.IEN(a,Ec); for (int i = 0; i < lX.nrows(); i++) { lX(i,a) = com_mod.x(i,Ac); - //std::cout << lX(i,a) << " "; } - //std::cout << std::endl; if (com_mod.mvMsh) { for (int i = 0; i < lX.nrows(); i++) { @@ -663,7 +641,6 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, } } - auto v = utils::cross(xXi); for (int i = 0; i < nsd; i++) { v(i) = v(i) / sqrt(utils::norm(v)); @@ -707,43 +684,16 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, Array xXi(nsd,insd); - //std::cout << "Nx: " << std::endl; - for (int i = 0; isource_term.value(); break; - - case PhysicalProperyType::permeability: - rtmp = domain_params->permeability.value(); - break; - - case PhysicalProperyType::porosity: - rtmp = domain_params->porosity.value(); - break; - - case PhysicalProperyType::porosity_pressure: - rtmp = domain_params->porosity_pressure.value(); - break; - - case PhysicalProperyType::media_compressibility: - rtmp = domain_params->media_compressibility.value(); - break; - - case PhysicalProperyType::fluid_compressibility: - rtmp = domain_params->fluid_compressibility.value(); - break; - - case PhysicalProperyType::fluid_viscosity: - rtmp = domain_params->fluid_viscosity.value(); - break; - - case PhysicalProperyType::density_pressure: - rtmp = domain_params->density_pressure.value(); - break; } lEq.dmn[iDmn].prop[prop] = rtmp; @@ -1646,40 +1618,7 @@ void read_files(Simulation* simulation, const std::string& file_name) if (!com_mod.mvMsh) { throw std::runtime_error("mesh equation can only be specified after FSI equation"); } - } - - if (eq.phys == EquationType::phys_darcy) { - com_mod.perfusion_pressure_source.resize(com_mod.gtnNo); - com_mod.perfusion_pressure_sink.resize(com_mod.gtnNo); - com_mod.perfusion_beta0.resize(com_mod.gtnNo); - com_mod.perfusion_beta1.resize(com_mod.gtnNo); - std::ifstream perfusion_pressure_source_file; - std::ifstream perfusion_pressure_sink_file; - std::ifstream perfusion_pressure_b0_file; - std::ifstream perfusion_pressure_b1_file; - perfusion_pressure_source_file.open("perfusionSource"); - perfusion_pressure_sink_file.open("perfusionSink"); - perfusion_pressure_b0_file.open("perfusionB0"); - perfusion_pressure_b1_file.open("perfusionB1"); - if (!perfusion_pressure_source_file.is_open()) { - throw std::runtime_error("Failed to open the perfusion pressure source file 'perfusionSource'."); - } - if (!perfusion_pressure_sink_file.is_open()) { - throw std::runtime_error("Failed to open the perfusion pressure source file 'perfusionSource'."); - } - if (!perfusion_pressure_b0_file.is_open()) { - throw std::runtime_error("Failed to open the perfusion pressure source file 'perfusionSource'."); - } - if (!perfusion_pressure_b1_file.is_open()) { - throw std::runtime_error("Failed to open the perfusion pressure source file 'perfusionSource'."); - } - for (int i = 0; i < com_mod.gtnNo; i++) { - perfusion_pressure_source_file >> com_mod.perfusion_pressure_source[i]; - perfusion_pressure_sink_file >> com_mod.perfusion_pressure_sink[i]; - perfusion_pressure_b0_file >> com_mod.perfusion_beta0[i]; - perfusion_pressure_b1_file >> com_mod.perfusion_beta1[i]; - } - } + } } #ifdef debug_read_files dmsg << "Done Read equations " << " "; diff --git a/Code/Source/svFSI/read_msh.cpp b/Code/Source/svFSI/read_msh.cpp index ca016ce8..67e5a18f 100644 --- a/Code/Source/svFSI/read_msh.cpp +++ b/Code/Source/svFSI/read_msh.cpp @@ -679,7 +679,7 @@ void check_quad4_conn(mshType& mesh) // void check_tet_conn(mshType& mesh) { - std::cout << "[check_tet_conn] ========== check_tet_conn ========== " << std::endl; + //std::cout << "[check_tet_conn] ========== check_tet_conn ========== " << std::endl; Vector v1, v2, v3; int num_elems = mesh.gnEl; int num_reorder = 0; @@ -1117,7 +1117,6 @@ void read_msh(Simulation* simulation) // // READMSH - CALL READSV(lPM, msh(iM)) // - std::cout << "Number of Mesh Elements: " << com_mod.nMsh << std::endl; for (int iM = 0; iM < com_mod.nMsh; iM++) { auto param = simulation->parameters.mesh_parameters[iM]; auto& mesh = com_mod.msh[iM]; @@ -1154,11 +1153,9 @@ void read_msh(Simulation* simulation) #ifdef debug_read_msh dmsg << "Check for unique face names " << " ...."; #endif - std::cout << "Number of Faces: " << mesh.nFa << std::endl; for (int iFa = 0; iFa < mesh.nFa; iFa++) { mesh.fa[iFa].iM = iM; auto ctmp = mesh.fa[iFa].name; - std::cout << "Face Name: " << ctmp << std::endl; for (int i = 0; i < iM; i++) { for (int j = 0; j < com_mod.msh[i].nFa; j++) { if ((ctmp == com_mod.msh[i].fa[j].name) && ((i != iM || j != iFa))) { @@ -1167,7 +1164,7 @@ void read_msh(Simulation* simulation) } } } - std::cout << "Done Face Name Checking." << std::endl; + // Global total number of nodes. com_mod.gtnNo += mesh.gnNo; #ifdef debug_read_msh @@ -1197,7 +1194,7 @@ void read_msh(Simulation* simulation) } mesh.x.clear(); } - std::cout << "Done creating global nodal coordinate array." << std::endl; + com_mod.x = gX; // Check for shell elements. @@ -1244,7 +1241,6 @@ void read_msh(Simulation* simulation) dmsg << "Allocate new mesh nodes or something " << " ..."; dmsg << "com_mod.gtnNo: " << com_mod.gtnNo; #endif - std::cout << "gX size: " << com_mod.nsd << ", " << com_mod.gtnNo << std::endl; gX.resize(com_mod.nsd,com_mod.gtnNo); gX = com_mod.x; @@ -1261,8 +1257,7 @@ void read_msh(Simulation* simulation) msh.gN.resize(msh.gnNo); msh.gN = -1; } - } - std::cout << "Done allocating new mesh nodes or something." << std::endl; + } // Examining the existance of projection faces and setting %gN. // Reseting gtnNo and recounting nodes that are not duplicated @@ -1289,7 +1284,6 @@ void read_msh(Simulation* simulation) } } } - std::cout << "Done examining the existance of projection faces and setting %gN." << std::endl; #ifdef debug_read_msh dmsg << "Allocate com_mod.x " << " ..."; @@ -1319,7 +1313,7 @@ void read_msh(Simulation* simulation) com_mod.msh[iM].lN[Ac] = a; } } - std::cout << "Done temporarily allocating msh%lN array." << std::endl; + // Re-arranging x and finding the size of the entire domain // First rearrange 2D/3D mesh and then, 1D fiber mesh // @@ -1362,12 +1356,13 @@ void read_msh(Simulation* simulation) } b = b + 1; - } + } + #ifdef debug_read_msh dmsg << "b: " << b+1; #endif - } - std::cout << "Done rearranging x." << std::endl; + } + b = 0; for (int iM = 0; iM < com_mod.nMsh; iM++) { @@ -1391,23 +1386,24 @@ void read_msh(Simulation* simulation) } for (int i = 0; i < com_mod.nsd; i++) { - //com_mod.x(i,Ac) = gX(i,b); + com_mod.x(i,Ac) = gX(i,b); if (maxX[i] < com_mod.x(i,Ac)) { - //maxX[i] = com_mod.x(i,Ac); + maxX[i] = com_mod.x(i,Ac); } if (minX[i] > com_mod.x(i,Ac)) { - //minX[i] = com_mod.x(i,Ac); + minX[i] = com_mod.x(i,Ac); } } b = b + 1; } + #ifdef debug_read_msh dmsg << "b: " << b+1; #endif } - std::cout << "Done rearranging x." << std::endl; + #ifdef debug_read_msh dmsg << "minX: " << minX; dmsg << "maxX: " << maxX; @@ -1419,32 +1415,24 @@ void read_msh(Simulation* simulation) dmsg << "Rearrange " << " fa ... "; #endif - - std::cout << "Number of Meshes: " << com_mod.nMsh << std::endl; for (int iM = 0; iM < com_mod.nMsh; iM++) { for (int iFa = 0; iFa < com_mod.msh[iM].nFa; iFa++) { for (int a = 0; a < com_mod.msh[iM].fa[iFa].nNo; a++) { - //int Ac = com_mod.msh[iM].fa[iFa].gN[a]; - //std::cout << "Ac: " << Ac << std::endl; - //Ac = com_mod.msh[iM].gN[Ac]; - //com_mod.msh[iM].fa[iFa].gN[a] = Ac; + int Ac = com_mod.msh[iM].fa[iFa].gN[a]; + Ac = com_mod.msh[iM].gN[Ac]; + com_mod.msh[iM].fa[iFa].gN[a] = Ac; } - std::cout << "Done rearranging fa inner most loop 1." << std::endl; + for (int e = 0; e < com_mod.msh[iM].fa[iFa].nEl; e++) { for (int a = 0; a < com_mod.msh[iM].fa[iFa].eNoN; a++) { - //int Ac = com_mod.msh[iM].fa[iFa].IEN(a,e); - //Ac = com_mod.msh[iM].gN[Ac]; - //std::cout << "IEN(" << a << "," << e << "): " << Ac << std::endl; - //com_mod.msh[iM].fa[iFa].IEN(a,e) = Ac; - //std::cout << "IEN(" << a << "," << e << "): " << Ac << " Done." << std::endl; + int Ac = com_mod.msh[iM].fa[iFa].IEN(a,e); + Ac = com_mod.msh[iM].gN[Ac]; + com_mod.msh[iM].fa[iFa].IEN(a,e) = Ac; } - } - std::cout << "Done rearranging fa inner most loop 2." << std::endl; + } } } - std::cout << "Done rearranging fa." << std::endl; - if (com_mod.resetSim) { auto& rmsh = com_mod.rmsh; int gtnNo = com_mod.gtnNo; @@ -1482,7 +1470,6 @@ void read_msh(Simulation* simulation) } } - std::cout << "Done resetSim." << std::endl; // Setting dmnId parameter here, if there is at least one mesh that // has defined eId. // diff --git a/Code/Source/svFSI/set_equation_dof.h b/Code/Source/svFSI/set_equation_dof.h index 7033adec..a7cddb85 100644 --- a/Code/Source/svFSI/set_equation_dof.h +++ b/Code/Source/svFSI/set_equation_dof.h @@ -41,7 +41,6 @@ std::map equation_dof_map = {EquationType::phys_fluid, std::make_tuple(nsd+1, "NS") }, {EquationType::phys_heatF, std::make_tuple(1, "HF") }, {EquationType::phys_heatS, std::make_tuple(1, "HS") }, - {EquationType::phys_darcy, std::make_tuple(1, "HS") }, {EquationType::phys_lElas, std::make_tuple(nsd, "LE") }, {EquationType::phys_struct, std::make_tuple(nsd, "ST") }, {EquationType::phys_ustruct, std::make_tuple(nsd+1, "ST") }, diff --git a/Code/Source/svFSI/set_equation_props.h b/Code/Source/svFSI/set_equation_props.h index db89c75e..74aeea31 100644 --- a/Code/Source/svFSI/set_equation_props.h +++ b/Code/Source/svFSI/set_equation_props.h @@ -298,38 +298,7 @@ SetEquationPropertiesMapType set_equation_props = { } }, //---------------------------// -// phys_darcy // -//---------------------------// -{consts::EquationType::phys_darcy, [](Simulation* simulation, EquationParameters* eq_params, eqType& lEq, EquationProps& propL, - EquationOutputs& outPuts, EquationNdop& nDOP) -> void -{ - using namespace consts; - auto& com_mod = simulation->get_com_mod(); - lEq.phys = consts::EquationType::phys_darcy; - - propL[0][0] = PhysicalProperyType::permeability; - propL[1][0] = PhysicalProperyType::source_term; - propL[2][0] = PhysicalProperyType::solid_density; - propL[3][0] = PhysicalProperyType::porosity; - propL[4][0] = PhysicalProperyType::fluid_density; - propL[5][0] = PhysicalProperyType::porosity_pressure; - propL[6][0] = PhysicalProperyType::media_compressibility; - propL[7][0] = PhysicalProperyType::fluid_compressibility; - propL[8][0] = PhysicalProperyType::fluid_viscosity; - propL[9][0] = PhysicalProperyType::density_pressure; - - read_domain(simulation, eq_params, lEq, propL); - - nDOP = {2,1,1,0}; - outPuts = {OutputType::out_temperature, OutputType::out_heatFlux}; - - // Set solver parameters. - read_ls(simulation, eq_params, SolverType::lSolver_CG, lEq); -} }, - - -//---------------------------// -// phys_FSI //const +// phys_FSI // //---------------------------// // {consts::EquationType::phys_FSI, [](Simulation* simulation, EquationParameters* eq_params, eqType& lEq, EquationProps& propL, diff --git a/Code/Source/svFSI/set_output_props.h b/Code/Source/svFSI/set_output_props.h index 0076e177..5c128880 100644 --- a/Code/Source/svFSI/set_output_props.h +++ b/Code/Source/svFSI/set_output_props.h @@ -82,7 +82,6 @@ std::map output_props_map = {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")}, - {OutputType::out_MBF, std::make_tuple(OutputType::outGrp_MBF, 0, 1, "MBF")} + {OutputType::out_WSS, std::make_tuple(OutputType::outGrp_WSS, 0, maxNSD, "WSS") } }; diff --git a/Code/Source/svFSI/vtk_xml.cpp b/Code/Source/svFSI/vtk_xml.cpp index c38e3d37..c8626338 100644 --- a/Code/Source/svFSI/vtk_xml.cpp +++ b/Code/Source/svFSI/vtk_xml.cpp @@ -479,7 +479,6 @@ void read_vtp(const std::string& file_name, faceType& face) //std::cout << std::endl; } } - } //---------------- diff --git a/Code/Source/svFSI/vtk_xml_parser.cpp b/Code/Source/svFSI/vtk_xml_parser.cpp index cde8e70c..fff079e9 100644 --- a/Code/Source/svFSI/vtk_xml_parser.cpp +++ b/Code/Source/svFSI/vtk_xml_parser.cpp @@ -137,22 +137,14 @@ void store_element_conn(vtkSmartPointer vtk_polydata, faceType& fac face.eNoN = np_elem; face.IEN = Array(np_elem, num_elems); - auto global_node_ar = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); for (int i = 0; i < num_elems; i++) { vtk_polydata->GetCell(i, cell); - // Map the Extracted Local Cell Point Ids to the Global Volume Mesh Node Ids - //std::cout << "Begin Safe Down Cast from the Data Storage" << std::endl; - //auto global_node_ar = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); auto num_cell_pts = cell->GetNumberOfPoints(); - //std::cout << "Success on DownCast Operation" << std::endl; for (int j = 0; j < num_cell_pts; j++) { - auto id = global_node_ar->GetValue(cell->PointIds->GetId(j)); - std::cout << "Id type: " << static_cast(id) <<" id: " << id << std::endl; + auto id = cell->PointIds->GetId(j); face.IEN(j,i) = id; - //std::cout << "face_" < vtk_polydata, mshType& mesh) @@ -177,15 +169,12 @@ void store_element_conn(vtkSmartPointer vtk_polydata, mshType& mesh std::cout << "[store_element_conn(polydata,mesh)] np_elem: " << np_elem << std::endl; #endif - auto global_node_ar = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); for (int i = 0; i < num_elems; i++) { - // int elem_id = elem_ids->GetValue(i); - // Node id's within the total mesh are the same since interior elements are not neglected + //int elem_id = elem_ids->GetValue(i); vtk_polydata->GetCell(i, cell); auto num_cell_pts = cell->GetNumberOfPoints(); for (int j = 0; j < num_cell_pts; j++) { - //auto id = cell->PointIds->GetId(j); - auto id = global_node_ar->GetValue(cell->PointIds->GetId(j)); + auto id = cell->PointIds->GetId(j); mesh.gIEN(j,i) = id; } } @@ -206,7 +195,6 @@ void store_element_conn(vtkSmartPointer vtk_ugrid, mshType& std::cout << "[store_element_conn(ugrid)] " << std::endl; std::cout << "[store_element_conn(ugrid)] ========== store_element_conn(ugrid) =========" << std::endl; #endif - std::cout << "[store_element_conn(polydata,mesh)] connecting mesh" << std::endl; vtkSmartPointer cell_types = vtk_ugrid->GetCellTypesArray(); auto num_elems = vtk_ugrid->GetNumberOfCells(); #ifdef debug_store_element_conn @@ -302,7 +290,7 @@ void store_element_conn(vtkSmartPointer vtk_ugrid, mshType& std::cout << "[store_element_conn] Number of elements: " << num_elems << std::endl; std::cout << "[store_element_conn] Number of nodes per element: " << np_elem << std::endl; #endif - auto global_node_ar = vtkIntArray::SafeDownCast(vtk_ugrid->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); + auto cell = vtkGenericCell::New(); for (int i = 0; i < num_elems; i++) { vtk_ugrid->GetCell(i, cell); @@ -313,9 +301,6 @@ void store_element_conn(vtkSmartPointer vtk_ugrid, mshType& } for (int j = 0; j < num_cell_pts; j++) { auto id = cell->PointIds->GetId(j); - //auto id2 = global_node_ar->GetValue(id); - //std::cout << "mesh.gIEN("<< j << "," << i << ") = " << id << std::endl; - // std::cout << "mesh.gIEN("<< j << "," << i << ") = " << id2 << std::endl; //mesh.gIEN(i,j) = id; mesh.gIEN(j,i) = id; } @@ -356,8 +341,7 @@ void store_element_ids(vtkSmartPointer vtk_polydata, faceType& face // [NOTE] It is not clear how these IDs are used but if they // index into arrays or vectors then they need to be offset by -1. for (int i = 0; i < num_elem_ids; i++) { - face.gE(i) = elem_ids->GetValue(i); // - 1; assume that the face meshes are zero indexed - //std::cout << "[store_element_ids] face.gE(" << i << "): " << face.gE(i) << std::endl; + face.gE(i) = elem_ids->GetValue(i) - 1; } } @@ -433,24 +417,14 @@ void store_nodal_ids(vtkSmartPointer vtk_ugrid, mshType& me void store_nodal_ids(vtkSmartPointer vtk_polydata, faceType& face) { vtkIdType num_nodes = vtk_polydata->GetNumberOfPoints(); - vtkIdType num_pt_arrays = vtk_polydata->GetPointData()->GetNumberOfArrays(); - vtkIdType num_cell_arrays = vtk_polydata->GetCellData()->GetNumberOfArrays(); auto node_ids = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(NODE_IDS_NAME.c_str())); - //for (vtkIdType i = 0; i < num_pt_arrays; i++) { - // if (vtk_polydata->GetPointData()->GetArrayName(i) == NODE_IDS_NAME){ - // //std::cout << "GlobalNodeID Array Number: " << i << std::endl; - // node_ids = vtkIntArray::SafeDownCast(vtk_polydata->GetPointData()->GetArray(i)); - // //std::cout << node_ids << std::endl; - // } - //} if (node_ids == nullptr) { return; } face.gN = Vector(num_nodes); for (int i = 0; i < num_nodes; i++) { - //std::cout << "[store_nodal_ids] node_ids->GetValue(" << i << "): " << std::endl; // [NOTE] It seems that face node IDs are 1 based. - face.gN(i) = node_ids->GetValue(i); //- 1; + face.gN(i) = node_ids->GetValue(i) - 1; //std::cout << "[store_nodal_ids] face.gN(" << i << "): " << face.gN(i) << std::endl; } } @@ -465,7 +439,7 @@ void store_nodal_ids(vtkSmartPointer vtk_polydata, mshType& mesh) mesh.gN = Vector(num_nodes); for (int i = 0; i < num_nodes; i++) { // [NOTE] It seems that face node IDs are 1 based. - mesh.gN(i) = node_ids->GetValue(i); // - 1; + mesh.gN(i) = node_ids->GetValue(i) - 1; //std::cout << "[store_nodal_ids] face.gN(" << i << "): " << face.gN(i) << std::endl; } } @@ -552,10 +526,8 @@ void load_vtp(const std::string& file_name, faceType& face) std::cout << "[load_vtp] ===== vtk_xml_parser.cpp::load_vtp ===== " << std::endl; #endif auto reader = vtkSmartPointer::New(); - std::cout << "[load_vtp] Reading VTK file '" << file_name << "' ... " << std::endl; reader->SetFileName(file_name.c_str()); reader->Update(); - std::cout << "[load_vtp] Done. " << std::endl; vtkSmartPointer vtk_polydata = reader->GetOutput(); vtkIdType num_nodes = vtk_polydata->GetNumberOfPoints(); @@ -564,8 +536,6 @@ void load_vtp(const std::string& file_name, faceType& face) } vtkIdType num_elems = vtk_polydata->GetNumberOfCells(); - std::cout << "[load_vtp] Number of nodes: " << num_nodes << std::endl; - std::cout << "[load_vtp] Number of elements: " << num_elems << std::endl; #ifdef debug_load_vtp std::cout << "[load_vtp] Number of nodes: " << num_nodes << std::endl; std::cout << "[load_vtp] Number of elements: " << num_elems << std::endl; @@ -580,27 +550,10 @@ void load_vtp(const std::string& file_name, faceType& face) // Store element connectivity. store_element_conn(vtk_polydata, face); - //for (int e = 0; e < face.nEl; e++) { - // for (int a = 0; a < face.eNoN; a++) { - //int Ac = face.IEN(a,e); - //Ac = face.gN(Ac); - //face.IEN(a,e) = Ac; - // std::cout << "[read_vtp] IEN(" << a << "," << e << "): " << face.IEN(a,e) << std::endl; - // } - //} + // Store element IDs. store_element_ids(vtk_polydata, face); - /* - for (int e = 0; e < face.nEl; e++) { - for (int a = 0; a < face.eNoN; a++) { - int Ac = face.IEN(a,e); - Ac = face.gN(Ac); - face.IEN(a,e) = Ac; - std::cout << "[read_vtp] IEN(" << a << "," << e << "): " << face.IEN(a,e) << std::endl; - } - } - */ - std::cout << "[store vtp] Complete. " << std::endl; + #ifdef debug_load_vtp std::cout << "[load_vtp] Done. " << std::endl; #endif @@ -616,9 +569,7 @@ void load_vtp(const std::string& file_name, mshType& mesh) #endif auto reader = vtkSmartPointer::New(); reader->SetFileName(file_name.c_str()); - std::cout << "[load_vtp] Reading VTK file '" << file_name << "' ... " << std::endl; reader->Update(); - std::cout << "[load_vtp] Done. " << std::endl; vtkSmartPointer vtk_polydata = reader->GetOutput(); vtkIdType num_nodes = vtk_polydata->GetNumberOfPoints(); @@ -701,7 +652,6 @@ void load_vtu(const std::string& file_name, mshType& mesh) // Store element connectivity. store_element_conn(vtk_ugrid, mesh); - } } // namespace vtk_utils From ca9ff014de8c055e8f941498b560c593e8c0452c Mon Sep 17 00:00:00 2001 From: Zachary Sexton <47196674+zasexton@users.noreply.github.com> Date: Fri, 3 Nov 2023 12:21:20 -0700 Subject: [PATCH 13/15] encapsulating hash map generation (#114) --- Code/Source/svFSI/load_msh.cpp | 49 ++++++++++++++++++++++++++++++++-- 1 file changed, 47 insertions(+), 2 deletions(-) diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index be95781e..6d3a5ea3 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -50,6 +50,46 @@ namespace load_msh { + /// @brief Generate hash maps for mesh nodes and elements + +class MeshHashMaps { +public: + MeshHashMaps() {} + + std::unordered_map createNodeHashMap(const mshType& mesh, const int& nsd) { + std::unordered_map mesh_node_map; + for (int i = 0; i < mesh.gnNo; i++) { + std::ostringstream key; + key << std::scientific << std::setprecision(16); + for (int j = 0; j < nsd; j++) { + key << mesh.x(j,i) <<","; + } + mesh_node_map[key.str()] = i; + } + return mesh_node_map; + } + + std::unordered_map createElementHashMap(const mshType& mesh) { + std::unordered_map mesh_element_set; + for (int i = 0; i < mesh.gnEl; i++) { + for (unsigned int j = 0; j < mesh.ordering.size(); j++) { + std::vector element_nodes; + for (unsigned int k = 0; k < mesh.ordering[j].size(); k++) { + element_nodes.push_back(mesh.gIEN(mesh.ordering[j][k], i)); + } + std::sort(element_nodes.begin(), element_nodes.end()); + std::string key = ""; + for (int node: element_nodes) { + key += std::to_string(node) + ","; + } + mesh_element_set[key] = i; + } + } + return mesh_element_set; + } +}; + + #define ndbg_load_msh /// @brief Read mesh position coordinates and connectivity. @@ -60,7 +100,7 @@ void read_ccne(Simulation* simulation, mshType& mesh, const MeshParameters* mesh { auto mesh_path = mesh_param->mesh_file_path(); auto mesh_name = mesh_param->name(); - throw std::runtime_error("[read_ccne] read_ccne() is not implemented."); + throw std::runtime_error("[read_ccne] read_ccne() is not implemented."); } /// @brief Read list of end nodes from a file into the face data structure. @@ -201,6 +241,10 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p } if (!mesh.lFib) { // Create a hash map for nodes and elements. + MeshHashMaps mesh_hash_maps; + auto mesh_node_map = mesh_hash_maps.createNodeHashMap(mesh, com_mod.nsd); + auto mesh_element_set = mesh_hash_maps.createElementHashMap(mesh); + /* std::unordered_map mesh_node_map; for (int i = 0; i < mesh.gnNo; i++) { std::ostringstream key; @@ -233,7 +277,8 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p } mesh_element_set[key] = i; } - } + */ + //} for (int i = 0; i < mesh.nFa; i++) { auto &face = mesh.fa[i]; From 7ed7fd026998f0cbfaf1a557eb6139196aa29243 Mon Sep 17 00:00:00 2001 From: Zachary Sexton <47196674+zasexton@users.noreply.github.com> Date: Fri, 3 Nov 2023 14:09:11 -0700 Subject: [PATCH 14/15] remove commented lines performing un-encapsulated hash map generation --- Code/Source/svFSI/load_msh.cpp | 48 +++++++++------------------------- 1 file changed, 13 insertions(+), 35 deletions(-) diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index 6d3a5ea3..2d06df42 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -56,6 +56,13 @@ class MeshHashMaps { public: MeshHashMaps() {} + /// @brief Create a hash map for the global mesh nodes + /// + /// \code {.cpp} + /// mesh + /// nsd + /// \endcode + /// std::unordered_map createNodeHashMap(const mshType& mesh, const int& nsd) { std::unordered_map mesh_node_map; for (int i = 0; i < mesh.gnNo; i++) { @@ -69,6 +76,12 @@ class MeshHashMaps { return mesh_node_map; } + /// @brief Create a hash map for the global mesh elements + /// + /// \code {.cpp} + /// mesh + /// \endcode + /// std::unordered_map createElementHashMap(const mshType& mesh) { std::unordered_map mesh_element_set; for (int i = 0; i < mesh.gnEl; i++) { @@ -244,41 +257,6 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p MeshHashMaps mesh_hash_maps; auto mesh_node_map = mesh_hash_maps.createNodeHashMap(mesh, com_mod.nsd); auto mesh_element_set = mesh_hash_maps.createElementHashMap(mesh); - /* - std::unordered_map mesh_node_map; - for (int i = 0; i < mesh.gnNo; i++) { - std::ostringstream key; - key << std::scientific << std::setprecision(16); - for (int j = 0; j < com_mod.nsd; j++) { - key << mesh.x(j,i) <<","; - } - mesh_node_map[key.str()] = i; - } - // Create a hash map for elements IEN - std::unordered_map mesh_element_set; - // Generate all possible combinations of element nodes for faces - std::vector mask(mesh.eNoN); - Vector face_nodes(mesh.eNoN); - for (int i = 0; i < mesh.eNoN; i++) { - face_nodes(i) = i; - mask[i] = 0; - } - - for (int i = 0; i < mesh.gnEl; i++) { - for (unsigned int j = 0; j < mesh.ordering.size(); j++) { - std::vector element_nodes; - for (unsigned int k = 0; k < mesh.ordering[j].size(); k++) { - element_nodes.push_back(mesh.gIEN(mesh.ordering[j][k], i)); - } - std::sort(element_nodes.begin(), element_nodes.end()); - std::string key = ""; - for (int node: element_nodes) { - key += std::to_string(node) + ","; - } - mesh_element_set[key] = i; - } - */ - //} for (int i = 0; i < mesh.nFa; i++) { auto &face = mesh.fa[i]; From 3555c550690132b4c442abe1cf0962330fadad89 Mon Sep 17 00:00:00 2001 From: Zachary Sexton <47196674+zasexton@users.noreply.github.com> Date: Mon, 6 Nov 2023 09:29:29 -0800 Subject: [PATCH 15/15] adding error checking for duplicate nodes during hash map generation (#114) --- Code/Source/svFSI/load_msh.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index 2d06df42..f46b60c1 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -73,6 +73,9 @@ class MeshHashMaps { } mesh_node_map[key.str()] = i; } + if (mesh_node_map.size() != mesh.gnNo) { + throw std::runtime_error("There may be a duplicate nodes within the mesh."); + } return mesh_node_map; }