From e1bec2f63c8f107d95fcb28405aa62f7329b0edc Mon Sep 17 00:00:00 2001 From: KB Date: Thu, 14 May 2026 13:53:11 -0700 Subject: [PATCH 01/15] Added comments for fib_strech for readability --- Code/Source/solver/post.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 57bbe6194..80f967187 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -794,8 +794,9 @@ void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const double w = lM.w(g)*Jac; auto N = lM.N.col(g); - F = mat_fun::mat_id(nsd); + // Compute Deformation Gradient: F = I + grad(u) + F = mat_fun::mat_id(nsd); for (int a = 0; a < eNoN; a++) { if (nsd == 3) { F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); @@ -815,9 +816,11 @@ void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const } } + // Compute fiber stretch based on 4th invariant: I_{4,f} = F.fN.F.fN auto fl = mat_fun::mat_mul(F, lM.fN.rows(0,nsd-1,e)); double I4f = utils::norm(fl); + // L2 projection of I4f from integration points to nodes for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); sA(Ac) = sA(Ac) + w*N(a); From 1b7e8271b2d8822e199578d449d7ee87ff746569 Mon Sep 17 00:00:00 2001 From: KB Date: Thu, 14 May 2026 15:58:10 -0700 Subject: [PATCH 02/15] added stretch rate and minor fixes --- Code/Source/solver/cep_ion.cpp | 2 +- Code/Source/solver/post.cpp | 158 ++++++++++++++++++++++++++++++++- Code/Source/solver/post.h | 4 +- 3 files changed, 159 insertions(+), 5 deletions(-) diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 45c427817..fb7075b43 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -157,7 +157,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt if (msh.nFn != 0) { Vector sA(msh.nNo); - post::fib_strech(simulation, iEq, msh, solutions, sA); + post::fib_stretch(simulation, iEq, msh, solutions, sA); for (int a = 0; a < msh.nNo; a++) { int Ac = msh.gN(a); I4f(Ac) = sA(a); diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 80f967187..5152e66be 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -742,7 +742,7 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra /// @brief Compute fiber stretch based on 4th invariant: I_{4,f} // -void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) +void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) { using namespace consts; @@ -787,9 +787,10 @@ void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const for (int g = 0; g < lM.nG; g++) { double Jac = 0.0; Array F(nsd,nsd); + Array ksix(nsd,nsd); if (g == 0 || !lM.lShpF) { - auto Nx = lM.Nx.slice(g); - nn::gnn(eNoN, nsd, nsd, Nx, xl, Nx, Jac, F); + auto Nxi = lM.Nx.slice(g); + nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, ksix); } double w = lM.w(g)*Jac; @@ -843,6 +844,157 @@ void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const } +/// @brief Compute fiber stretch rate dλ/dt = (1/2λ) fN^T (dC/dt) fN +// +void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) +{ + using namespace consts; + + auto& com_mod = simulation->com_mod; + auto& cm = com_mod.cm; + auto& cm_mod = simulation->cm_mod; + const auto& lD = solutions.old.get_displacement(); + const auto& lY = solutions.current.get_velocity(); + auto& eq = com_mod.eq[iEq]; + + int nsd = com_mod.nsd; + int tnNo = com_mod.tnNo; + int tDof = com_mod.tDof; + + // [NOTE] Setting gobal variable 'dof'. + com_mod.dof = eq.dof; + + int eNoN = lM.eNoN; + int i = eq.s; + int j = i + 1; + int k = j + 1; + + Vector sA(tnNo); + Vector sF(tnNo); + Array xl(nsd,eNoN); + Array dl(tDof,eNoN); + Array yl(tDof,eNoN); + Array Nx(nsd,eNoN); + + for (int e = 0; e < lM.nEl; e++) { + int cDmn = all_fun::domain(com_mod, lM, iEq, e); + auto cPhys = eq.dmn[cDmn].phys; + if (cPhys != EquationType::phys_struct && cPhys != EquationType::phys_ustruct) { + continue; + } + if (lM.eType == ElementType::NRB) { + //CALL NRBNNX(lM, e) + } + + for (int a = 0; a < eNoN; a++) { + int Ac = lM.IEN(a,e); + xl.set_col(a, com_mod.x.col(Ac)); + dl.set_col(a, lD.col(Ac)); + for (int i = 0; i < tDof; i++) { + yl(i,a) = lY(i,Ac); + } + } + + for (int g = 0; g < lM.nG; g++) { + double Jac = 0.0; + Array F(nsd,nsd); + Array vx(nsd,nsd); + Array ksix(nsd,nsd); + + if (g == 0 || !lM.lShpF) { + auto Nxi = lM.Nx.slice(g); + nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, ksix); + } + + double w = lM.w(g)*Jac; + auto N = lM.N.col(g); + F = mat_fun::mat_id(nsd); + + for (int a = 0; a < eNoN; a++) { + if (nsd == 3) { + F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); + F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); + F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); + F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); + F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); + F(1,2) = F(1,2) + Nx(2,a)*dl(j,a); + F(2,0) = F(2,0) + Nx(0,a)*dl(k,a); + F(2,1) = F(2,1) + Nx(1,a)*dl(k,a); + F(2,2) = F(2,2) + Nx(2,a)*dl(k,a); + + vx(0,0) = vx(0,0) + Nx(0,a)*yl(i,a); + vx(0,1) = vx(0,1) + Nx(1,a)*yl(i,a); + vx(0,2) = vx(0,2) + Nx(2,a)*yl(i,a); + vx(1,0) = vx(1,0) + Nx(0,a)*yl(j,a); + vx(1,1) = vx(1,1) + Nx(1,a)*yl(j,a); + vx(1,2) = vx(1,2) + Nx(2,a)*yl(j,a); + vx(2,0) = vx(2,0) + Nx(0,a)*yl(k,a); + vx(2,1) = vx(2,1) + Nx(1,a)*yl(k,a); + vx(2,2) = vx(2,2) + Nx(2,a)*yl(k,a); + } else { + F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); + F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); + F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); + F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); + + vx(0,0) = vx(0,0) + Nx(0,a)*yl(i,a); + vx(0,1) = vx(0,1) + Nx(1,a)*yl(i,a); + vx(1,0) = vx(1,0) + Nx(0,a)*yl(j,a); + vx(1,1) = vx(1,1) + Nx(1,a)*yl(j,a); + } + } + + auto fl = mat_fun::mat_mul(F, lM.fN.rows(0,nsd-1,e)); + double I4f = utils::norm(fl); + double lambda = sqrt(I4f); + if (lambda < 1e-12) { + lambda = 1.0; + } + + // dC/dt = (vx)^T F + F^T vx + Array Vx_T = mat_fun::transpose(vx); + Array F_T = mat_fun::transpose(F); + Array dC_dt(nsd, nsd); + auto Vx_T_F = mat_fun::mat_mul(Vx_T, F); + auto F_T_Vx = mat_fun::mat_mul(F_T, vx); + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + dC_dt(i,j) = Vx_T_F(i,j) + F_T_Vx(i,j); + } + } + + auto fN = lM.fN.rows(0,nsd-1,e); + double fT_dC_dt_f = 0.0; + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + fT_dC_dt_f += fN(i) * dC_dt(i,j) * fN(j); + } + } + + double dlambda_dt = fT_dC_dt_f / (2.0 * lambda); + + for (int a = 0; a < eNoN; a++) { + int Ac = lM.IEN(a,e); + sA(Ac) = sA(Ac) + w*N(a); + sF(Ac) = sF(Ac) + w*N(a)*dlambda_dt; + } + } + } + + all_fun::commu(com_mod, sF); + all_fun::commu(com_mod, sA); + + res = 0.0; + + for (int a = 0; a < lM.nNo; a++) { + int Ac = lM.gN(a); + if (!utils::is_zero(sA(Ac))) { + res(a) = res(a) + sF(Ac) / sA(Ac); + } + } + +} + void post(Simulation* simulation, const mshType& lM, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp, const int iEq) { diff --git a/Code/Source/solver/post.h b/Code/Source/solver/post.h index b40daae87..ee9616f12 100644 --- a/Code/Source/solver/post.h +++ b/Code/Source/solver/post.h @@ -22,7 +22,9 @@ void fib_algn_post(Simulation* simulation, const mshType& lM, Array& res void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Array& res, const SolutionStates& solutions, const int iEq); -void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); +void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); + +void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); void post(Simulation* simulation, const mshType& lM, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp, const int iEq); From ae9d115b90e13d1f29896ebdc981bc70213fa86d Mon Sep 17 00:00:00 2001 From: KB Date: Fri, 15 May 2026 15:57:42 -0700 Subject: [PATCH 03/15] fiber stretch test case added --- Code/Source/solver/consts.h | 8 +- Code/Source/solver/set_equation_props.h | 13 +- Code/Source/solver/set_output_props.h | 2 + Code/Source/solver/vtk_xml.cpp | 22 ++- .../blocks-mesh-complete/block_domains.vtu | 3 + .../blocks-mesh-complete/domain_ids.dat | 3 + .../mesh-complete.exterior.vtp | 3 + .../mesh-complete.mesh.vtu | 3 + .../blocks-mesh-complete/mesh-surfaces/X0.vtp | 3 + .../blocks-mesh-complete/mesh-surfaces/X1.vtp | 3 + .../blocks-mesh-complete/mesh-surfaces/Y0.vtp | 3 + .../blocks-mesh-complete/mesh-surfaces/Y1.vtp | 3 + .../blocks-mesh-complete/mesh-surfaces/Z0.vtp | 3 + .../blocks-mesh-complete/mesh-surfaces/Z1.vtp | 3 + .../blocks-mesh-complete/walls_combined.vtp | 3 + .../struct/uniaxial_block_stretch/load.dat | 3 + .../mesh/mesh-complete.mesh.vtu | 3 + .../mesh/mesh-surfaces/X0.vtp | 3 + .../mesh/mesh-surfaces/X1.vtp | 3 + .../mesh/mesh-surfaces/Y0.vtp | 3 + .../mesh/mesh-surfaces/Y1.vtp | 3 + .../mesh/mesh-surfaces/Z0.vtp | 3 + .../mesh/mesh-surfaces/Z1.vtp | 3 + .../struct/uniaxial_block_stretch/solver.xml | 126 ++++++++++++++++++ 24 files changed, 220 insertions(+), 8 deletions(-) create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/block_domains.vtu create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/domain_ids.dat create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.exterior.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.mesh.vtu create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X0.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X1.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y0.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y1.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z0.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z1.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/walls_combined.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/load.dat create mode 100644 tests/cases/struct/uniaxial_block_stretch/mesh/mesh-complete.mesh.vtu create mode 100644 tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/X0.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/X1.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Y0.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Y1.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Z0.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Z1.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/solver.xml diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 5743a13a3..2f1d214b4 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -342,9 +342,11 @@ enum class OutputNameType outGrp_divV = 521, outGrp_Visc = 522, outGrp_fS = 523, - outGrp_C = 524, + outGrp_C = 524, outGrp_I1 = 525, outGrp_ionicState = 526, + outGrp_I4f = 527, + outGrp_I4fRate = 528, out_velocity = 599, out_pressure = 598, @@ -373,7 +375,9 @@ enum class OutputNameType out_viscosity = 575, out_fibStrn = 574, out_CGstrain = 573, - out_CGInv1 = 572 + out_CGInv1 = 572, + out_fibStretch = 571, + out_fibStretchRate = 570 }; /// @brief Simulation output file types. diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index 950a5939f..b797c4a02 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -556,12 +556,13 @@ SetEquationPropertiesMapType set_equation_props = { outPuts = {OutputNameType::out_displacement, OutputNameType::out_stress, OutputNameType::out_cauchy, OutputNameType::out_strain}; //simulation->com_mod.pstEq = true; } else { - nDOP = {12,2,0,0}; - outPuts = { + nDOP = {14,2,0,0}; + outPuts = { OutputNameType::out_displacement, OutputNameType::out_mises, OutputNameType::out_stress, OutputNameType::out_cauchy, OutputNameType::out_strain, OutputNameType::out_jacobian, OutputNameType::out_defGrad, OutputNameType::out_integ, OutputNameType::out_fibDir, - OutputNameType::out_fibAlign, OutputNameType::out_velocity, OutputNameType::out_acceleration + OutputNameType::out_fibAlign, OutputNameType::out_velocity, OutputNameType::out_acceleration, + OutputNameType::out_fibStretch, OutputNameType::out_fibStretchRate }; } @@ -596,7 +597,7 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); - nDOP = {14, 2, 0, 0}; + nDOP = {16, 2, 0, 0}; outPuts = { OutputNameType::out_displacement, OutputNameType::out_mises, @@ -611,7 +612,9 @@ SetEquationPropertiesMapType set_equation_props = { OutputNameType::out_velocity, OutputNameType::out_pressure, OutputNameType::out_acceleration, - OutputNameType::out_divergence + OutputNameType::out_divergence, + OutputNameType::out_fibStretch, + OutputNameType::out_fibStretchRate }; // Set solver parameters. diff --git a/Code/Source/solver/set_output_props.h b/Code/Source/solver/set_output_props.h index e9b314494..e34dac4c4 100644 --- a/Code/Source/solver/set_output_props.h +++ b/Code/Source/solver/set_output_props.h @@ -38,6 +38,8 @@ std::map output_props_map = {OutputNameType::out_fibAlign, std::make_tuple(OutputNameType::outGrp_fA, 0, 1, "Fiber_alignment") }, {OutputNameType::out_fibDir, std::make_tuple(OutputNameType::outGrp_fN, 0, nsd, "Fiber_direction") }, {OutputNameType::out_fibStrn, std::make_tuple(OutputNameType::outGrp_fS, 0, 1, "Fiber_shortening") }, + {OutputNameType::out_fibStretch, std::make_tuple(OutputNameType::outGrp_I4f, 0, 1, "Fiber_stretch") }, + {OutputNameType::out_fibStretchRate, std::make_tuple(OutputNameType::outGrp_I4fRate, 0, 1, "Fiber_stretch_rate") }, {OutputNameType::out_heatFlux, std::make_tuple(OutputNameType::outGrp_hFlx, 0, nsd, "Heat_flux") }, {OutputNameType::out_integ, std::make_tuple(OutputNameType::outGrp_I, 0, 1, nsd == 2 ? "Area" : "Volume") }, diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 818549b3c..6da073e81 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1293,7 +1293,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b break; case OutputNameType::outGrp_divV: - tmpV.resize(l,msh.nNo); + tmpV.resize(l,msh.nNo); post::div_post(simulation, msh, tmpV, solutions, iEq); for (int a = 0; a < msh.nNo; a++) { d[iM].x(is,a) = tmpV(0,a); @@ -1301,6 +1301,26 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b tmpV.resize(consts::maxNSD,msh.nNo); break; + case OutputNameType::outGrp_I4f: { + Vector res(msh.nNo); + if (msh.nFn != 0) { + post::fib_stretch(simulation, iEq, msh, solutions, res); + } + for (int a = 0; a < msh.nNo; a++) { + d[iM].x(is,a) = res(a); + } + } break; + + case OutputNameType::outGrp_I4fRate: { + Vector res(msh.nNo); + if (msh.nFn != 0) { + post::fib_stretch_rate(simulation, iEq, msh, solutions, res); + } + for (int a = 0; a < msh.nNo; a++) { + d[iM].x(is,a) = res(a); + } + } break; + default: throw std::runtime_error("Undefined output"); break; diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/block_domains.vtu b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/block_domains.vtu new file mode 100644 index 000000000..c6cb6d1e2 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/block_domains.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:a0cb409f170784be249c44f4bcfe79fc2091896fbb168df08ad8b7794c9d83e4 +size 532192 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/domain_ids.dat b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/domain_ids.dat new file mode 100644 index 000000000..008f628fd --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/domain_ids.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9b48b35d7e8ab963bbe30e79cb31fdbcf53378f04efdd73f76d849e42d167d67 +size 56292 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.exterior.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.exterior.vtp new file mode 100644 index 000000000..0715a1204 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:a068c4c7059084eb73e847b160232272caeb6757c4df14e9352ab937e9ad40bd +size 101403 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.mesh.vtu b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.mesh.vtu new file mode 100644 index 000000000..fef5731de --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f6acc26ca3faadd43dfc6b04362a30576fe649b093675c92035f34eab0f4a5bd +size 393738 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X0.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X0.vtp new file mode 100644 index 000000000..49a8a698e --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1419ac505ac92e9a26ca3f2866b658f1fe730ff838e767e93b8fc0c83688e1b3 +size 10393 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X1.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X1.vtp new file mode 100644 index 000000000..93f8b4450 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f4aa3f141aff31b65b4db9f729a1320ca20d2660413ee414358394a104187bf1 +size 10064 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y0.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y0.vtp new file mode 100644 index 000000000..827761b5c --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6d10f855cf7d09492bca423bdde8863d8fce7651ed77fabd9888619ab2396648 +size 44146 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y1.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y1.vtp new file mode 100644 index 000000000..5334fb594 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bea8c981cc4627a80dcd1a45509a6ea7c213cfa2c7673381feade49a89aa9717 +size 40392 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z0.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z0.vtp new file mode 100644 index 000000000..ce46f8870 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:626f66a867ac3118c34f15b1ac82699df775be01d260703ea9804e25f6f16c09 +size 4846 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z1.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z1.vtp new file mode 100644 index 000000000..bd80afca8 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d3a80a2c57b1b8e362d7b0111cf458500b1043e282174a3c34377f4747b2d7a1 +size 5062 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/walls_combined.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/walls_combined.vtp new file mode 100644 index 000000000..1cf403c5c --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/walls_combined.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ccd0154deef7668e4243672246ce03c5cf856d04f89d73211e2acd3e9b300ed9 +size 94430 diff --git a/tests/cases/struct/uniaxial_block_stretch/load.dat b/tests/cases/struct/uniaxial_block_stretch/load.dat new file mode 100644 index 000000000..b40bc41c8 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/load.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:09f6500d9c93686d7c623a47890b244cadea9bae9a6546bd27004a6259b5ffb0 +size 22 diff --git a/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-complete.mesh.vtu new file mode 100644 index 000000000..f09242eb6 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d74418f571657908b5c4e73f7cc5ce040fd300ac444b1c893ab23001bb6e09ae +size 7499 diff --git a/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/X0.vtp b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/X0.vtp new file mode 100644 index 000000000..904401ad3 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f71480fecee6f177719ab68b14bb2b61bbbce8fc745c87576fd51f3fb0b3ef8 +size 3655 diff --git a/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/X1.vtp b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/X1.vtp new file mode 100644 index 000000000..20b2d0e58 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd46bcf90705c0e062aa658846aacb9091334a56fffac81715631628cd1f097f +size 3669 diff --git a/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Y0.vtp b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Y0.vtp new file mode 100644 index 000000000..19c7f0081 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Y0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ef44a9a441563692f58baccf576b3c10e6df4277046d1054db5ab821d498bf0 +size 4795 diff --git a/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Y1.vtp b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Y1.vtp new file mode 100644 index 000000000..7ac40cb0a --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Y1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4d2ad3c14f29583d9d2015b0c0ba681b66c0d2380e4b8851ca30a2cb81773183 +size 4629 diff --git a/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Z0.vtp b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Z0.vtp new file mode 100644 index 000000000..49112c41e --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Z0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6ee78716770c3ed486e5e81d32f90226e215bc095b5cda4f5c708218b84ee407 +size 3188 diff --git a/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Z1.vtp b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Z1.vtp new file mode 100644 index 000000000..22122b6fb --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/mesh/mesh-surfaces/Z1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c2bdc7813853e641c5a71c20f1cea215c9c0215523f51272c5650077614dcb29 +size 3245 diff --git a/tests/cases/struct/uniaxial_block_stretch/solver.xml b/tests/cases/struct/uniaxial_block_stretch/solver.xml new file mode 100644 index 000000000..21d436dc2 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/solver.xml @@ -0,0 +1,126 @@ + + + + + 0 + 3 + 100 + 0.01 + 0.50 + STOP_SIM + + 1 + result + 10 + 1 + + 100 + 0 + + 1 + 0 + 0 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/X0.vtp + + + + mesh/mesh-surfaces/X1.vtp + + + + mesh/mesh-surfaces/Y0.vtp + + + + mesh/mesh-surfaces/Y1.vtp + + + + mesh/mesh-surfaces/Z0.vtp + + + + mesh/mesh-surfaces/Z1.vtp + + + (0.0, 0.0, 1.0) + (1.0, 0.0, 0.0) + + 1.0 + + + + + + + true + 1 + 10 + 1e-12 + + 240e6 + 0.45 + + + ST91 + + + true + true + true + true + true + true + true + + + + + fsils + + 1e-12 + 600 + + + + Dir + 0.0 + (1, 0, 0) + + + + Dir + 0.0 + (0, 1, 0) + + + + Dir + 0.0 + (0, 0, 1) + false + + + + Dir + Unsteady + load.dat + true + (0, 0, 1) + false + + + + + + + From 23617fd9f8aaeb4214c4a33a642ded361d8886c0 Mon Sep 17 00:00:00 2001 From: KB Date: Mon, 18 May 2026 16:23:02 -0700 Subject: [PATCH 04/15] changed fib stretch rate calc to FD method --- Code/Source/solver/post.cpp | 198 ++++-------------- .../struct/uniaxial_block_stretch/solver.xml | 1 + 2 files changed, 38 insertions(+), 161 deletions(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 5152e66be..c2fd88f8a 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -740,45 +740,37 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra } -/// @brief Compute fiber stretch based on 4th invariant: I_{4,f} +/// @brief L2-project nodal fiber stretch λ = sqrt(I_{4,f}) from a given displacement field. // -void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) +static void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, + const Array& lD, Vector& res) { using namespace consts; auto& com_mod = simulation->com_mod; - auto& cm = com_mod.cm; - auto& cm_mod = simulation->cm_mod; - const auto& lD = solutions.old.get_displacement(); auto& eq = com_mod.eq[iEq]; - int nsd = com_mod.nsd; + int nsd = com_mod.nsd; int tnNo = com_mod.tnNo; int tDof = com_mod.tDof; - - // [NOTE] Setting gobal variable 'dof'. - com_mod.dof = eq.dof; - int eNoN = lM.eNoN; int i = eq.s; int j = i + 1; int k = j + 1; - Vector sA(tnNo); - Vector sF(tnNo); - Array xl(nsd,eNoN); + Vector sA(tnNo); + Vector sF(tnNo); + Array xl(nsd,eNoN); Array dl(tDof,eNoN); - Array Nx(nsd,eNoN); - Vector N(eNoN); + Array Nx(nsd,eNoN); for (int e = 0; e < lM.nEl; e++) { - int cDmn = all_fun::domain(com_mod, lM, iEq, e); - auto cPhys = eq.dmn[cDmn].phys; + all_fun::domain(com_mod, lM, iEq, e); if (lM.eType == ElementType::NRB) { //CALL NRBNNX(lM, e) } - for (int a = 0; a < eNoN; a++) { + for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); xl.set_col(a, com_mod.x.col(Ac)); dl.set_col(a, lD.col(Ac)); @@ -797,8 +789,8 @@ void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const auto N = lM.N.col(g); // Compute Deformation Gradient: F = I + grad(u) - F = mat_fun::mat_id(nsd); - for (int a = 0; a < eNoN; a++) { + F = mat_fun::mat_id(nsd); + for (int a = 0; a < eNoN; a++) { if (nsd == 3) { F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); @@ -819,13 +811,13 @@ void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const // Compute fiber stretch based on 4th invariant: I_{4,f} = F.fN.F.fN auto fl = mat_fun::mat_mul(F, lM.fN.rows(0,nsd-1,e)); - double I4f = utils::norm(fl); + double lambda = sqrt(utils::norm(fl)); - // L2 projection of I4f from integration points to nodes - for (int a = 0; a < eNoN; a++) { + // L2 projection from integration points to nodes + for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); sA(Ac) = sA(Ac) + w*N(a); - sF(Ac) = sF(Ac) + w*N(a)*I4f; + sF(Ac) = sF(Ac) + w*N(a)*lambda; } } } @@ -841,160 +833,44 @@ void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const res(a) = res(a) + sF(Ac) / sA(Ac); } } - } -/// @brief Compute fiber stretch rate dλ/dt = (1/2λ) fN^T (dC/dt) fN + +/// @brief Compute fiber stretch based on 4th invariant: λ = sqrt(I_{4,f}) // -void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) +void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) { - using namespace consts; - auto& com_mod = simulation->com_mod; - auto& cm = com_mod.cm; - auto& cm_mod = simulation->cm_mod; - const auto& lD = solutions.old.get_displacement(); - const auto& lY = solutions.current.get_velocity(); auto& eq = com_mod.eq[iEq]; - - int nsd = com_mod.nsd; - int tnNo = com_mod.tnNo; - int tDof = com_mod.tDof; - - // [NOTE] Setting gobal variable 'dof'. com_mod.dof = eq.dof; + compute_fib_stretch(simulation, iEq, lM, solutions.current.get_displacement(), res); +} - int eNoN = lM.eNoN; - int i = eq.s; - int j = i + 1; - int k = j + 1; - - Vector sA(tnNo); - Vector sF(tnNo); - Array xl(nsd,eNoN); - Array dl(tDof,eNoN); - Array yl(tDof,eNoN); - Array Nx(nsd,eNoN); - - for (int e = 0; e < lM.nEl; e++) { - int cDmn = all_fun::domain(com_mod, lM, iEq, e); - auto cPhys = eq.dmn[cDmn].phys; - if (cPhys != EquationType::phys_struct && cPhys != EquationType::phys_ustruct) { - continue; - } - if (lM.eType == ElementType::NRB) { - //CALL NRBNNX(lM, e) - } - - for (int a = 0; a < eNoN; a++) { - int Ac = lM.IEN(a,e); - xl.set_col(a, com_mod.x.col(Ac)); - dl.set_col(a, lD.col(Ac)); - for (int i = 0; i < tDof; i++) { - yl(i,a) = lY(i,Ac); - } - } - - for (int g = 0; g < lM.nG; g++) { - double Jac = 0.0; - Array F(nsd,nsd); - Array vx(nsd,nsd); - Array ksix(nsd,nsd); - - if (g == 0 || !lM.lShpF) { - auto Nxi = lM.Nx.slice(g); - nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, ksix); - } - - double w = lM.w(g)*Jac; - auto N = lM.N.col(g); - F = mat_fun::mat_id(nsd); - - for (int a = 0; a < eNoN; a++) { - if (nsd == 3) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(j,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(k,a); - - vx(0,0) = vx(0,0) + Nx(0,a)*yl(i,a); - vx(0,1) = vx(0,1) + Nx(1,a)*yl(i,a); - vx(0,2) = vx(0,2) + Nx(2,a)*yl(i,a); - vx(1,0) = vx(1,0) + Nx(0,a)*yl(j,a); - vx(1,1) = vx(1,1) + Nx(1,a)*yl(j,a); - vx(1,2) = vx(1,2) + Nx(2,a)*yl(j,a); - vx(2,0) = vx(2,0) + Nx(0,a)*yl(k,a); - vx(2,1) = vx(2,1) + Nx(1,a)*yl(k,a); - vx(2,2) = vx(2,2) + Nx(2,a)*yl(k,a); - } else { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - - vx(0,0) = vx(0,0) + Nx(0,a)*yl(i,a); - vx(0,1) = vx(0,1) + Nx(1,a)*yl(i,a); - vx(1,0) = vx(1,0) + Nx(0,a)*yl(j,a); - vx(1,1) = vx(1,1) + Nx(1,a)*yl(j,a); - } - } - auto fl = mat_fun::mat_mul(F, lM.fN.rows(0,nsd-1,e)); - double I4f = utils::norm(fl); - double lambda = sqrt(I4f); - if (lambda < 1e-12) { - lambda = 1.0; - } - - // dC/dt = (vx)^T F + F^T vx - Array Vx_T = mat_fun::transpose(vx); - Array F_T = mat_fun::transpose(F); - Array dC_dt(nsd, nsd); - auto Vx_T_F = mat_fun::mat_mul(Vx_T, F); - auto F_T_Vx = mat_fun::mat_mul(F_T, vx); - for (int i = 0; i < nsd; i++) { - for (int j = 0; j < nsd; j++) { - dC_dt(i,j) = Vx_T_F(i,j) + F_T_Vx(i,j); - } - } +/// @brief Compute fiber stretch rate dλ/dt via backward finite difference. +// +void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) +{ + auto& com_mod = simulation->com_mod; + auto& eq = com_mod.eq[iEq]; + com_mod.dof = eq.dof; - auto fN = lM.fN.rows(0,nsd-1,e); - double fT_dC_dt_f = 0.0; - for (int i = 0; i < nsd; i++) { - for (int j = 0; j < nsd; j++) { - fT_dC_dt_f += fN(i) * dC_dt(i,j) * fN(j); - } - } + const double dt = com_mod.dt; + int nNo = lM.nNo; - double dlambda_dt = fT_dC_dt_f / (2.0 * lambda); + Vector lambda_curr(nNo); + Vector lambda_old(nNo); - for (int a = 0; a < eNoN; a++) { - int Ac = lM.IEN(a,e); - sA(Ac) = sA(Ac) + w*N(a); - sF(Ac) = sF(Ac) + w*N(a)*dlambda_dt; - } - } - } - - all_fun::commu(com_mod, sF); - all_fun::commu(com_mod, sA); + compute_fib_stretch(simulation, iEq, lM, solutions.current.get_displacement(), lambda_curr); + compute_fib_stretch(simulation, iEq, lM, solutions.old.get_displacement(), lambda_old); res = 0.0; - - for (int a = 0; a < lM.nNo; a++) { - int Ac = lM.gN(a); - if (!utils::is_zero(sA(Ac))) { - res(a) = res(a) + sF(Ac) / sA(Ac); - } + for (int a = 0; a < nNo; a++) { + res(a) = (lambda_curr(a) - lambda_old(a)) / dt; } - } + void post(Simulation* simulation, const mshType& lM, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp, const int iEq) { diff --git a/tests/cases/struct/uniaxial_block_stretch/solver.xml b/tests/cases/struct/uniaxial_block_stretch/solver.xml index 21d436dc2..a45270a8d 100644 --- a/tests/cases/struct/uniaxial_block_stretch/solver.xml +++ b/tests/cases/struct/uniaxial_block_stretch/solver.xml @@ -75,6 +75,7 @@ true + true true true true From a626e668ab9250d3bd1b84a449c3398c1197e65a Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 19 May 2026 10:14:27 -0700 Subject: [PATCH 05/15] added reference solution for test case --- .../blocks-mesh-complete/block_domains.vtu | 3 --- .../blocks-mesh-complete/domain_ids.dat | 3 --- .../blocks-mesh-complete/mesh-complete.exterior.vtp | 3 --- .../blocks-mesh-complete/mesh-complete.mesh.vtu | 3 --- .../blocks-mesh-complete/mesh-surfaces/X0.vtp | 3 --- .../blocks-mesh-complete/mesh-surfaces/X1.vtp | 3 --- .../blocks-mesh-complete/mesh-surfaces/Y0.vtp | 3 --- .../blocks-mesh-complete/mesh-surfaces/Y1.vtp | 3 --- .../blocks-mesh-complete/mesh-surfaces/Z0.vtp | 3 --- .../blocks-mesh-complete/mesh-surfaces/Z1.vtp | 3 --- .../blocks-mesh-complete/walls_combined.vtp | 3 --- tests/cases/struct/uniaxial_block_stretch/result_100.vtu | 3 +++ tests/conftest.py | 2 ++ tests/test_struct.py | 4 ++++ 14 files changed, 9 insertions(+), 33 deletions(-) delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/block_domains.vtu delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/domain_ids.dat delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.exterior.vtp delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.mesh.vtu delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X0.vtp delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X1.vtp delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y0.vtp delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y1.vtp delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z0.vtp delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z1.vtp delete mode 100644 tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/walls_combined.vtp create mode 100644 tests/cases/struct/uniaxial_block_stretch/result_100.vtu diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/block_domains.vtu b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/block_domains.vtu deleted file mode 100644 index c6cb6d1e2..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/block_domains.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:a0cb409f170784be249c44f4bcfe79fc2091896fbb168df08ad8b7794c9d83e4 -size 532192 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/domain_ids.dat b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/domain_ids.dat deleted file mode 100644 index 008f628fd..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/domain_ids.dat +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:9b48b35d7e8ab963bbe30e79cb31fdbcf53378f04efdd73f76d849e42d167d67 -size 56292 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.exterior.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.exterior.vtp deleted file mode 100644 index 0715a1204..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.exterior.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:a068c4c7059084eb73e847b160232272caeb6757c4df14e9352ab937e9ad40bd -size 101403 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.mesh.vtu b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.mesh.vtu deleted file mode 100644 index fef5731de..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-complete.mesh.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:f6acc26ca3faadd43dfc6b04362a30576fe649b093675c92035f34eab0f4a5bd -size 393738 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X0.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X0.vtp deleted file mode 100644 index 49a8a698e..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X0.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:1419ac505ac92e9a26ca3f2866b658f1fe730ff838e767e93b8fc0c83688e1b3 -size 10393 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X1.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X1.vtp deleted file mode 100644 index 93f8b4450..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/X1.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:f4aa3f141aff31b65b4db9f729a1320ca20d2660413ee414358394a104187bf1 -size 10064 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y0.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y0.vtp deleted file mode 100644 index 827761b5c..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y0.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6d10f855cf7d09492bca423bdde8863d8fce7651ed77fabd9888619ab2396648 -size 44146 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y1.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y1.vtp deleted file mode 100644 index 5334fb594..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Y1.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:bea8c981cc4627a80dcd1a45509a6ea7c213cfa2c7673381feade49a89aa9717 -size 40392 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z0.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z0.vtp deleted file mode 100644 index ce46f8870..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z0.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:626f66a867ac3118c34f15b1ac82699df775be01d260703ea9804e25f6f16c09 -size 4846 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z1.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z1.vtp deleted file mode 100644 index bd80afca8..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/mesh-surfaces/Z1.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:d3a80a2c57b1b8e362d7b0111cf458500b1043e282174a3c34377f4747b2d7a1 -size 5062 diff --git a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/walls_combined.vtp b/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/walls_combined.vtp deleted file mode 100644 index 1cf403c5c..000000000 --- a/tests/cases/struct/uniaxial_block_stretch/blocks-mesh-complete/walls_combined.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ccd0154deef7668e4243672246ce03c5cf856d04f89d73211e2acd3e9b300ed9 -size 94430 diff --git a/tests/cases/struct/uniaxial_block_stretch/result_100.vtu b/tests/cases/struct/uniaxial_block_stretch/result_100.vtu new file mode 100644 index 000000000..1f3632205 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/result_100.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1f954966510318f36048b1b7c44fd9877fa398a497d769c25e8479cad12ea3c9 +size 9844 diff --git a/tests/conftest.py b/tests/conftest.py index 7c4a36cab..d2f884c96 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -33,6 +33,8 @@ "VonMises_stress": 1.0e-3, "Vorticity": 1.0e-7, "WSS": 1.0e-8, + "Fiber_stretch": 1.0e-10, + "Fiber_stretch_rate": 1.0e-10, } # Number of processors to test diff --git a/tests/test_struct.py b/tests/test_struct.py index 10188253a..0fecd6425 100644 --- a/tests/test_struct.py +++ b/tests/test_struct.py @@ -110,3 +110,7 @@ def test_tensile_adventitia_Newtonian_viscosity(n_proc): def test_tensile_adventitia_Potential_viscosity(n_proc): test_folder = "tensile_adventitia_Potential_viscosity" run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) + +def test_uniaxial_block_stretch(n_proc): + test_folder = "uniaxial_block_stretch" + run_with_reference(base_folder, test_folder, ["Fiber_stretch", "Fiber_stretch_rate"], n_proc, t_max=100) \ No newline at end of file From 41c82622314215d3838343c9616436a75ec8cc73 Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 19 May 2026 13:35:14 -0700 Subject: [PATCH 06/15] added documentation to test case --- Code/Source/solver/post.cpp | 2 +- .../struct/uniaxial_block_stretch/BlockStretch.gif | 3 +++ tests/cases/struct/uniaxial_block_stretch/README.md | 11 +++++++++++ .../struct/uniaxial_block_stretch/fiberstretch.png | 3 +++ 4 files changed, 18 insertions(+), 1 deletion(-) create mode 100644 tests/cases/struct/uniaxial_block_stretch/BlockStretch.gif create mode 100644 tests/cases/struct/uniaxial_block_stretch/README.md create mode 100644 tests/cases/struct/uniaxial_block_stretch/fiberstretch.png diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index c2fd88f8a..199fdfd50 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -750,7 +750,7 @@ static void compute_fib_stretch(Simulation* simulation, const int iEq, const msh auto& com_mod = simulation->com_mod; auto& eq = com_mod.eq[iEq]; - int nsd = com_mod.nsd; + int nsd = com_mod.nsd; int tnNo = com_mod.tnNo; int tDof = com_mod.tDof; int eNoN = lM.eNoN; diff --git a/tests/cases/struct/uniaxial_block_stretch/BlockStretch.gif b/tests/cases/struct/uniaxial_block_stretch/BlockStretch.gif new file mode 100644 index 000000000..73f7c4717 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/BlockStretch.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5db6a5c985750cced6062e9903dc89ec0bbc57aa2165a806b3c8c36f579775e1 +size 81267 diff --git a/tests/cases/struct/uniaxial_block_stretch/README.md b/tests/cases/struct/uniaxial_block_stretch/README.md new file mode 100644 index 000000000..8f7eb2743 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/README.md @@ -0,0 +1,11 @@ +This test case simulates uniaxial stretch on a slab of material. The reference solution contains the fiber stretch and fiber stretch rate fields that were obtained analytically by computing $\lambda = 1 + \frac{\Delta u}{L_0}$ and $\frac{d \lambda}{dt} = \frac{1}{L_0}\frac{\Delta u}{\Delta T}$. Note that $L_0$ is the initial length, $\Delta u$ is the displacement, and $\Delta T$ is the time during which the deformation is applied. + +The plot below shows the analytical and computed fiber stretch. +

+ +

+ +The resulting deformation and the fiber stretch values are shown in the video below: +

+ +

diff --git a/tests/cases/struct/uniaxial_block_stretch/fiberstretch.png b/tests/cases/struct/uniaxial_block_stretch/fiberstretch.png new file mode 100644 index 000000000..af52f3af3 --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/fiberstretch.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6858da7b986421b77b92ccf3972d8a325704a8b0f13f0e3084b047792f99112c +size 331355 From ea047984763578a9f08c5b6435e3b92afda66c51 Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 19 May 2026 13:37:22 -0700 Subject: [PATCH 07/15] minor formatting --- tests/cases/struct/uniaxial_block_stretch/README.md | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/tests/cases/struct/uniaxial_block_stretch/README.md b/tests/cases/struct/uniaxial_block_stretch/README.md index 8f7eb2743..0bbea78a5 100644 --- a/tests/cases/struct/uniaxial_block_stretch/README.md +++ b/tests/cases/struct/uniaxial_block_stretch/README.md @@ -1,4 +1,10 @@ -This test case simulates uniaxial stretch on a slab of material. The reference solution contains the fiber stretch and fiber stretch rate fields that were obtained analytically by computing $\lambda = 1 + \frac{\Delta u}{L_0}$ and $\frac{d \lambda}{dt} = \frac{1}{L_0}\frac{\Delta u}{\Delta T}$. Note that $L_0$ is the initial length, $\Delta u$ is the displacement, and $\Delta T$ is the time during which the deformation is applied. +This test case simulates uniaxial stretch on a slab of material. The reference solution contains the fiber stretch and fiber stretch rate fields that were obtained analytically. + +The fiber stretch is computed as $\lambda = 1 + \frac{\Delta u}{L_0}$. + +The fiber stretch rate is given by $\frac{d \lambda}{dt} = \frac{1}{L_0}\frac{\Delta u}{\Delta T}$. + +Note that $L_0$ is the initial length, $\Delta u$ is the displacement, and $\Delta T$ is the time during which the deformation is applied. The plot below shows the analytical and computed fiber stretch.

From a70e028d30ca5ab97cac778f07186568cc419f00 Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 19 May 2026 13:38:59 -0700 Subject: [PATCH 08/15] latex formatting --- tests/cases/struct/uniaxial_block_stretch/README.md | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/tests/cases/struct/uniaxial_block_stretch/README.md b/tests/cases/struct/uniaxial_block_stretch/README.md index 0bbea78a5..a6764c167 100644 --- a/tests/cases/struct/uniaxial_block_stretch/README.md +++ b/tests/cases/struct/uniaxial_block_stretch/README.md @@ -1,10 +1,4 @@ -This test case simulates uniaxial stretch on a slab of material. The reference solution contains the fiber stretch and fiber stretch rate fields that were obtained analytically. - -The fiber stretch is computed as $\lambda = 1 + \frac{\Delta u}{L_0}$. - -The fiber stretch rate is given by $\frac{d \lambda}{dt} = \frac{1}{L_0}\frac{\Delta u}{\Delta T}$. - -Note that $L_0$ is the initial length, $\Delta u$ is the displacement, and $\Delta T$ is the time during which the deformation is applied. +This test case simulates uniaxial stretch on a slab of material. The reference solution contains the fiber stretch and fiber stretch rate fields that were obtained analytically. The fiber stretch is computed as $\small \lambda = 1 + \frac{\Delta u}{L_0}$ and the fiber stretch rate is given by $\small \frac{d \lambda}{dt} = \frac{1}{L_0}\frac{\Delta u}{\Delta T}$. Note that $L_0$ is the initial length, $\Delta u$ is the displacement, and $\Delta T$ is the time during which the deformation is applied. The plot below shows the analytical and computed fiber stretch.

From cc476400a362e6e8f0f0d7611cd92db2edeb73c0 Mon Sep 17 00:00:00 2001 From: KB Date: Wed, 20 May 2026 08:50:50 -0700 Subject: [PATCH 09/15] explicitly expose compute_fib_stretch function --- Code/Source/solver/post.cpp | 2 +- Code/Source/solver/post.h | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 199fdfd50..7bcce5e1c 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -742,7 +742,7 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra /// @brief L2-project nodal fiber stretch λ = sqrt(I_{4,f}) from a given displacement field. // -static void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, +void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const Array& lD, Vector& res) { using namespace consts; diff --git a/Code/Source/solver/post.h b/Code/Source/solver/post.h index ee9616f12..bfa7fbe81 100644 --- a/Code/Source/solver/post.h +++ b/Code/Source/solver/post.h @@ -22,6 +22,8 @@ void fib_algn_post(Simulation* simulation, const mshType& lM, Array& res void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Array& res, const SolutionStates& solutions, const int iEq); +void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const Array& lD, Vector& res); + void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); From 2c9164a5bfce0d0c5ba9164e33a27ba6e396b5a4 Mon Sep 17 00:00:00 2001 From: KB Date: Thu, 21 May 2026 09:15:33 -0700 Subject: [PATCH 10/15] renamed I4f to fibStretch --- Code/Source/solver/consts.h | 5 ++--- Code/Source/solver/post.cpp | 2 +- Code/Source/solver/set_output_props.h | 5 ++--- Code/Source/solver/vtk_xml.cpp | 5 ++--- tests/cases/struct/uniaxial_block_stretch/README.md | 5 +++++ .../cases/struct/uniaxial_block_stretch/fiberstretchrate.png | 3 +++ 6 files changed, 15 insertions(+), 10 deletions(-) create mode 100644 tests/cases/struct/uniaxial_block_stretch/fiberstretchrate.png diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 2f1d214b4..40d81a7d3 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -345,8 +345,8 @@ enum class OutputNameType outGrp_C = 524, outGrp_I1 = 525, outGrp_ionicState = 526, - outGrp_I4f = 527, - outGrp_I4fRate = 528, + outGrp_fibStretch = 527, + outGrp_fibStretchRate = 528, out_velocity = 599, out_pressure = 598, @@ -501,4 +501,3 @@ enum class LinearAlgebraType { }; #endif - diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 7bcce5e1c..d1ae03fba 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -740,7 +740,7 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra } -/// @brief L2-project nodal fiber stretch λ = sqrt(I_{4,f}) from a given displacement field. +/// @brief Compute fiber stretch λ = sqrt(I_{4,f}) from a given displacement field. // void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const Array& lD, Vector& res) diff --git a/Code/Source/solver/set_output_props.h b/Code/Source/solver/set_output_props.h index e34dac4c4..59660ce6f 100644 --- a/Code/Source/solver/set_output_props.h +++ b/Code/Source/solver/set_output_props.h @@ -38,8 +38,8 @@ std::map output_props_map = {OutputNameType::out_fibAlign, std::make_tuple(OutputNameType::outGrp_fA, 0, 1, "Fiber_alignment") }, {OutputNameType::out_fibDir, std::make_tuple(OutputNameType::outGrp_fN, 0, nsd, "Fiber_direction") }, {OutputNameType::out_fibStrn, std::make_tuple(OutputNameType::outGrp_fS, 0, 1, "Fiber_shortening") }, - {OutputNameType::out_fibStretch, std::make_tuple(OutputNameType::outGrp_I4f, 0, 1, "Fiber_stretch") }, - {OutputNameType::out_fibStretchRate, std::make_tuple(OutputNameType::outGrp_I4fRate, 0, 1, "Fiber_stretch_rate") }, + {OutputNameType::out_fibStretch, std::make_tuple(OutputNameType::outGrp_fibStretch, 0, 1, "Fiber_stretch") }, + {OutputNameType::out_fibStretchRate, std::make_tuple(OutputNameType::outGrp_fibStretchRate, 0, 1, "Fiber_stretch_rate") }, {OutputNameType::out_heatFlux, std::make_tuple(OutputNameType::outGrp_hFlx, 0, nsd, "Heat_flux") }, {OutputNameType::out_integ, std::make_tuple(OutputNameType::outGrp_I, 0, 1, nsd == 2 ? "Area" : "Volume") }, @@ -58,4 +58,3 @@ std::map output_props_map = {OutputNameType::out_vorticity, std::make_tuple(OutputNameType::outGrp_vort, 0, maxNSD, "Vorticity") }, {OutputNameType::out_WSS, std::make_tuple(OutputNameType::outGrp_WSS, 0, maxNSD, "WSS") } }; - diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 6da073e81..7d8bb499f 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1301,7 +1301,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b tmpV.resize(consts::maxNSD,msh.nNo); break; - case OutputNameType::outGrp_I4f: { + case OutputNameType::outGrp_fibStretch: { Vector res(msh.nNo); if (msh.nFn != 0) { post::fib_stretch(simulation, iEq, msh, solutions, res); @@ -1311,7 +1311,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b } } break; - case OutputNameType::outGrp_I4fRate: { + case OutputNameType::outGrp_fibStretchRate: { Vector res(msh.nNo); if (msh.nFn != 0) { post::fib_stretch_rate(simulation, iEq, msh, solutions, res); @@ -1553,4 +1553,3 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b } }; - diff --git a/tests/cases/struct/uniaxial_block_stretch/README.md b/tests/cases/struct/uniaxial_block_stretch/README.md index a6764c167..d0d34b84f 100644 --- a/tests/cases/struct/uniaxial_block_stretch/README.md +++ b/tests/cases/struct/uniaxial_block_stretch/README.md @@ -5,6 +5,11 @@ The plot below shows the analytical and computed fiber stretch.

+The analytical and computed fiber stretch rates are provided below. +

+ +

+ The resulting deformation and the fiber stretch values are shown in the video below:

diff --git a/tests/cases/struct/uniaxial_block_stretch/fiberstretchrate.png b/tests/cases/struct/uniaxial_block_stretch/fiberstretchrate.png new file mode 100644 index 000000000..5da043dfa --- /dev/null +++ b/tests/cases/struct/uniaxial_block_stretch/fiberstretchrate.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5c56a5f884a243e765f27d0522cc27d5ac1de167f4005be92ba8aa180f79e684 +size 309911 From 8a5f5cef384641b2d0039cb4713698c7318810d9 Mon Sep 17 00:00:00 2001 From: KB Date: Thu, 21 May 2026 12:43:17 -0700 Subject: [PATCH 11/15] Refactor deformation gradient in post.cpp --- Code/Source/solver/post.cpp | 128 ++++++++++-------------------------- 1 file changed, 34 insertions(+), 94 deletions(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index d1ae03fba..5a6c1b6ee 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -17,6 +17,34 @@ namespace post { + namespace { + + Array deformation_gradient(const Array& Nx, const Array& dl, + int nsd, int nNo, int eq_start) + { + auto F = mat_fun::mat_id(nsd); + for (int a = 0; a < nNo; a++) { + if (nsd == 3) { + F(0,0) = F(0,0) + Nx(0,a)*dl(eq_start,a); + F(0,1) = F(0,1) + Nx(1,a)*dl(eq_start,a); + F(0,2) = F(0,2) + Nx(2,a)*dl(eq_start,a); + F(1,0) = F(1,0) + Nx(0,a)*dl(eq_start+1,a); + F(1,1) = F(1,1) + Nx(1,a)*dl(eq_start+1,a); + F(1,2) = F(1,2) + Nx(2,a)*dl(eq_start+1,a); + F(2,0) = F(2,0) + Nx(0,a)*dl(eq_start+2,a); + F(2,1) = F(2,1) + Nx(1,a)*dl(eq_start+2,a); + F(2,2) = F(2,2) + Nx(2,a)*dl(eq_start+2,a); + } else { + F(0,0) = F(0,0) + Nx(0,a)*dl(eq_start,a); + F(0,1) = F(0,1) + Nx(1,a)*dl(eq_start,a); + F(1,0) = F(1,0) + Nx(0,a)*dl(eq_start+1,a); + F(1,1) = F(1,1) + Nx(1,a)*dl(eq_start+1,a); + } + } + return F; + } + } + void all_post(Simulation* simulation, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp, const int iEq) { @@ -443,17 +471,8 @@ void div_post(Simulation* simulation, const mshType& lM, Array& res, con vx(2,0) = vx(2,0) + Nx(0,a)*yl(k,a); vx(2,1) = vx(2,1) + Nx(1,a)*yl(k,a); vx(2,2) = vx(2,2) + Nx(2,a)*yl(k,a); - - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(j,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(k,a); } + F = deformation_gradient(Nx, dl, nsd, eNoN, i); auto Fi = mat_fun::mat_inv(F,3); @@ -468,12 +487,8 @@ void div_post(Simulation* simulation, const mshType& lM, Array& res, con vx(0,1) = vx(0,1) + Nx(1,a)*yl(i,a); vx(1,0) = vx(1,0) + Nx(0,a)*yl(j,a); vx(1,1) = vx(1,1) + Nx(1,a)*yl(j,a); - - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); } + F = deformation_gradient(Nx, dl, nsd, eNoN, i); auto Fi = mat_fun::mat_inv(F,2); VxFi(0) = vx(0,0)*Fi(0,0) + vx(0,1)*Fi(1,0); @@ -567,26 +582,7 @@ void fib_algn_post(Simulation* simulation, const mshType& lM, Array& res } double w = lM.w(g)*Jac; - auto F = mat_fun::mat_id(nsd); - - for (int a = 0; a < eNoN; a++) { - if (nsd == 3) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(j,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(k,a); - } else { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - } - } + auto F = deformation_gradient(Nx, dl, nsd, eNoN, i); for (int iFn = 0; iFn < 2; iFn++) { for (int i = 0; i < nsd; i++) { auto fN_col = fN.col(iFn); @@ -684,26 +680,7 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra double w = lM.w(g) * Jac; N = lM.N.col(g); - F = mat_fun::mat_id(nsd); - - for (int a = 0; a < eNoN; a++) { - if (nsd == 3) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(j,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(k,a); - } else { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - } - } + F = deformation_gradient(Nx, dl, nsd, eNoN, i); for (int iFn = 0; iFn < lM.nFn; iFn++) { for (int i = 0; i < nsd; i++) { @@ -789,25 +766,7 @@ void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& l auto N = lM.N.col(g); // Compute Deformation Gradient: F = I + grad(u) - F = mat_fun::mat_id(nsd); - for (int a = 0; a < eNoN; a++) { - if (nsd == 3) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(j,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(k,a); - } else { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - } - } + F = deformation_gradient(Nx, dl, nsd, eNoN, i); // Compute fiber stretch based on 4th invariant: I_{4,f} = F.fN.F.fN auto fl = mat_fun::mat_mul(F, lM.fN.rows(0,nsd-1,e)); @@ -1834,26 +1793,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array Je = Je + w; auto Im = mat_fun::mat_id(nsd); - auto F = Im; - - for (int a = 0; a < fs.eNoN; a++) { - if (nsd == 3) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(j,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(k,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(k,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(k,a); - } else { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(j,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(j,a); - } - } + auto F = deformation_gradient(Nx, dl, nsd, fs.eNoN, i); double detF = mat_fun::mat_det(F, nsd); From 081789624b4eefee32216c9c23f256026d64abe0 Mon Sep 17 00:00:00 2001 From: KB Date: Thu, 21 May 2026 14:18:09 -0700 Subject: [PATCH 12/15] code cleanup suggested by PR comments --- Code/Source/solver/post.cpp | 81 +++++++++++++++++-------------------- Code/Source/solver/post.h | 3 -- 2 files changed, 37 insertions(+), 47 deletions(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 5a6c1b6ee..7d4724931 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -3,6 +3,7 @@ #include "post.h" +#include "FE/Common/FEException.h" #include "all_fun.h" #include "fluid.h" #include "fs.h" @@ -43,6 +44,9 @@ namespace post { } return F; } + + void compute_fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, + const Array& lD, Vector& res); } void all_post(Simulation* simulation, Array& res, const SolutionStates& solutions, @@ -717,23 +721,48 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra } -/// @brief Compute fiber stretch λ = sqrt(I_{4,f}) from a given displacement field. +/// @brief Compute fiber stretch based on 4th invariant: λ = sqrt(I_{4,f}) +// +void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) +{ + compute_fib_stretch(simulation->com_mod, iEq, lM, solutions.current.get_displacement(), res); +} + +/// @brief Compute fiber stretch rate dλ/dt via backward finite difference. // -void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, +void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) +{ + auto& com_mod = simulation->com_mod; + const double dt = com_mod.dt; + int nNo = lM.nNo; + + if (res.size() != nNo) { + svmp::raise( + SVMP_HERE, + "[fib_stretch_rate] Expected res size " + std::to_string(nNo) + ", but got " + std::to_string(res.size()) + "."); + } + + Vector lambda_old(nNo); + + compute_fib_stretch(simulation->com_mod, iEq, lM, solutions.current.get_displacement(), res); + compute_fib_stretch(simulation->com_mod, iEq, lM, solutions.old.get_displacement(), lambda_old); + + res = (res - lambda_old) / dt; +} + +namespace { +/// @brief Compute fiber stretch lambda = sqrt(I_{4,f}) from a given displacement field. +void compute_fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, const Array& lD, Vector& res) { using namespace consts; - auto& com_mod = simulation->com_mod; auto& eq = com_mod.eq[iEq]; - int nsd = com_mod.nsd; int tnNo = com_mod.tnNo; int tDof = com_mod.tDof; int eNoN = lM.eNoN; int i = eq.s; - int j = i + 1; - int k = j + 1; Vector sA(tnNo); Vector sF(tnNo); @@ -762,9 +791,6 @@ void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& l nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, ksix); } - double w = lM.w(g)*Jac; - auto N = lM.N.col(g); - // Compute Deformation Gradient: F = I + grad(u) F = deformation_gradient(Nx, dl, nsd, eNoN, i); @@ -773,6 +799,8 @@ void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& l double lambda = sqrt(utils::norm(fl)); // L2 projection from integration points to nodes + double w = lM.w(g)*Jac; + auto N = lM.N.col(g); for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); sA(Ac) = sA(Ac) + w*N(a); @@ -793,43 +821,8 @@ void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& l } } } - - -/// @brief Compute fiber stretch based on 4th invariant: λ = sqrt(I_{4,f}) -// -void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) -{ - auto& com_mod = simulation->com_mod; - auto& eq = com_mod.eq[iEq]; - com_mod.dof = eq.dof; - compute_fib_stretch(simulation, iEq, lM, solutions.current.get_displacement(), res); } - -/// @brief Compute fiber stretch rate dλ/dt via backward finite difference. -// -void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) -{ - auto& com_mod = simulation->com_mod; - auto& eq = com_mod.eq[iEq]; - com_mod.dof = eq.dof; - - const double dt = com_mod.dt; - int nNo = lM.nNo; - - Vector lambda_curr(nNo); - Vector lambda_old(nNo); - - compute_fib_stretch(simulation, iEq, lM, solutions.current.get_displacement(), lambda_curr); - compute_fib_stretch(simulation, iEq, lM, solutions.old.get_displacement(), lambda_old); - - res = 0.0; - for (int a = 0; a < nNo; a++) { - res(a) = (lambda_curr(a) - lambda_old(a)) / dt; - } -} - - void post(Simulation* simulation, const mshType& lM, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp, const int iEq) { diff --git a/Code/Source/solver/post.h b/Code/Source/solver/post.h index bfa7fbe81..3eb6c73d4 100644 --- a/Code/Source/solver/post.h +++ b/Code/Source/solver/post.h @@ -22,8 +22,6 @@ void fib_algn_post(Simulation* simulation, const mshType& lM, Array& res void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Array& res, const SolutionStates& solutions, const int iEq); -void compute_fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const Array& lD, Vector& res); - void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); @@ -42,4 +40,3 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array }; #endif - From 66bdb9b3aa28a9bebfab153d51957d86df201bdc Mon Sep 17 00:00:00 2001 From: KB Date: Fri, 22 May 2026 09:10:38 -0700 Subject: [PATCH 13/15] simplified function calls and arguments --- Code/Source/solver/cep_ion.cpp | 2 +- Code/Source/solver/post.cpp | 53 ++++++++++++++-------------------- Code/Source/solver/post.h | 4 +-- Code/Source/solver/vtk_xml.cpp | 4 +-- 4 files changed, 26 insertions(+), 37 deletions(-) diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index fb7075b43..8c91a54fd 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -157,7 +157,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt if (msh.nFn != 0) { Vector sA(msh.nNo); - post::fib_stretch(simulation, iEq, msh, solutions, sA); + post::fib_stretch(com_mod, iEq, msh, solutions.current.get_displacement(), sA); for (int a = 0; a < msh.nNo; a++) { int Ac = msh.gN(a); I4f(Ac) = sA(a); diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 7d4724931..18e7eb280 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -45,8 +45,6 @@ namespace post { return F; } - void compute_fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, - const Array& lD, Vector& res); } void all_post(Simulation* simulation, Array& res, const SolutionStates& solutions, @@ -723,36 +721,7 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra /// @brief Compute fiber stretch based on 4th invariant: λ = sqrt(I_{4,f}) // -void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) -{ - compute_fib_stretch(simulation->com_mod, iEq, lM, solutions.current.get_displacement(), res); -} - -/// @brief Compute fiber stretch rate dλ/dt via backward finite difference. -// -void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) -{ - auto& com_mod = simulation->com_mod; - const double dt = com_mod.dt; - int nNo = lM.nNo; - - if (res.size() != nNo) { - svmp::raise( - SVMP_HERE, - "[fib_stretch_rate] Expected res size " + std::to_string(nNo) + ", but got " + std::to_string(res.size()) + "."); - } - - Vector lambda_old(nNo); - - compute_fib_stretch(simulation->com_mod, iEq, lM, solutions.current.get_displacement(), res); - compute_fib_stretch(simulation->com_mod, iEq, lM, solutions.old.get_displacement(), lambda_old); - - res = (res - lambda_old) / dt; -} - -namespace { -/// @brief Compute fiber stretch lambda = sqrt(I_{4,f}) from a given displacement field. -void compute_fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, +void fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, const Array& lD, Vector& res) { using namespace consts; @@ -821,6 +790,26 @@ void compute_fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM } } } + +/// @brief Compute fiber stretch rate dλ/dt via backward finite difference. +// +void fib_stretch_rate(const ComMod& com_mod, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) +{ + const double dt = com_mod.dt; + int nNo = lM.nNo; + + if (res.size() != nNo) { + svmp::raise( + SVMP_HERE, + "[fib_stretch_rate] Expected res size " + std::to_string(nNo) + ", but got " + std::to_string(res.size()) + "."); + } + + Vector lambda_old(nNo); + + fib_stretch(com_mod, iEq, lM, solutions.current.get_displacement(), res); + fib_stretch(com_mod, iEq, lM, solutions.old.get_displacement(), lambda_old); + + res = (res - lambda_old) / dt; } void post(Simulation* simulation, const mshType& lM, Array& res, const SolutionStates& solutions, diff --git a/Code/Source/solver/post.h b/Code/Source/solver/post.h index 3eb6c73d4..88982d867 100644 --- a/Code/Source/solver/post.h +++ b/Code/Source/solver/post.h @@ -22,9 +22,9 @@ void fib_algn_post(Simulation* simulation, const mshType& lM, Array& res void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Array& res, const SolutionStates& solutions, const int iEq); -void fib_stretch(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); +void fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, const Array& lD, Vector& res); -void fib_stretch_rate(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); +void fib_stretch_rate(const ComMod& com_mod, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res); void post(Simulation* simulation, const mshType& lM, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp, const int iEq); diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 7d8bb499f..8450f7aaa 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1304,7 +1304,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b case OutputNameType::outGrp_fibStretch: { Vector res(msh.nNo); if (msh.nFn != 0) { - post::fib_stretch(simulation, iEq, msh, solutions, res); + post::fib_stretch(simulation->com_mod, iEq, msh, solutions.current.get_displacement(), res); } for (int a = 0; a < msh.nNo; a++) { d[iM].x(is,a) = res(a); @@ -1314,7 +1314,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b case OutputNameType::outGrp_fibStretchRate: { Vector res(msh.nNo); if (msh.nFn != 0) { - post::fib_stretch_rate(simulation, iEq, msh, solutions, res); + post::fib_stretch_rate(simulation->com_mod, iEq, msh, solutions, res); } for (int a = 0; a < msh.nNo; a++) { d[iM].x(is,a) = res(a); From 4726d54afe458c70e4a2996db0a090513fea1334 Mon Sep 17 00:00:00 2001 From: KB Date: Fri, 22 May 2026 09:54:01 -0700 Subject: [PATCH 14/15] addresses copilot comments on dt failure and inconsistent gnn call --- Code/Source/solver/post.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 18e7eb280..96f4aa520 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -579,8 +579,8 @@ void fib_algn_post(Simulation* simulation, const mshType& lM, Array& res double Jac = 0.0; if (g == 0 || !lM.lShpF) { - auto Nx = lM.Nx.slice(g); - nn::gnn(eNoN, nsd, nsd, Nx, xl, Nx, Jac, F); + auto Nxi = lM.Nx.slice(g); + nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, F); } double w = lM.w(g)*Jac; @@ -798,6 +798,12 @@ void fib_stretch_rate(const ComMod& com_mod, const int iEq, const mshType& lM, c const double dt = com_mod.dt; int nNo = lM.nNo; + if (dt <= 0.0) { + svmp::raise( + SVMP_HERE, + "[fib_stretch_rate] Expected com_mod.dt > 0, but got " + std::to_string(dt) + "."); + } + if (res.size() != nNo) { svmp::raise( SVMP_HERE, From b0cdc87e19e914578b0fc856e4e59044f23870f8 Mon Sep 17 00:00:00 2001 From: KB Date: Fri, 22 May 2026 10:25:57 -0700 Subject: [PATCH 15/15] Addresses element domain check and dummy variable for gnn call comments --- Code/Source/solver/post.cpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 96f4aa520..fbf8278e1 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -740,7 +740,12 @@ void fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, Array Nx(nsd,eNoN); for (int e = 0; e < lM.nEl; e++) { - all_fun::domain(com_mod, lM, iEq, e); + int cDmn = all_fun::domain(com_mod, lM, iEq, e); + auto cPhys = eq.dmn[cDmn].phys; + if (cPhys != EquationType::phys_struct && cPhys != EquationType::phys_ustruct) { + continue; + } + if (lM.eType == ElementType::NRB) { //CALL NRBNNX(lM, e) } @@ -754,10 +759,10 @@ void fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, for (int g = 0; g < lM.nG; g++) { double Jac = 0.0; Array F(nsd,nsd); - Array ksix(nsd,nsd); if (g == 0 || !lM.lShpF) { auto Nxi = lM.Nx.slice(g); - nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, ksix); + Array dummy_ksix(nsd,nsd); + nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, dummy_ksix); } // Compute Deformation Gradient: F = I + grad(u)