From c2ca96212a352daae6cafcf7fa3ae29ead62dc0a Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 22 Jul 2026 09:14:47 +0200 Subject: [PATCH 1/9] Coupling to 0D solver supported for equations other than the first --- Code/Source/solver/ComMod.h | 4 ++ .../solver/CoupledBoundaryCondition.cpp | 31 ++++++++----- Code/Source/solver/Integrator.cpp | 2 +- Code/Source/solver/baf_ini.cpp | 2 +- Code/Source/solver/initialize.cpp | 25 ++++++++-- Code/Source/solver/set_bc.cpp | 46 +++++++++++++++---- 6 files changed, 81 insertions(+), 29 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 2e9f60068..53e7b7eb3 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -790,6 +790,10 @@ class cplBCType { public: cplBCType(); + + /// @brief Index of the equation that this condition is associated with. + unsigned int equationIndex = 0; + /// @brief Is multi-domain active bool coupled = false; diff --git a/Code/Source/solver/CoupledBoundaryCondition.cpp b/Code/Source/solver/CoupledBoundaryCondition.cpp index 2d2ed353d..012e1c647 100644 --- a/Code/Source/solver/CoupledBoundaryCondition.cpp +++ b/Code/Source/solver/CoupledBoundaryCondition.cpp @@ -318,12 +318,16 @@ void CoupledBoundaryCondition::compute_flowrates(ComMod& com_mod, const CmMod& c int nsd = com_mod.nsd; const auto& Yo = solutions.old.get_velocity(); const auto& Yn = solutions.current.get_velocity(); - - Qo_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, 0, solutions, - std::optional(nsd - 1), false, flowrate_cfg_o_); - Qn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, 0, solutions, - std::optional(nsd - 1), false, flowrate_cfg_n_); - + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + + Qo_ = + all_fun::integ(com_mod, cm_mod, *face_, Yo, equation_offset, solutions, + equation_offset + nsd - 1, false, flowrate_cfg_o_); + Qn_ = + all_fun::integ(com_mod, cm_mod, *face_, Yn, equation_offset, solutions, + equation_offset + nsd - 1, false, flowrate_cfg_n_); + if (has_cap_) { const auto [Qo_cap, Qn_cap] = calculate_cap_contribution(com_mod, cm_mod, solutions, flowrate_cfg_o_, flowrate_cfg_n_); @@ -347,12 +351,15 @@ void CoupledBoundaryCondition::compute_pressures(ComMod& com_mod, const CmMod& c double area = face_->area; const auto& Yo = solutions.old.get_velocity(); const auto& Yn = solutions.current.get_velocity(); - - Po_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, nsd, solutions, - std::nullopt, false, flowrate_cfg_o_) / area; - Pn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, nsd, solutions, - std::nullopt, false, flowrate_cfg_n_) / area; - + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + + Po_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, equation_offset + nsd, + solutions, std::nullopt, false, flowrate_cfg_o_) / + area; + Pn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, equation_offset + nsd, + solutions, std::nullopt, false, flowrate_cfg_n_) / + area; } double CoupledBoundaryCondition::get_Qo() const diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 75b9b5447..5fd4557c2 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -92,7 +92,7 @@ bool Integrator::step() { iEqOld = cEq; auto& eq = com_mod.eq[cEq]; - if (com_mod.cplBC.coupled && cEq == 0) { + if (cEq == com_mod.cplBC.equationIndex && com_mod.cplBC.coupled) { #ifdef debug_integrator_step dmsg << "Set coupled BCs " << std::endl; #endif diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index 9355348f2..e3a341636 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -98,7 +98,7 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) // cplBC faces are initialized here // - int iEq = 0; + int iEq = com_mod.cplBC.equationIndex; com_mod.cplBC.fa.resize(com_mod.cplBC.nFa); com_mod.cplBC.xn.resize(com_mod.cplBC.nX); diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 7baf4f330..8decef715 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -3,6 +3,8 @@ // The code here replicates the Fortran code in DISTRIBUTE.f. +#include "Core/Exception.h" + #include "initialize.h" #include "distribute.h" @@ -417,11 +419,24 @@ void initialize(Simulation* simulation, Vector& timeP) nFacesLS = nFacesLS + 1; } - for (auto& bc : com_mod.eq[0].bc) { - // Check for coupled faces (Dir, Neu via cplBC) or Coupled BCs - if (bc.cplBCptr != -1 || utils::btest(bc.bType, static_cast(consts::BoundaryConditionType::bType_Coupled))) { - com_mod.cplBC.coupled = true; - break; + // Check for coupled faces (Dir, Neu via cplBC) or Coupled BCs + for (unsigned int i = 0; i < com_mod.eq.size(); ++i) { + for (auto &bc : com_mod.eq[i].bc) { + + if (bc.cplBCptr != -1 || + utils::btest( + bc.bType, + static_cast(consts::BoundaryConditionType::bType_Coupled))) { + svmp::throw_if( + com_mod.cplBC.coupled && com_mod.cplBC.equationIndex != i, + "Coupled boundary conditions can only be assigned in one equation, " + "but they were assigned in equations " + + std::to_string(com_mod.cplBC.equationIndex) + " and " + + std::to_string(i) + "."); + + com_mod.cplBC.equationIndex = i; + com_mod.cplBC.coupled = true; + } } } diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 65f744727..86c5d1eb5 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -44,7 +44,7 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& dmsg.banner(); #endif - const int iEq = 0; + const int iEq = com_mod.cplBC.equationIndex; // NOTE: For coupling with svZeroDPlus, absTol needs to be > 1e-8 to be compatible with the default convergence tolerance of svZeroDPlus (1e-8) // If this is not true, the finite difference computation of bc.r below results in zero because the perturbation is below the svZeroDPlus tolerance const double absTol = 1.0e-7; @@ -144,8 +144,16 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& else { throw std::runtime_error("[calc_der_cpl_bc] Invalid physics type for 0D coupling"); } - cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, fa, Yo, 0, solutions, nsd-1, false, cfg_o); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, fa, Yn, 0, solutions, nsd-1, false, cfg_n); + + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + cplBC.fa[ptr].Qo = + all_fun::integ(com_mod, cm_mod, fa, Yo, equation_offset, solutions, + equation_offset + nsd - 1, false, cfg_o); + cplBC.fa[ptr].Qn = + all_fun::integ(com_mod, cm_mod, fa, Yn, equation_offset, solutions, + equation_offset + nsd - 1, false, cfg_n); + cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; #ifdef debug_calc_der_cpl_bc @@ -616,7 +624,7 @@ void rcr_init(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& soluti using namespace consts; - const int iEq = 0; + const int iEq = com_mod.cplBC.equationIndex; int nsd = com_mod.nsd; auto& eq = com_mod.eq[iEq]; auto& cplBC = com_mod.cplBC; @@ -751,7 +759,7 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) const int tnNo = com_mod.tnNo; auto& cplBC = com_mod.cplBC; // Yo, Ao, Do now passed as parameters - const int iEq = 0; + const int iEq = com_mod.cplBC.equationIndex; auto& eq = com_mod.eq[iEq]; // Determine current physics @@ -825,16 +833,34 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) else { throw std::runtime_error("[set_bc_cpl] Invalid physics type for 0D coupling"); } - - cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, 0, solutions, nsd-1, false, cfg_o); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, 0, solutions, nsd-1, false, cfg_n); + + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + cplBC.fa[ptr].Qo = all_fun::integ( + com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, equation_offset, + solutions, equation_offset + nsd - 1, false, cfg_o); + cplBC.fa[ptr].Qn = all_fun::integ( + com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, equation_offset, + solutions, equation_offset + nsd - 1, false, cfg_n); + cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; } // Compute avg pressures at 3D Dirichlet boundaries at timesteps n and n+1 else if (utils::btest(bc.bType,iBC_Dir)) { - cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; - cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + cplBC.fa[ptr].Po = + all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, + equation_offset + nsd, solutions, std::nullopt, + false, MechanicalConfigurationType::reference) / + area; + cplBC.fa[ptr].Pn = + all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, + equation_offset + nsd, solutions, std::nullopt, + false, MechanicalConfigurationType::reference) / + area; + cplBC.fa[ptr].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; } From 2ddf52ab15819d05874b2c0586f1379e59ecf670 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 28 Jul 2026 12:11:29 +0200 Subject: [PATCH 2/9] Rename active stress and ionic model files with upper camel case --- .../{active_stress.cpp => ActiveStress.cpp} | 2 +- .../{active_stress.h => ActiveStress.h} | 0 ...filov.cpp => ActiveStressNashPanfilov.cpp} | 2 +- ..._panfilov.h => ActiveStressNashPanfilov.h} | 2 +- ...ive_stress_ode.cpp => ActiveStressODE.cpp} | 2 +- ...{active_stress_ode.h => ActiveStressODE.h} | 2 +- ...eady.cpp => ActiveStressUniformSteady.cpp} | 2 +- ...m_steady.h => ActiveStressUniformSteady.h} | 2 +- ...dy.cpp => ActiveStressUniformUnsteady.cpp} | 2 +- ...steady.h => ActiveStressUniformUnsteady.h} | 4 ++-- Code/Source/solver/CMakeLists.txt | 22 +++++++++---------- Code/Source/solver/CepMod.h | 2 +- Code/Source/solver/ComMod.h | 2 +- .../{ionic_model.cpp => IonicModel.cpp} | 2 +- .../solver/{ionic_model.h => IonicModel.h} | 0 ...nfilov.cpp => IonicModelAlievPanfilov.cpp} | 2 +- ...v_panfilov.h => IonicModelAlievPanfilov.h} | 2 +- ...o_orovio.cpp => IonicModelBuenoOrovio.cpp} | 2 +- ...bueno_orovio.h => IonicModelBuenoOrovio.h} | 2 +- ...agumo.cpp => IonicModelFitzHughNagumo.cpp} | 2 +- ...gh_nagumo.h => IonicModelFitzHughNagumo.h} | 2 +- .../{ionic_ttp.cpp => IonicModelTTP.cpp} | 3 +-- .../solver/{ionic_ttp.h => IonicModelTTP.h} | 2 +- Code/Source/solver/distribute.cpp | 2 +- Code/Source/solver/read_files.cpp | 4 ++-- 25 files changed, 35 insertions(+), 36 deletions(-) rename Code/Source/solver/{active_stress.cpp => ActiveStress.cpp} (98%) rename Code/Source/solver/{active_stress.h => ActiveStress.h} (100%) rename Code/Source/solver/{active_stress_nash_panfilov.cpp => ActiveStressNashPanfilov.cpp} (97%) rename Code/Source/solver/{active_stress_nash_panfilov.h => ActiveStressNashPanfilov.h} (99%) rename Code/Source/solver/{active_stress_ode.cpp => ActiveStressODE.cpp} (97%) rename Code/Source/solver/{active_stress_ode.h => ActiveStressODE.h} (99%) rename Code/Source/solver/{active_stress_uniform_steady.cpp => ActiveStressUniformSteady.cpp} (92%) rename Code/Source/solver/{active_stress_uniform_steady.h => ActiveStressUniformSteady.h} (98%) rename Code/Source/solver/{active_stress_uniform_unsteady.cpp => ActiveStressUniformUnsteady.cpp} (96%) rename Code/Source/solver/{active_stress_uniform_unsteady.h => ActiveStressUniformUnsteady.h} (99%) rename Code/Source/solver/{ionic_model.cpp => IonicModel.cpp} (99%) rename Code/Source/solver/{ionic_model.h => IonicModel.h} (100%) rename Code/Source/solver/{ionic_aliev_panfilov.cpp => IonicModelAlievPanfilov.cpp} (98%) rename Code/Source/solver/{ionic_aliev_panfilov.h => IonicModelAlievPanfilov.h} (99%) rename Code/Source/solver/{ionic_bueno_orovio.cpp => IonicModelBuenoOrovio.cpp} (99%) rename Code/Source/solver/{ionic_bueno_orovio.h => IonicModelBuenoOrovio.h} (99%) rename Code/Source/solver/{ionic_fitzhugh_nagumo.cpp => IonicModelFitzHughNagumo.cpp} (97%) rename Code/Source/solver/{ionic_fitzhugh_nagumo.h => IonicModelFitzHughNagumo.h} (99%) rename Code/Source/solver/{ionic_ttp.cpp => IonicModelTTP.cpp} (99%) rename Code/Source/solver/{ionic_ttp.h => IonicModelTTP.h} (99%) diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/ActiveStress.cpp similarity index 98% rename from Code/Source/solver/active_stress.cpp rename to Code/Source/solver/ActiveStress.cpp index d8630565d..9bf16a043 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress.h" +#include "ActiveStress.h" bool supports_active_stress(const consts::EquationType eq_type) { return eq_type == consts::EquationType::phys_struct || diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/ActiveStress.h similarity index 100% rename from Code/Source/solver/active_stress.h rename to Code/Source/solver/ActiveStress.h diff --git a/Code/Source/solver/active_stress_nash_panfilov.cpp b/Code/Source/solver/ActiveStressNashPanfilov.cpp similarity index 97% rename from Code/Source/solver/active_stress_nash_panfilov.cpp rename to Code/Source/solver/ActiveStressNashPanfilov.cpp index 3e23f21b2..f0e72b617 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.cpp +++ b/Code/Source/solver/ActiveStressNashPanfilov.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_nash_panfilov.h" +#include "ActiveStressNashPanfilov.h" void NashPanfilov::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/ActiveStressNashPanfilov.h similarity index 99% rename from Code/Source/solver/active_stress_nash_panfilov.h rename to Code/Source/solver/ActiveStressNashPanfilov.h index a164be7c6..6dc1f0907 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.h +++ b/Code/Source/solver/ActiveStressNashPanfilov.h @@ -4,7 +4,7 @@ #ifndef ACTIVE_STRESS_NASH_PANFILOV_H #define ACTIVE_STRESS_NASH_PANFILOV_H -#include "active_stress_ode.h" +#include "ActiveStressODE.h" /** * @brief Nash-Panfilov active stress model. diff --git a/Code/Source/solver/active_stress_ode.cpp b/Code/Source/solver/ActiveStressODE.cpp similarity index 97% rename from Code/Source/solver/active_stress_ode.cpp rename to Code/Source/solver/ActiveStressODE.cpp index 99909cf7f..9b9ca9e17 100644 --- a/Code/Source/solver/active_stress_ode.cpp +++ b/Code/Source/solver/ActiveStressODE.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_ode.h" +#include "ActiveStressODE.h" void ActiveStressODE::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { diff --git a/Code/Source/solver/active_stress_ode.h b/Code/Source/solver/ActiveStressODE.h similarity index 99% rename from Code/Source/solver/active_stress_ode.h rename to Code/Source/solver/ActiveStressODE.h index f0696eabe..dd038dcd0 100644 --- a/Code/Source/solver/active_stress_ode.h +++ b/Code/Source/solver/ActiveStressODE.h @@ -4,7 +4,7 @@ #ifndef ACTIVE_STRESS_ODE_H #define ACTIVE_STRESS_ODE_H -#include "active_stress.h" +#include "ActiveStress.h" /** * @brief Abstract ODE-based active stress model. diff --git a/Code/Source/solver/active_stress_uniform_steady.cpp b/Code/Source/solver/ActiveStressUniformSteady.cpp similarity index 92% rename from Code/Source/solver/active_stress_uniform_steady.cpp rename to Code/Source/solver/ActiveStressUniformSteady.cpp index 3c72855de..46ece414f 100644 --- a/Code/Source/solver/active_stress_uniform_steady.cpp +++ b/Code/Source/solver/ActiveStressUniformSteady.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_uniform_steady.h" +#include "ActiveStressUniformSteady.h" void UniformSteadyActiveStress::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { diff --git a/Code/Source/solver/active_stress_uniform_steady.h b/Code/Source/solver/ActiveStressUniformSteady.h similarity index 98% rename from Code/Source/solver/active_stress_uniform_steady.h rename to Code/Source/solver/ActiveStressUniformSteady.h index 93676f9b0..c50e81fd7 100644 --- a/Code/Source/solver/active_stress_uniform_steady.h +++ b/Code/Source/solver/ActiveStressUniformSteady.h @@ -4,7 +4,7 @@ #ifndef ACTIVE_STRESS_UNIFORM_STEADY_H #define ACTIVE_STRESS_UNIFORM_STEADY_H -#include "active_stress.h" +#include "ActiveStress.h" /** * @brief Uniform and steady active stress model. diff --git a/Code/Source/solver/active_stress_uniform_unsteady.cpp b/Code/Source/solver/ActiveStressUniformUnsteady.cpp similarity index 96% rename from Code/Source/solver/active_stress_uniform_unsteady.cpp rename to Code/Source/solver/ActiveStressUniformUnsteady.cpp index 5d02d76b2..48a070d80 100644 --- a/Code/Source/solver/active_stress_uniform_unsteady.cpp +++ b/Code/Source/solver/ActiveStressUniformUnsteady.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_uniform_unsteady.h" +#include "ActiveStressUniformUnsteady.h" #include #include diff --git a/Code/Source/solver/active_stress_uniform_unsteady.h b/Code/Source/solver/ActiveStressUniformUnsteady.h similarity index 99% rename from Code/Source/solver/active_stress_uniform_unsteady.h rename to Code/Source/solver/ActiveStressUniformUnsteady.h index 696be51ab..98438fde0 100644 --- a/Code/Source/solver/active_stress_uniform_unsteady.h +++ b/Code/Source/solver/ActiveStressUniformUnsteady.h @@ -4,9 +4,9 @@ #ifndef ACTIVE_STRESS_UNIFORM_UNSTEADY_H #define ACTIVE_STRESS_UNIFORM_UNSTEADY_H -#include "FourierInterpolation.h" +#include "ActiveStress.h" -#include "active_stress.h" +#include "FourierInterpolation.h" /** * @brief Uniform and time dependent active stress model. diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c5ab81146..f78cc4114 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -223,17 +223,17 @@ set(CSRCS Timer.h ArtificialNeuralNetMaterial.h ArtificialNeuralNetMaterial.cpp - ionic_model.cpp - ionic_aliev_panfilov.cpp - ionic_bueno_orovio.cpp - ionic_fitzhugh_nagumo.cpp - ionic_ttp.cpp - - active_stress.cpp - active_stress_uniform_steady.cpp - active_stress_uniform_unsteady.cpp - active_stress_ode.cpp - active_stress_nash_panfilov.cpp + IonicModel.cpp + IonicModelAlievPanfilov.cpp + IonicModelBuenoOrovio.cpp + IonicModelFitzHughNagumo.cpp + IonicModelTTP.cpp + + ActiveStress.cpp + ActiveStressUniformSteady.cpp + ActiveStressUniformUnsteady.cpp + ActiveStressODE.cpp + ActiveStressNashPanfilov.cpp SPLIT.c diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index 6d280228b..950b94896 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -12,7 +12,7 @@ #define CEP_MOD_H #include "consts.h" -#include "ionic_model.h" +#include "IonicModel.h" #include "Array.h" #include "Vector.h" diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 53e7b7eb3..1d378c2c6 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -22,7 +22,7 @@ #include "SolutionStates.h" #include "Timer.h" #include "Vector.h" -#include "active_stress.h" +#include "ActiveStress.h" #include "DebugMsg.h" diff --git a/Code/Source/solver/ionic_model.cpp b/Code/Source/solver/IonicModel.cpp similarity index 99% rename from Code/Source/solver/ionic_model.cpp rename to Code/Source/solver/IonicModel.cpp index 6bcfdd0bd..1472d1772 100644 --- a/Code/Source/solver/ionic_model.cpp +++ b/Code/Source/solver/IonicModel.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_model.h" +#include "IonicModel.h" #include "ComMod.h" #include "Parameters.h" diff --git a/Code/Source/solver/ionic_model.h b/Code/Source/solver/IonicModel.h similarity index 100% rename from Code/Source/solver/ionic_model.h rename to Code/Source/solver/IonicModel.h diff --git a/Code/Source/solver/ionic_aliev_panfilov.cpp b/Code/Source/solver/IonicModelAlievPanfilov.cpp similarity index 98% rename from Code/Source/solver/ionic_aliev_panfilov.cpp rename to Code/Source/solver/IonicModelAlievPanfilov.cpp index 2d0476f41..9914ddbce 100644 --- a/Code/Source/solver/ionic_aliev_panfilov.cpp +++ b/Code/Source/solver/IonicModelAlievPanfilov.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_aliev_panfilov.h" +#include "IonicModelAlievPanfilov.h" void AlievPanfilov::read_parameters(const IonicModelParameters ¶ms) { IonicModel::read_parameters(params); diff --git a/Code/Source/solver/ionic_aliev_panfilov.h b/Code/Source/solver/IonicModelAlievPanfilov.h similarity index 99% rename from Code/Source/solver/ionic_aliev_panfilov.h rename to Code/Source/solver/IonicModelAlievPanfilov.h index 9f0c956ee..5e1edef8a 100644 --- a/Code/Source/solver/ionic_aliev_panfilov.h +++ b/Code/Source/solver/IonicModelAlievPanfilov.h @@ -4,7 +4,7 @@ #ifndef IONIC_ALIEV_PANFILOV_H #define IONIC_ALIEV_PANFILOV_H -#include "ionic_model.h" +#include "IonicModel.h" #include "Vector.h" diff --git a/Code/Source/solver/ionic_bueno_orovio.cpp b/Code/Source/solver/IonicModelBuenoOrovio.cpp similarity index 99% rename from Code/Source/solver/ionic_bueno_orovio.cpp rename to Code/Source/solver/IonicModelBuenoOrovio.cpp index fd8001dcb..88617071e 100644 --- a/Code/Source/solver/ionic_bueno_orovio.cpp +++ b/Code/Source/solver/IonicModelBuenoOrovio.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_bueno_orovio.h" +#include "IonicModelBuenoOrovio.h" void BuenoOrovio::read_parameters(const IonicModelParameters ¶ms) { IonicModel::read_parameters(params); diff --git a/Code/Source/solver/ionic_bueno_orovio.h b/Code/Source/solver/IonicModelBuenoOrovio.h similarity index 99% rename from Code/Source/solver/ionic_bueno_orovio.h rename to Code/Source/solver/IonicModelBuenoOrovio.h index fc8868f4f..e59134945 100644 --- a/Code/Source/solver/ionic_bueno_orovio.h +++ b/Code/Source/solver/IonicModelBuenoOrovio.h @@ -4,7 +4,7 @@ #ifndef IONIC_BUENO_OROVIO_H #define IONIC_BUENO_OROVIO_H -#include "ionic_model.h" +#include "IonicModel.h" #include "Vector.h" #include "utils.h" diff --git a/Code/Source/solver/ionic_fitzhugh_nagumo.cpp b/Code/Source/solver/IonicModelFitzHughNagumo.cpp similarity index 97% rename from Code/Source/solver/ionic_fitzhugh_nagumo.cpp rename to Code/Source/solver/IonicModelFitzHughNagumo.cpp index 608196377..06c6585e2 100644 --- a/Code/Source/solver/ionic_fitzhugh_nagumo.cpp +++ b/Code/Source/solver/IonicModelFitzHughNagumo.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_fitzhugh_nagumo.h" +#include "IonicModelFitzHughNagumo.h" void FitzHughNagumo::read_parameters(const IonicModelParameters ¶ms) { IonicModel::read_parameters(params); diff --git a/Code/Source/solver/ionic_fitzhugh_nagumo.h b/Code/Source/solver/IonicModelFitzHughNagumo.h similarity index 99% rename from Code/Source/solver/ionic_fitzhugh_nagumo.h rename to Code/Source/solver/IonicModelFitzHughNagumo.h index 572aa2ad4..2ea6b4c77 100644 --- a/Code/Source/solver/ionic_fitzhugh_nagumo.h +++ b/Code/Source/solver/IonicModelFitzHughNagumo.h @@ -4,7 +4,7 @@ #ifndef IONIC_FITZHUGH_NAGUMO_H #define IONIC_FITZHUGH_NAGUMO_H -#include "ionic_model.h" +#include "IonicModel.h" #include "Vector.h" #include "utils.h" diff --git a/Code/Source/solver/ionic_ttp.cpp b/Code/Source/solver/IonicModelTTP.cpp similarity index 99% rename from Code/Source/solver/ionic_ttp.cpp rename to Code/Source/solver/IonicModelTTP.cpp index c7cf58152..15e6e9979 100644 --- a/Code/Source/solver/ionic_ttp.cpp +++ b/Code/Source/solver/IonicModelTTP.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_ttp.h" +#include "IonicModelTTP.h" void TTP::read_parameters(const IonicModelParameters ¶ms) { IonicModel::read_parameters(params); @@ -456,5 +456,4 @@ Vector TTP::getf(const unsigned int zone_id, const Vector &X, return dX; } - REGISTER_IONIC_MODEL("TTP", TTP); \ No newline at end of file diff --git a/Code/Source/solver/ionic_ttp.h b/Code/Source/solver/IonicModelTTP.h similarity index 99% rename from Code/Source/solver/ionic_ttp.h rename to Code/Source/solver/IonicModelTTP.h index 1a90408b2..2602c036c 100644 --- a/Code/Source/solver/ionic_ttp.h +++ b/Code/Source/solver/IonicModelTTP.h @@ -4,7 +4,7 @@ #ifndef IONIC_TTP_H #define IONIC_TTP_H -#include "ionic_model.h" +#include "IonicModel.h" #include "Parameters.h" diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 00986613e..244dd1667 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -18,7 +18,7 @@ #include #include -#include "ionic_model.h" +#include "IonicModel.h" extern "C" { diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index befd3b9e4..7c4f770eb 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -7,10 +7,10 @@ #include "Core/Exception.h" #include "FE/Common/FEException.h" -#include "active_stress.h" +#include "ActiveStress.h" #include "all_fun.h" #include "consts.h" -#include "ionic_model.h" +#include "IonicModel.h" #include "read_msh.h" #include "vtk_xml.h" From 9cb185b0eada2e16e58e574636373c548484d575 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 4 Aug 2026 15:50:19 -0500 Subject: [PATCH 3/9] nn::gnnb and all_fun::integ take displacement index as input to compute deformed configurations without assuming equation order --- .../solver/CoupledBoundaryCondition.cpp | 18 +-- Code/Source/solver/all_fun.cpp | 52 ++++--- Code/Source/solver/all_fun.h | 127 ++++++++++++++++- Code/Source/solver/baf_ini.cpp | 24 +++- Code/Source/solver/cmm.cpp | 4 +- Code/Source/solver/eq_assem.cpp | 13 +- Code/Source/solver/nn.cpp | 26 ++-- Code/Source/solver/nn.h | 47 ++++++- Code/Source/solver/ris.cpp | 27 +++- Code/Source/solver/set_bc.cpp | 131 +++++++++++------- Code/Source/solver/txt.cpp | 21 ++- 11 files changed, 365 insertions(+), 125 deletions(-) diff --git a/Code/Source/solver/CoupledBoundaryCondition.cpp b/Code/Source/solver/CoupledBoundaryCondition.cpp index 012e1c647..c147799b6 100644 --- a/Code/Source/solver/CoupledBoundaryCondition.cpp +++ b/Code/Source/solver/CoupledBoundaryCondition.cpp @@ -321,12 +321,12 @@ void CoupledBoundaryCondition::compute_flowrates(ComMod& com_mod, const CmMod& c const unsigned int equation_offset = com_mod.eq[com_mod.cplBC.equationIndex].s; - Qo_ = - all_fun::integ(com_mod, cm_mod, *face_, Yo, equation_offset, solutions, - equation_offset + nsd - 1, false, flowrate_cfg_o_); - Qn_ = - all_fun::integ(com_mod, cm_mod, *face_, Yn, equation_offset, solutions, - equation_offset + nsd - 1, false, flowrate_cfg_n_); + Qo_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, equation_offset, + solutions, equation_offset + nsd - 1, false, + flowrate_cfg_o_, equation_offset); + Qn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, equation_offset, + solutions, equation_offset + nsd - 1, false, + flowrate_cfg_n_, equation_offset); if (has_cap_) { const auto [Qo_cap, Qn_cap] = @@ -355,10 +355,12 @@ void CoupledBoundaryCondition::compute_pressures(ComMod& com_mod, const CmMod& c com_mod.eq[com_mod.cplBC.equationIndex].s; Po_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, equation_offset + nsd, - solutions, std::nullopt, false, flowrate_cfg_o_) / + solutions, std::nullopt, false, flowrate_cfg_o_, + equation_offset) / area; Pn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, equation_offset + nsd, - solutions, std::nullopt, false, flowrate_cfg_n_) / + solutions, std::nullopt, false, flowrate_cfg_n_, + equation_offset) / area; } diff --git a/Code/Source/solver/all_fun.cpp b/Code/Source/solver/all_fun.cpp index 623b77412..fedc7b3c0 100644 --- a/Code/Source/solver/all_fun.cpp +++ b/Code/Source/solver/all_fun.cpp @@ -676,18 +676,10 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, - const SolutionStates& solutions, bool pFlag, MechanicalConfigurationType cfg) -{ +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Vector &s, const SolutionStates &solutions, + bool pFlag, MechanicalConfigurationType cfg, + const unsigned int displacement_index) { using namespace consts; #define n_debug_integ_s #ifdef debug_integ_s @@ -811,7 +803,8 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co if (!isIB) { // Get normal vector in cfg configuration auto Nx = fs.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, insd, fs.eNoN, Nx, n, solutions, cfg); + nn::gnnb(com_mod, lFa, e, g, nsd, insd, fs.eNoN, Nx, n, solutions, cfg, + displacement_index); } // Calculating the Jacobian (encodes area of face element) @@ -835,7 +828,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co } result = com_mod.cm.reduce(cm_mod, result); - return result; + return result; } /// @brief This routine integrates vector field s dotted with the face normal n @@ -849,9 +842,10 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co /// @param pFlag flag for using Taylor-Hood function space for pressure /// @param cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference. // -double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, - const Array& s, const SolutionStates& solutions, MechanicalConfigurationType cfg) -{ +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const SolutionStates &solutions, + MechanicalConfigurationType cfg, + const unsigned int displacement_index) { using namespace consts; #define n_debug_integ_V @@ -931,7 +925,8 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, if (!isIB) { // Get normal vector in cfg configuration auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, cfg); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, lFa.eNoN, Nx, n, solutions, + cfg, displacement_index); //CALL GNNB(lFa, e, g, nsd-1, lFa.eNoN, lFa.Nx(:,:,g), n) } else { //CALL GNNIB(lFa, e, g, n) @@ -963,7 +958,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, result = cm.reduce(cm_mod, result); - return result; + return result; } /// @brief This routine integrate s(l:u,:) over the surface faId, where s is an @@ -983,9 +978,11 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, /// @param cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference. /// // -double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, - const Array& s, const int l, const SolutionStates& solutions, std::optional uo, bool THflag, MechanicalConfigurationType cfg) -{ +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const int l, + const SolutionStates &solutions, std::optional uo, + bool THflag, MechanicalConfigurationType cfg, + const unsigned int displacement_index) { using namespace consts; #define n_debug_integ_g @@ -1040,22 +1037,23 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, vec(n,a) = s(i,a); } } - result = integ(com_mod, cm_mod, lFa, vec, solutions, cfg); - // If s scalar, integrate as scalar + result = + integ(com_mod, cm_mod, lFa, vec, solutions, cfg, displacement_index); + // If s scalar, integrate as scalar } else if (l == u) { Vector sclr(nNo); for (int a = 0; a < nNo; a++) { sclr(a) = s(l,a); } - result = integ(com_mod, cm_mod, lFa, sclr, solutions, flag, cfg); + result = integ(com_mod, cm_mod, lFa, sclr, solutions, flag, cfg, + displacement_index); } else { throw std::runtime_error("Unexpected dof in integ"); } - return result; + return result; } - bool is_domain(const ComMod& com_mod, const eqType& eq, const int node, const consts::EquationType phys) { bool result = false; diff --git a/Code/Source/solver/all_fun.h b/Code/Source/solver/all_fun.h index 2aa588cfa..7dd3637bb 100644 --- a/Code/Source/solver/all_fun.h +++ b/Code/Source/solver/all_fun.h @@ -35,13 +35,126 @@ namespace all_fun { double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, int l, int u, const SolutionStates& solutions, bool pFlag=false); - double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Vector& s, - const SolutionStates& solutions, bool pFlag=false, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); - - double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s, - const int l, const SolutionStates& solutions, std::optional uo=std::nullopt, bool THflag=false, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); - - double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s, const SolutionStates& solutions, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); + /** + * @brief Integrate a scalar field over a boundary face. + * + * Reproduces 'FUNCTION IntegS(lFa, s, pflag)'. The scalar field s is + * integrated over the face, i.e. this computes + * \f[ + * \int_{\Gamma} s \, d\Gamma, + * \f] + * where \f$\Gamma\f$ is the boundary face. + * + * For simulations involving structural displacement, this function allows + * computing the integral in any of the following configurations: + * - reference configuration (the mesh is not displaced); + * - current configuration (the mesh is displaced by the current displacement + * field); + * - old configuration (the mesh is displaced by the displacement field from + * previous time step). + * + * @param[in] com_mod The common module. + * @param[in] cm_mod The communication module containing MPI data. + * @param[in] lFa The boundary face over which the integral is computed. + * @param[in] s The scalar value at each node of the mesh. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @param[in] pFlag Whether to use the Taylor-Hood function space for the + * pressure field. + * @param[in] cfg The configuration in which the integral is computed + * (reference, old or current). + * @param[in] displacement_index The index of the displacement field in the + * solution arrays. This should correspond to the start index of the + * equation that solves for the displacement. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Vector &s, const SolutionStates &solutions, + bool pFlag, consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); + + /** + * @brief Integrate one or more components of a field over a boundary face. + * + * Reproduces 'FUNCTION IntegG(lFa, s, l, u, THflag)'. Rows l to u of s are + * integrated over the face. When they span the spatial dimensions the field + * is dotted with the outward surface normal (i.e. a flux is computed), + * \f[ + * \int_{\Gamma} \sum_{i=l}^{u} s_i \, n_i \, d\Gamma, + * \f] + * otherwise the single row l is integrated as a scalar, + * \f[ + * \int_{\Gamma} s_l \, d\Gamma, + * \f] + * where \f$\Gamma\f$ is the boundary face and \f$\mathbf{n}\f$ its outward + * unit normal. + * + * For simulations involving structural displacement, this function allows + * computing the integral in any of the following configurations: + * - reference configuration (the mesh is not displaced); + * - current configuration (the mesh is displaced by the current displacement + * field); + * - old configuration (the mesh is displaced by the displacement field from + * previous time step). + * + * @param[in] com_mod The common module. + * @param[in] cm_mod The communication module containing MPI data. + * @param[in] lFa The boundary face over which the integral is computed. + * @param[in] s The field value at each node of the mesh. + * @param[in] l The first row of s to integrate. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @param[in] uo The last row of s to integrate. Defaults to l, i.e. a single + * component. + * @param[in] THflag Whether to use the Taylor-Hood function space for the + * pressure field. + * @param[in] cfg The configuration in which the integral is computed + * (reference, old or current). + * @param[in] displacement_index The index of the displacement field in the + * solution arrays. This should correspond to the start index of the + * equation that solves for the displacement. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const int l, + const SolutionStates &solutions, std::optional uo, + bool THflag, consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); + + /** + * @brief Integrate the flux of a vector field over a boundary face. + * + * Reproduces 'FUNCTION IntegV(lFa, s)'. The vector field s (one component per + * spatial dimension at each node) is dotted with the outward surface normal + * and integrated over the face, i.e. this computes + * \f[ + * \int_{\Gamma} \mathbf{s} \cdot \mathbf{n} \, d\Gamma, + * \f] + * where \f$\Gamma\f$ is the boundary face and \f$\mathbf{n}\f$ its outward + * unit normal. + * + * For simulations involving structural displacement, this function allows + * computing the integral in any of the following configurations: + * - reference configuration (the mesh is not displaced); + * - current configuration (the mesh is displaced by the current displacement + * field); + * - old configuration (the mesh is displaced by the displacement field from + * previous time step). + * + * @param[in] com_mod The common module. + * @param[in] cm_mod The communication module containing MPI data. + * @param[in] lFa The boundary face over which the integral is computed. + * @param[in] s The vector value at each node of the mesh. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @param[in] cfg The configuration in which the integral is computed + * (reference, old or current). + * @param[in] displacement_index The index of the displacement field in the + * solution arrays. This should correspond to the start index of the + * equation that solves for the displacement. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const SolutionStates &solutions, + consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); bool is_domain(const ComMod& com_mod, const eqType& eq, const int node, const consts::EquationType phys); diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index e3a341636..4d9cf136b 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -337,7 +337,11 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l } else if (btest(lBc.bType, iBC_para)) { Vector center(3); for (int i = 0; i < nsd; i++) { - center(i) = all_fun::integ(com_mod, cm_mod, lFa, com_mod.x, i, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference) / lFa.area; + center(i) = all_fun::integ(com_mod, cm_mod, lFa, com_mod.x, i, solutions, + std::nullopt, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / + lFa.area; } // gNodes is one if a node located on the boundary (beside iFa) @@ -445,8 +449,10 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l // Normalizing the profile for flux // double tmp = 1.0; - if (btest(lBc.bType, enum_int(BoundaryConditionType::bType_flx))) { - tmp = all_fun::integ(com_mod, cm_mod, lFa, s, solutions, false, consts::MechanicalConfigurationType::reference); + if (btest(lBc.bType, enum_int(BoundaryConditionType::bType_flx))) { + tmp = all_fun::integ(com_mod, cm_mod, lFa, s, solutions, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); if (is_zero(tmp)) { tmp = 1.0; throw std::runtime_error("Face '" + lFa.name + "' used for a BC has no non-zero node."); @@ -488,7 +494,9 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution // Vector sA(com_mod.tnNo); sA = 1.0; - double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, consts::MechanicalConfigurationType::reference); + double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); #ifdef debug_face_ini dmsg << "Face '" << lFa.name << "' area: " << area; #endif @@ -530,7 +538,9 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, lFa.eNoN, Nx, nV, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); @@ -794,7 +804,9 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp for (int g = 0; g < lFa.nG; g++) { Vector n(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, consts::MechanicalConfigurationType::reference); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, lFa.eNoN, Nx, n, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); diff --git a/Code/Source/solver/cmm.cpp b/Code/Source/solver/cmm.cpp index 3d78b5b96..c1c120ffb 100644 --- a/Code/Source/solver/cmm.cpp +++ b/Code/Source/solver/cmm.cpp @@ -282,7 +282,9 @@ void cmm_b(ComMod& com_mod, const faceType& lFa, const int e, const Array nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, 3, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, 3, Nx, nV, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; diff --git a/Code/Source/solver/eq_assem.cpp b/Code/Source/solver/eq_assem.cpp index 9cea03b30..78e2829ec 100644 --- a/Code/Source/solver/eq_assem.cpp +++ b/Code/Source/solver/eq_assem.cpp @@ -80,7 +80,9 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& for (int g = 0; g < lFa.nG; g++) { Vector nV(nsd); auto Nx = lFa.Nx.rslice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoN, Nx, nV, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -253,7 +255,9 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const // Get surface normal vector Vector nV(nsd); auto Nx_g = lFa.Nx.rslice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoNb, Nx_g, nV, solutions, consts::MechanicalConfigurationType::reference); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoNb, Nx_g, nV, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -330,8 +334,9 @@ void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const S auto cfg = MechanicalConfigurationType::new_timestep; - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, cfg); - // + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, lFa.eNoN, Nx, n, solutions, + cfg, com_mod.eq[com_mod.cEq].s); + for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); for (int i = 0; i < nsd; i++) { diff --git a/Code/Source/solver/nn.cpp b/Code/Source/solver/nn.cpp index 78c6e2207..f73dabd30 100644 --- a/Code/Source/solver/nn.cpp +++ b/Code/Source/solver/nn.cpp @@ -899,18 +899,12 @@ void gnn(const int eNoN, const int nsd, const int insd, Array& Nxi, Arra } } -/// @brief This routine returns a surface normal vector at element "e" and Gauss point -/// 'g' of face 'lFa' that is the normal weighted by Jac, i.e. -/// Jac = norm(n), the Jacobian of the mapping from parent surface element to -/// reference/old/new configuration. -/// -/// cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference -/// -/// Reproduce Fortran 'GNNB'. -// -void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, const int nsd, const int insd, - const int eNoNb, const Array& Nx, Vector& n, const SolutionStates& solutions, MechanicalConfigurationType cfg) -{ +void gnnb(const ComMod &com_mod, const faceType &lFa, const int e, const int g, + const int nsd, const int insd, const int eNoNb, + const Array &Nx, Vector &n, + const SolutionStates &solutions, + consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index) { // Local aliases for displacement arrays const auto& Dn = solutions.current.get_displacement(); const auto& Do = solutions.old.get_displacement(); @@ -998,6 +992,10 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, if (com_mod.mvMsh) { for (int i = 0; i < lX.nrows(); i++) { // Add mesh displacement + // Notice that this assumes that the mesh displacement is stored + // starting at the nsd+1 index of the solution array. This is enforced + // in read_files, by throwing an exception if the equations are not + // ordered correctly. lX(i,a) = lX(i,a) + Do(i+nsd+1,Ac); } } @@ -1009,13 +1007,13 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, case MechanicalConfigurationType::old_timestep: for (int i = 0; i < lX.nrows(); i++) { // Add displacement at timestep n - lX(i,a) = lX(i,a) + Do(i,Ac); + lX(i, a) = lX(i, a) + Do(displacement_index + i, Ac); } break; case MechanicalConfigurationType::new_timestep: for (int i = 0; i < lX.nrows(); i++) { // Add displacement at timestep n+1 - lX(i,a) = lX(i,a) + Dn(i,Ac); + lX(i, a) = lX(i, a) + Dn(displacement_index + i, Ac); } break; default: diff --git a/Code/Source/solver/nn.h b/Code/Source/solver/nn.h index 0bd862dd8..ecf432867 100644 --- a/Code/Source/solver/nn.h +++ b/Code/Source/solver/nn.h @@ -38,8 +38,51 @@ namespace nn { void gnn(const int eNoN, const int nsd, const int insd, Array& Nxi, Array& x, Array& Nx, double& Jac, Array& ks); - void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, const int nsd, const int insd, - const int eNoNb, const Array& Nx, Vector& n, const SolutionStates& solutions, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); + /** + * @brief Return the surface normal vector at a given element and Gauss point. + * + * 'g' of face 'lFa' that is the normal weighted by Jac, i.e. + * Jac = norm(n), the Jacobian of the mapping from parent surface element to + * reference/old/new configuration. + * + * For simulations involving structural displacement, this function allows + * computing the normal vector in any of the following configurations: + * - reference configuration (the mesh is not displaced); + * - current configuration (the mesh is displaced by the current displacement + * field); + * - old configuration (the mesh is displaced by the displacement field from + * previous time step). + * + * @todo[michelebucelli] Many of the arguments passed to this function are + * redundant, as they can be inferred from the others. For example, insd and + * eNoNb can be inferred from the size of Nx. Those arguments should be + * removed to simplify the function signature and make it less error-prone. + * + * @param[in] com_mod The common module. + * @param[in] lFa The boundary face for which the normal vector is computed. + * @param[in] e The face-local index of the element for which the normal + * vector is computed. + * @param[in] g The Gauss point index for which the normal vector is computed. + * @param[in] nsd The number of spatial dimensions. + * @param[in] insd The intrinsic dimension of the boundary surface (typically + * nsd - 1). + * @param[in] eNoNb Number of nodes on a boundary face element. + * @param[in] Nx The shape function derivatives at the Gauss point. + * @param[out] n The computed normal vector. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @param[in] cfg The configuration in which the normal vector is computed + * (reference, old or current). + * @param[in] displacement_index The index of the displacement field in the + * solution arrays. This should correspond to the start index of the + * equation that solves for the displacement. + */ + void gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const int nsd, const int insd, const int eNoNb, + const Array &Nx, Vector &n, + const SolutionStates &solutions, + consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); void gnns(const int nsd, const int eNoN, const Array& Nxi, Array& xl, Vector& nV, Array& gCov, Array& gCnv); diff --git a/Code/Source/solver/ris.cpp b/Code/Source/solver/ris.cpp index cfa5b81d8..6057eaed5 100644 --- a/Code/Source/solver/ris.cpp +++ b/Code/Source/solver/ris.cpp @@ -59,7 +59,12 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) int iM = RIS.lst(i,0,iProj); int iFa = RIS.lst(i,1,iProj); double tmp = msh[iM].fa[iFa].area; - RIS.meanP(iProj,i) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference)/tmp; + RIS.meanP(iProj, i) = + all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, + std::nullopt, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / + tmp; } } @@ -77,7 +82,10 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) } int iM = RIS.lst(0,0,iProj); int iFa = RIS.lst(0,1,iProj); - RIS.meanFl(iProj) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); + RIS.meanFl(iProj) = all_fun::integ( + com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m - 1, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); if (cm.mas(cm_mod)) { std::cout << "For RIS projection: " << iProj << std::endl; @@ -480,8 +488,14 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solution sA = 1.0; lFa = msh[iM].fa[iFa]; // such update may be not correct - tmp_new = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, consts::MechanicalConfigurationType::reference); - meanP = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference)/tmp_new; + tmp_new = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); + meanP = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, + m - 1, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / + tmp_new; // For the velocity m = nsd; @@ -496,7 +510,10 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solution } } - meanFl = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); + meanFl = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, + solutions, m - 1, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); std::cout << "The average pressure is: " << meanP << std::endl; std::cout << "The pressure from 0D is: " << eq[cEq].bc[iBc].g << std::endl; diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 86c5d1eb5..931114db3 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -147,12 +147,12 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& const unsigned int equation_offset = com_mod.eq[com_mod.cplBC.equationIndex].s; - cplBC.fa[ptr].Qo = - all_fun::integ(com_mod, cm_mod, fa, Yo, equation_offset, solutions, - equation_offset + nsd - 1, false, cfg_o); - cplBC.fa[ptr].Qn = - all_fun::integ(com_mod, cm_mod, fa, Yn, equation_offset, solutions, - equation_offset + nsd - 1, false, cfg_n); + cplBC.fa[ptr].Qo = all_fun::integ( + com_mod, cm_mod, fa, Yo, equation_offset, solutions, + equation_offset + nsd - 1, false, cfg_o, equation_offset); + cplBC.fa[ptr].Qn = all_fun::integ( + com_mod, cm_mod, fa, Yn, equation_offset, solutions, + equation_offset + nsd - 1, false, cfg_n, equation_offset); cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; @@ -166,8 +166,19 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& // Compute avg pressures at 3D Dirichlet boundaries at timesteps n and n+1 else if (utils::btest(bc.bType, iBC_Dir)) { double area = fa.area; - cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; - cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, fa, Yn, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; + cplBC.fa[ptr].Po = + all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, + std::nullopt, false, + MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / + area; + cplBC.fa[ptr].Pn = + all_fun::integ(com_mod, cm_mod, fa, Yn, nsd, solutions, + std::nullopt, false, + MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / + area; + cplBC.fa[ptr].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; #ifdef debug_calc_der_cpl_bc @@ -643,8 +654,15 @@ void rcr_init(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& soluti if (cplBC.initRCR) { auto& fa = com_mod.msh[iM].fa[iFa]; double area = fa.area; - double Qo = all_fun::integ(com_mod, cm_mod, fa, Yo, 0, solutions, nsd-1, false, MechanicalConfigurationType::reference); - double Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; + double Qo = + all_fun::integ(com_mod, cm_mod, fa, Yo, 0, solutions, nsd - 1, + false, MechanicalConfigurationType::reference, + /* displacement_index = */ 0); + double Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, + std::nullopt, false, + MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / + area; cplBC.xo[ptr] = Po - (Qo * cplBC.fa[ptr].RCR.Rp); } else { cplBC.xo[ptr] = cplBC.fa[ptr].RCR.Xo; @@ -788,7 +806,10 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) faceType& lFa = com_mod.msh[iM].fa[iFa]; Vector sA(com_mod.tnNo); sA = 1.0; - double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, consts::MechanicalConfigurationType::reference); + double area = + all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); baf_ini_ns::bc_ini(com_mod, cm_mod, eq.bc[iBc], lFa, solutions); int ptr = bc.cplBCptr; @@ -838,10 +859,12 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) com_mod.eq[com_mod.cplBC.equationIndex].s; cplBC.fa[ptr].Qo = all_fun::integ( com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, equation_offset, - solutions, equation_offset + nsd - 1, false, cfg_o); + solutions, equation_offset + nsd - 1, false, cfg_o, + equation_offset); cplBC.fa[ptr].Qn = all_fun::integ( com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, equation_offset, - solutions, equation_offset + nsd - 1, false, cfg_n); + solutions, equation_offset + nsd - 1, false, cfg_n, + equation_offset); cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; @@ -853,12 +876,14 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, equation_offset + nsd, solutions, std::nullopt, - false, MechanicalConfigurationType::reference) / + false, MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / area; cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, equation_offset + nsd, solutions, std::nullopt, - false, MechanicalConfigurationType::reference) / + false, MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / area; cplBC.fa[ptr].Qo = 0.0; @@ -1422,7 +1447,9 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const for (int g = 0; g < lFa.nG; g++) { Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoNb, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoNb, Nx, nV, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -1572,37 +1599,43 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const //h(0) = lBc.g; - double Q_3D = all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, - eq.s+nsd-1, false, - consts::MechanicalConfigurationType::reference); - - - h(0) = lBc.g; - //h(0) = lBc.g - lBc.r * std::abs(Q_3D); - - // Backflow kinetic energy correction: when backflow is detected - // (Q < 0), subtract the face-averaged dynamic pressure to further - // reduce the applied traction and damp the incoming flow. - if (Q_3D < 0.0) { - int iM = lFa.iM; - int cDmn_local = all_fun::domain(com_mod, com_mod.msh[iM], cEq, lFa.gE(0)); - double rho = eq.dmn[cDmn_local].prop.at( - consts::PhysicalProperyType::fluid_density); - double beta = eq.dmn[cDmn_local].prop.at( - consts::PhysicalProperyType::backflow_stab); - double A = lFa.area; - if (A > 0.0) { - double u_n = Q_3D / A; // face-averaged normal velocity (< 0) - h(0) -= 0.5 * beta * rho * u_n * u_n; - } - } - } else if (utils::btest(lBc.bType,iBC_res)) { - h(0) = lBc.r * all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, eq.s+nsd-1, false, consts::MechanicalConfigurationType::reference); - - } else if (utils::btest(lBc.bType,iBC_std)) { + double Q_3D = all_fun::integ( + com_mod, cm_mod, lFa, Yn, eq.s, solutions, eq.s + nsd - 1, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); + + h(0) = lBc.g; + // h(0) = lBc.g - lBc.r * std::abs(Q_3D); + + // Backflow kinetic energy correction: when backflow is detected + // (Q < 0), subtract the face-averaged dynamic pressure to further + // reduce the applied traction and damp the incoming flow. + if (Q_3D < 0.0) { + int iM = lFa.iM; + int cDmn_local = + all_fun::domain(com_mod, com_mod.msh[iM], cEq, lFa.gE(0)); + double rho = eq.dmn[cDmn_local].prop.at( + consts::PhysicalProperyType::fluid_density); + double beta = eq.dmn[cDmn_local].prop.at( + consts::PhysicalProperyType::backflow_stab); + double A = lFa.area; + if (A > 0.0) { + double u_n = Q_3D / A; // face-averaged normal velocity (< 0) + h(0) -= 0.5 * beta * rho * u_n * u_n; + } + } + + } else if (utils::btest(lBc.bType, iBC_res)) { + h(0) = lBc.r * + all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, + eq.s + nsd - 1, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); + + } else if (utils::btest(lBc.bType, iBC_std)) { h(0) = lBc.g; - } else if (utils::btest(lBc.bType,iBC_ustd)) { + } else if (utils::btest(lBc.bType, iBC_ustd)) { h = lBc.gt.value(com_mod.time); } else { @@ -1723,7 +1756,9 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondit for (int g = 0; g < lFa.nG; g++) { Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoN, Nx, nV, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -2005,7 +2040,9 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons for (int g = 0; g < lFa.nG; g++) { Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoN, Nx, nV, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); double w = lFa.w(g)*Jac; N = lFa.N.col(g); diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 89f03c907..2347da0e0 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -465,17 +465,30 @@ void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eq if (m == 1) { if (div) { tmp = fa.area; - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference) / tmp; + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, + std::nullopt, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0) / + tmp; } else { if (pFlag && lTH) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, true, consts::MechanicalConfigurationType::reference); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, + std::nullopt, true, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); } else { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, + std::nullopt, false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); } } } else if (m == nsd) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, m - 1, + false, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); } else { throw std::runtime_error("WTXT only accepts 1 and nsd"); } From b669531df5abbfdaf4330fd03eb88aa76fe1a09c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 4 Aug 2026 16:44:10 -0500 Subject: [PATCH 4/9] Simplify signature of nn::gnnb --- Code/Source/solver/all_fun.cpp | 6 ++-- Code/Source/solver/baf_ini.cpp | 14 ++++----- Code/Source/solver/cmm.cpp | 7 ++--- Code/Source/solver/eq_assem.cpp | 19 +++++------- Code/Source/solver/nn.cpp | 54 ++++++++++++++++++++------------- Code/Source/solver/nn.h | 36 +++++++++------------- Code/Source/solver/set_bc.cpp | 21 ++++++------- 7 files changed, 76 insertions(+), 81 deletions(-) diff --git a/Code/Source/solver/all_fun.cpp b/Code/Source/solver/all_fun.cpp index fedc7b3c0..ea6c43a1b 100644 --- a/Code/Source/solver/all_fun.cpp +++ b/Code/Source/solver/all_fun.cpp @@ -803,8 +803,7 @@ double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, if (!isIB) { // Get normal vector in cfg configuration auto Nx = fs.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, insd, fs.eNoN, Nx, n, solutions, cfg, - displacement_index); + n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, cfg, displacement_index); } // Calculating the Jacobian (encodes area of face element) @@ -925,8 +924,7 @@ double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, if (!isIB) { // Get normal vector in cfg configuration auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, lFa.eNoN, Nx, n, solutions, - cfg, displacement_index); + n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, cfg, displacement_index); //CALL GNNB(lFa, e, g, nsd-1, lFa.eNoN, lFa.Nx(:,:,g), n) } else { //CALL GNNIB(lFa, e, g, n) diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index 4d9cf136b..9b4136160 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -528,7 +528,6 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution // Compute integral of normal vector over surface element if (!flag) { - Vector nV(nsd); for (int e = 0; e < lFa.nEl; e++) { if (lFa.eType == ElementType::NRB) { @@ -538,9 +537,9 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, lFa.eNoN, Nx, nV, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + const Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); @@ -802,11 +801,10 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp // CALL NRBNNXB(msh(iM),lFa,e) } for (int g = 0; g < lFa.nG; g++) { - Vector n(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, lFa.eNoN, Nx, n, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + const Vector n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); diff --git a/Code/Source/solver/cmm.cpp b/Code/Source/solver/cmm.cpp index c1c120ffb..0749b4522 100644 --- a/Code/Source/solver/cmm.cpp +++ b/Code/Source/solver/cmm.cpp @@ -280,11 +280,10 @@ void cmm_b(ComMod& com_mod, const faceType& lFa, const int e, const Array nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, 3, Nx, nV, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; diff --git a/Code/Source/solver/eq_assem.cpp b/Code/Source/solver/eq_assem.cpp index 78e2829ec..6d491d11f 100644 --- a/Code/Source/solver/eq_assem.cpp +++ b/Code/Source/solver/eq_assem.cpp @@ -78,11 +78,10 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& } for (int g = 0; g < lFa.nG; g++) { - Vector nV(nsd); auto Nx = lFa.Nx.rslice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoN, Nx, nV, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -253,11 +252,10 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const } // Get surface normal vector - Vector nV(nsd); auto Nx_g = lFa.Nx.rslice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoNb, Nx_g, nV, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx_g, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -329,13 +327,12 @@ void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const S // CALL NRBNNXB(msh(iM),lFa,e) } for (int g = 0; g < lFa.nG; g++) { - Vector n(nsd); auto Nx = lFa.Nx.rslice(g); auto cfg = MechanicalConfigurationType::new_timestep; - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, lFa.eNoN, Nx, n, solutions, - cfg, com_mod.eq[com_mod.cEq].s); + const Vector n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, cfg, + com_mod.eq[com_mod.cEq].s); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); diff --git a/Code/Source/solver/nn.cpp b/Code/Source/solver/nn.cpp index f73dabd30..81b92dea1 100644 --- a/Code/Source/solver/nn.cpp +++ b/Code/Source/solver/nn.cpp @@ -899,17 +899,22 @@ void gnn(const int eNoN, const int nsd, const int insd, Array& Nxi, Arra } } -void gnnb(const ComMod &com_mod, const faceType &lFa, const int e, const int g, - const int nsd, const int insd, const int eNoNb, - const Array &Nx, Vector &n, - const SolutionStates &solutions, - consts::MechanicalConfigurationType cfg, - const unsigned int displacement_index) { +Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const Array &Nx, + const SolutionStates &solutions, + consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index) { // Local aliases for displacement arrays const auto& Dn = solutions.current.get_displacement(); const auto& Do = solutions.old.get_displacement(); auto& cm = com_mod.cm; + const int nsd = com_mod.nsd; + const int insd = Nx.nrows(); + const int eNoNb = Nx.ncols(); + + Vector n(nsd); + #define n_debug_gnnb #ifdef debug_gnnb DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -1034,10 +1039,15 @@ void gnnb(const ComMod &com_mod, const faceType &lFa, const int e, const int g, // Compute adjoining mesh element normal // - Array xXi(nsd,nsd-1); + // Note that here we use the shell element's shape function derivatives + // (msh.Nx), so the relevant intrinsic dimension is that of the shell + // element (nsd - 1, a surface), not that of the boundary edge lFa + // (which is what the local 'insd' holds). + const int msh_insd = nsd - 1; + Array xXi(nsd, msh_insd); for (int a = 0; a < eNoN; a++) { - for (int i = 0; i < insd; i++) { + for (int i = 0; i < msh_insd; i++) { for (int j = 0; j < nsd; j++) { xXi(j,i) = xXi(j,i) + lX(j,a)*msh.Nx(i,a,g); } @@ -1081,7 +1091,7 @@ void gnnb(const ComMod &com_mod, const faceType &lFa, const int e, const int g, n = -n; } - return; + return n; } else { @@ -1097,21 +1107,23 @@ void gnnb(const ComMod &com_mod, const faceType &lFa, const int e, const int g, } n = utils::cross(xXi); - } - // Changing the sign if neccessary. 'a' locates on the face and 'b' - // in the interior of the element. v points outward along ba - // - a = ptr(0); - int b = ptr(lFa.eNoN); - Vector v(nsd); + // Changing the sign if neccessary. 'a' locates on the face and 'b' + // in the interior of the element. v points outward along ba + // + a = ptr(0); + int b = ptr(lFa.eNoN); + Vector v(nsd); - for (int i = 0; i < nsd; i++) { - v(i) = lX(i,a) - lX(i,b); - } + for (int i = 0; i < nsd; i++) { + v(i) = lX(i, a) - lX(i, b); + } + + if (n * v < 0.0) { + n = -n; + } - if (n * v < 0.0) { - n = -n; + return n; } } diff --git a/Code/Source/solver/nn.h b/Code/Source/solver/nn.h index ecf432867..abe23741b 100644 --- a/Code/Source/solver/nn.h +++ b/Code/Source/solver/nn.h @@ -39,11 +39,13 @@ namespace nn { double& Jac, Array& ks); /** - * @brief Return the surface normal vector at a given element and Gauss point. + * @brief Return the area-weighted surface normal at a given element and Gauss + * point. * - * 'g' of face 'lFa' that is the normal weighted by Jac, i.e. - * Jac = norm(n), the Jacobian of the mapping from parent surface element to - * reference/old/new configuration. + * Returns the outward normal at element 'e', Gauss point 'g' of face 'lFa', + * weighted by the surface Jacobian: Jac = norm(n) is the Jacobian of the + * mapping from the parent surface element to the reference/old/new + * configuration. * * For simulations involving structural displacement, this function allows * computing the normal vector in any of the following configurations: @@ -53,22 +55,14 @@ namespace nn { * - old configuration (the mesh is displaced by the displacement field from * previous time step). * - * @todo[michelebucelli] Many of the arguments passed to this function are - * redundant, as they can be inferred from the others. For example, insd and - * eNoNb can be inferred from the size of Nx. Those arguments should be - * removed to simplify the function signature and make it less error-prone. - * * @param[in] com_mod The common module. * @param[in] lFa The boundary face for which the normal vector is computed. * @param[in] e The face-local index of the element for which the normal * vector is computed. * @param[in] g The Gauss point index for which the normal vector is computed. - * @param[in] nsd The number of spatial dimensions. - * @param[in] insd The intrinsic dimension of the boundary surface (typically - * nsd - 1). - * @param[in] eNoNb Number of nodes on a boundary face element. - * @param[in] Nx The shape function derivatives at the Gauss point. - * @param[out] n The computed normal vector. + * @param[in] Nx The shape function derivatives at the Gauss point. Its shape + * (insd x eNoNb) determines the surface's intrinsic dimension and the + * number of nodes per face element. * @param[in] solutions The solution states that the displacement fields are * extracted from. * @param[in] cfg The configuration in which the normal vector is computed @@ -76,13 +70,13 @@ namespace nn { * @param[in] displacement_index The index of the displacement field in the * solution arrays. This should correspond to the start index of the * equation that solves for the displacement. + * @return The area-weighted outward normal vector. */ - void gnnb(const ComMod &com_mod, const faceType &lFa, const int e, - const int g, const int nsd, const int insd, const int eNoNb, - const Array &Nx, Vector &n, - const SolutionStates &solutions, - consts::MechanicalConfigurationType cfg, - const unsigned int displacement_index); + Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const Array &Nx, + const SolutionStates &solutions, + consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); void gnns(const int nsd, const int eNoN, const Array& Nxi, Array& xl, Vector& nV, Array& gCov, Array& gCnv); diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 931114db3..4a8089deb 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1445,11 +1445,10 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const // Gauss integration 1 // for (int g = 0; g < lFa.nG; g++) { - Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoNb, Nx, nV, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -1754,11 +1753,10 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondit lKd = 0.0; for (int g = 0; g < lFa.nG; g++) { - Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoN, Nx, nV, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -2038,11 +2036,10 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons lR = 0.0; for (int g = 0; g < lFa.nG; g++) { - Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd - 1, eNoN, Nx, nV, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + const Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + /* displacement_index = */ 0); double Jac = utils::norm(nV); double w = lFa.w(g)*Jac; N = lFa.N.col(g); From c61a5de768df00e4628b12300037cf95da7ff14b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 4 Aug 2026 17:38:39 -0500 Subject: [PATCH 5/9] Add overloads to gnnb and integ for evaluating in reference configuration without passing meaningless arguments --- Code/Source/solver/all_fun.cpp | 58 +++++++++++++++--------------- Code/Source/solver/all_fun.h | 36 +++++++++++++++++++ Code/Source/solver/baf_ini.cpp | 20 +++-------- Code/Source/solver/cmm.cpp | 4 +-- Code/Source/solver/eq_assem.cpp | 8 ++--- Code/Source/solver/nn.cpp | 10 ++++++ Code/Source/solver/nn.h | 26 ++++++++++++++ Code/Source/solver/ris.cpp | 22 ++++-------- Code/Source/solver/set_bc.cpp | 64 ++++++++++----------------------- Code/Source/solver/txt.cpp | 16 +++------ 10 files changed, 139 insertions(+), 125 deletions(-) diff --git a/Code/Source/solver/all_fun.cpp b/Code/Source/solver/all_fun.cpp index ea6c43a1b..bdb15ab07 100644 --- a/Code/Source/solver/all_fun.cpp +++ b/Code/Source/solver/all_fun.cpp @@ -830,17 +830,16 @@ double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, return result; } -/// @brief This routine integrates vector field s dotted with the face normal n -/// over the face lFa. For example, if s contains the velocity at each node on -/// the face, this function computed the velocity flux through the face. -/// -/// Reproduces 'FUNCTION IntegV(lFa, s)' -/// -/// @param lFa face type, representing a face on the computational mesh -/// @param s an array containing a vector value for each node in the mesh -/// @param pFlag flag for using Taylor-Hood function space for pressure -/// @param cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference. -// +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Vector &s, const SolutionStates &solutions, + bool pFlag) { + // The displacement index is not used in the reference configuration. + constexpr unsigned int unused_displacement_index = 0; + return integ(com_mod, cm_mod, lFa, s, solutions, pFlag, + consts::MechanicalConfigurationType::reference, + unused_displacement_index); +} + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, const Array &s, const SolutionStates &solutions, MechanicalConfigurationType cfg, @@ -959,23 +958,15 @@ double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, return result; } -/// @brief This routine integrate s(l:u,:) over the surface faId, where s is an -/// array of scalars or an array of nsd-vectors. This routine calls overloaded -/// functions to integrate scalars, if s is scalar (i.e. l=u), or vectors if s -/// is vector (i.e. l &s, const SolutionStates &solutions) { + // The displacement index is not used in the reference configuration. + constexpr unsigned int unused_displacement_index = 0; + return integ(com_mod, cm_mod, lFa, s, solutions, + consts::MechanicalConfigurationType::reference, + unused_displacement_index); +} + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, const Array &s, const int l, const SolutionStates &solutions, std::optional uo, @@ -1052,6 +1043,17 @@ double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, return result; } +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const int l, + const SolutionStates &solutions, std::optional uo, + bool THflag) { + // The displacement index is not used in the reference configuration. + constexpr unsigned int unused_displacement_index = 0; + return integ(com_mod, cm_mod, lFa, s, l, solutions, uo, THflag, + consts::MechanicalConfigurationType::reference, + unused_displacement_index); +} + bool is_domain(const ComMod& com_mod, const eqType& eq, const int node, const consts::EquationType phys) { bool result = false; diff --git a/Code/Source/solver/all_fun.h b/Code/Source/solver/all_fun.h index 7dd3637bb..3cf5bb8a8 100644 --- a/Code/Source/solver/all_fun.h +++ b/Code/Source/solver/all_fun.h @@ -72,6 +72,18 @@ namespace all_fun { bool pFlag, consts::MechanicalConfigurationType cfg, const unsigned int displacement_index); + /** + * @brief Integrate a scalar field over a boundary face. + * + * This is the overload to use in the general case. The other overload adds + * the ability to integrate over a displaced configuration, which is only + * relevant for simulations involving structural displacement. See it for the + * meaning of the arguments. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Vector &s, const SolutionStates &solutions, + bool pFlag); + /** * @brief Integrate one or more components of a field over a boundary face. * @@ -119,6 +131,19 @@ namespace all_fun { bool THflag, consts::MechanicalConfigurationType cfg, const unsigned int displacement_index); + /** + * @brief Integrate one or more components of a field over a boundary face. + * + * This is the overload to use in the general case. The other overload adds + * the ability to integrate over a displaced configuration, which is only + * relevant for simulations involving structural displacement. See it for the + * meaning of the arguments. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const int l, + const SolutionStates &solutions, std::optional uo, + bool THflag); + /** * @brief Integrate the flux of a vector field over a boundary face. * @@ -156,6 +181,17 @@ namespace all_fun { consts::MechanicalConfigurationType cfg, const unsigned int displacement_index); + /** + * @brief Integrate the flux of a vector field over a boundary face. + * + * This is the overload to use in the general case. The other overload adds + * the ability to integrate over a displaced configuration, which is only + * relevant for simulations involving structural displacement. See it for the + * meaning of the arguments. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const SolutionStates &solutions); + bool is_domain(const ComMod& com_mod, const eqType& eq, const int node, const consts::EquationType phys); double jacobian(ComMod& com_mod, const int nDim, const int eNoN, const Array& x, const Array&Nxi); diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index 9b4136160..626c21761 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -338,9 +338,7 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l Vector center(3); for (int i = 0; i < nsd; i++) { center(i) = all_fun::integ(com_mod, cm_mod, lFa, com_mod.x, i, solutions, - std::nullopt, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / + std::nullopt, false) / lFa.area; } @@ -450,9 +448,7 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l // double tmp = 1.0; if (btest(lBc.bType, enum_int(BoundaryConditionType::bType_flx))) { - tmp = all_fun::integ(com_mod, cm_mod, lFa, s, solutions, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + tmp = all_fun::integ(com_mod, cm_mod, lFa, s, solutions, false); if (is_zero(tmp)) { tmp = 1.0; throw std::runtime_error("Face '" + lFa.name + "' used for a BC has no non-zero node."); @@ -494,9 +490,7 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution // Vector sA(com_mod.tnNo); sA = 1.0; - double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false); #ifdef debug_face_ini dmsg << "Face '" << lFa.name << "' area: " << area; #endif @@ -537,9 +531,7 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.slice(g); - const Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + const Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); @@ -802,9 +794,7 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp } for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.slice(g); - const Vector n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + const Vector n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); diff --git a/Code/Source/solver/cmm.cpp b/Code/Source/solver/cmm.cpp index 0749b4522..785a3429b 100644 --- a/Code/Source/solver/cmm.cpp +++ b/Code/Source/solver/cmm.cpp @@ -281,9 +281,7 @@ void cmm_b(ComMod& com_mod, const faceType& lFa, const int e, const Array nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; diff --git a/Code/Source/solver/eq_assem.cpp b/Code/Source/solver/eq_assem.cpp index 6d491d11f..463a178ca 100644 --- a/Code/Source/solver/eq_assem.cpp +++ b/Code/Source/solver/eq_assem.cpp @@ -79,9 +79,7 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.rslice(g); - Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -253,9 +251,7 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const // Get surface normal vector auto Nx_g = lFa.Nx.rslice(g); - Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx_g, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx_g, solutions); Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; diff --git a/Code/Source/solver/nn.cpp b/Code/Source/solver/nn.cpp index 81b92dea1..d3e57bb12 100644 --- a/Code/Source/solver/nn.cpp +++ b/Code/Source/solver/nn.cpp @@ -1127,6 +1127,16 @@ Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, } } +Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const Array &Nx, + const SolutionStates &solutions) { + // The displacement index is not used in the reference configuration. + constexpr unsigned int unused_displacement_index = 0; + return gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + unused_displacement_index); +} + /// @brief Compute shell kinematics: normal vector, covariant & contravariant basis vectors /// /// Replicates 'SUBROUTINE GNNS(eNoN, Nxi, xl, nV, gCov, gCnv)' defined in NN.f. diff --git a/Code/Source/solver/nn.h b/Code/Source/solver/nn.h index abe23741b..c8c8910f3 100644 --- a/Code/Source/solver/nn.h +++ b/Code/Source/solver/nn.h @@ -78,6 +78,32 @@ namespace nn { consts::MechanicalConfigurationType cfg, const unsigned int displacement_index); + /** + * @brief Return the area-weighted surface normal at a given element and Gauss + * point. + * + * Returns the outward normal at element 'e', Gauss point 'g' of face 'lFa', + * weighted by the surface Jacobian: Jac = norm(n) is the Jacobian of the + * mapping from the parent surface element to the mesh. + * + * This is the overload to use in the general case. The other overload adds + * the ability to compute the normal in a displaced configuration, which is + * only relevant for simulations involving structural displacement. + * + * @param[in] com_mod The common module. + * @param[in] lFa The boundary face for which the normal vector is computed. + * @param[in] e The face-local index of the element for which the normal + * vector is computed. + * @param[in] g The Gauss point index for which the normal vector is computed. + * @param[in] Nx The shape function derivatives at the Gauss point. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @return The area-weighted outward normal vector. + */ + Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const Array &Nx, + const SolutionStates &solutions); + void gnns(const int nsd, const int eNoN, const Array& Nxi, Array& xl, Vector& nV, Array& gCov, Array& gCnv); diff --git a/Code/Source/solver/ris.cpp b/Code/Source/solver/ris.cpp index 6057eaed5..896fb8bd6 100644 --- a/Code/Source/solver/ris.cpp +++ b/Code/Source/solver/ris.cpp @@ -61,9 +61,7 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) double tmp = msh[iM].fa[iFa].area; RIS.meanP(iProj, i) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, - std::nullopt, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / + std::nullopt, false) / tmp; } } @@ -82,10 +80,8 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) } int iM = RIS.lst(0,0,iProj); int iFa = RIS.lst(0,1,iProj); - RIS.meanFl(iProj) = all_fun::integ( - com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m - 1, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + RIS.meanFl(iProj) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, + 0, solutions, m - 1, false); if (cm.mas(cm_mod)) { std::cout << "For RIS projection: " << iProj << std::endl; @@ -488,13 +484,9 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solution sA = 1.0; lFa = msh[iM].fa[iFa]; // such update may be not correct - tmp_new = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + tmp_new = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false); meanP = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, - m - 1, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / + m - 1, false) / tmp_new; // For the velocity @@ -511,9 +503,7 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solution } meanFl = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, - solutions, m - 1, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + solutions, m - 1, false); std::cout << "The average pressure is: " << meanP << std::endl; std::cout << "The pressure from 0D is: " << eq[cEq].bc[iBc].g << std::endl; diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 4a8089deb..4272ea598 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -166,18 +166,12 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& // Compute avg pressures at 3D Dirichlet boundaries at timesteps n and n+1 else if (utils::btest(bc.bType, iBC_Dir)) { double area = fa.area; - cplBC.fa[ptr].Po = - all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, - std::nullopt, false, - MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / - area; - cplBC.fa[ptr].Pn = - all_fun::integ(com_mod, cm_mod, fa, Yn, nsd, solutions, - std::nullopt, false, - MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / - area; + cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, + solutions, std::nullopt, false) / + area; + cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, fa, Yn, nsd, + solutions, std::nullopt, false) / + area; cplBC.fa[ptr].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; @@ -654,14 +648,10 @@ void rcr_init(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& soluti if (cplBC.initRCR) { auto& fa = com_mod.msh[iM].fa[iFa]; double area = fa.area; - double Qo = - all_fun::integ(com_mod, cm_mod, fa, Yo, 0, solutions, nsd - 1, - false, MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + double Qo = all_fun::integ(com_mod, cm_mod, fa, Yo, 0, solutions, + nsd - 1, false); double Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, - std::nullopt, false, - MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / + std::nullopt, false) / area; cplBC.xo[ptr] = Po - (Qo * cplBC.fa[ptr].RCR.Rp); } else { @@ -806,10 +796,7 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) faceType& lFa = com_mod.msh[iM].fa[iFa]; Vector sA(com_mod.tnNo); sA = 1.0; - double area = - all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false); baf_ini_ns::bc_ini(com_mod, cm_mod, eq.bc[iBc], lFa, solutions); int ptr = bc.cplBCptr; @@ -876,14 +863,12 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, equation_offset + nsd, solutions, std::nullopt, - false, MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / + false) / area; cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, equation_offset + nsd, solutions, std::nullopt, - false, MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / + false) / area; cplBC.fa[ptr].Qo = 0.0; @@ -1446,9 +1431,7 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const // for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.slice(g); - Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -1598,10 +1581,8 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const //h(0) = lBc.g; - double Q_3D = all_fun::integ( - com_mod, cm_mod, lFa, Yn, eq.s, solutions, eq.s + nsd - 1, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + double Q_3D = all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, + eq.s + nsd - 1, false); h(0) = lBc.g; // h(0) = lBc.g - lBc.r * std::abs(Q_3D); @@ -1625,11 +1606,8 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const } } else if (utils::btest(lBc.bType, iBC_res)) { - h(0) = lBc.r * - all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, - eq.s + nsd - 1, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + h(0) = lBc.r * all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, + eq.s + nsd - 1, false); } else if (utils::btest(lBc.bType, iBC_std)) { h(0) = lBc.g; @@ -1754,9 +1732,7 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondit for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.slice(g); - Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -2037,9 +2013,7 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.slice(g); - const Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + const Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); double w = lFa.w(g)*Jac; N = lFa.N.col(g); diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 2347da0e0..6598ef6b8 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -466,29 +466,21 @@ void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eq if (div) { tmp = fa.area; tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, - std::nullopt, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0) / + std::nullopt, false) / tmp; } else { if (pFlag && lTH) { tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, - std::nullopt, true, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + std::nullopt, true); } else { tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, - std::nullopt, false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + std::nullopt, false); } } } else if (m == nsd) { tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, m - 1, - false, - consts::MechanicalConfigurationType::reference, - /* displacement_index = */ 0); + false); } else { throw std::runtime_error("WTXT only accepts 1 and nsd"); } From 6ad3d84ac5956ed15002b353972375e5d451f674 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 5 Aug 2026 16:05:11 -0500 Subject: [PATCH 6/9] Cap displacement is computed taking the coupled equation index into account --- .../solver/CoupledBoundaryCondition.cpp | 38 ++++++++++++------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/Code/Source/solver/CoupledBoundaryCondition.cpp b/Code/Source/solver/CoupledBoundaryCondition.cpp index c147799b6..5785035db 100644 --- a/Code/Source/solver/CoupledBoundaryCondition.cpp +++ b/Code/Source/solver/CoupledBoundaryCondition.cpp @@ -592,7 +592,8 @@ void CoupledBoundaryCondition::initialize_cap(ComMod& com_mod) namespace { /// @brief Gathers cap-node mesh state on the serial rank (columns 0..n_cap-1 match \p cap_gn order). -void gather_global_mesh_state_serial(ComMod& com_mod, const SolutionStates& solutions, bool gather_Y, int nsd, int tnNo, +void gather_global_mesh_state_serial(ComMod& com_mod, const SolutionStates& solutions, bool gather_Y, + int equation_offset, int nsd, int tnNo, const Vector& cap_gn, CapGlobalMeshState& out) { const auto& Do = solutions.old.get_displacement(); @@ -619,14 +620,16 @@ void gather_global_mesh_state_serial(ComMod& com_mod, const SolutionStates& solu continue; } for (int i = 0; i < nsd; i++) { + // x is geometry (rows 0..nsd-1); Do/Dn are solution rows of the + // coupled equation, hence the equation_offset shift. out.x(i, a) = com_mod.x(i, Ac); - out.Do(i, a) = Do(i, Ac); - out.Dn(i, a) = Dn(i, Ac); + out.Do(i, a) = Do(equation_offset + i, Ac); + out.Dn(i, a) = Dn(equation_offset + i, Ac); } if (gather_Y) { for (int i = 0; i < nsd; i++) { - out.Yo(i, a) = Yo(i, Ac); - out.Yn(i, a) = Yn(i, Ac); + out.Yo(i, a) = Yo(equation_offset + i, Ac); + out.Yn(i, a) = Yn(equation_offset + i, Ac); } } break; @@ -637,7 +640,8 @@ void gather_global_mesh_state_serial(ComMod& com_mod, const SolutionStates& solu /// @brief Gathers cap-node mesh state on the MPI root (columns 0..n_cap-1 match \p cap_gn order). void gather_global_mesh_state_parallel(ComMod& com_mod, const CmMod& cm_mod, cmType& cm, const SolutionStates& solutions, - bool gather_Y, int nsd, int tnNo, int root, int nProcs, const Vector& cap_gn, + bool gather_Y, int equation_offset, int nsd, int tnNo, int root, int nProcs, + const Vector& cap_gn, const std::unordered_map& g_to_cap_col, CapGlobalMeshState& out) { @@ -673,21 +677,23 @@ void gather_global_mesh_state_parallel(ComMod& com_mod, const CmMod& cm_mod, cmT continue; } send_buf(idx++) = static_cast(g); + // x is geometry (rows 0..nsd-1); Do/Dn/Yo/Yn are solution rows of the + // coupled equation, hence the equation_offset shift. for (int i = 0; i < nsd; i++) { send_buf(idx++) = com_mod.x(i, Ac); } for (int i = 0; i < nsd; i++) { - send_buf(idx++) = Do(i, Ac); + send_buf(idx++) = Do(equation_offset + i, Ac); } for (int i = 0; i < nsd; i++) { - send_buf(idx++) = Dn(i, Ac); + send_buf(idx++) = Dn(equation_offset + i, Ac); } if (gather_Y) { for (int i = 0; i < nsd; i++) { - send_buf(idx++) = Yo(i, Ac); + send_buf(idx++) = Yo(equation_offset + i, Ac); } for (int i = 0; i < nsd; i++) { - send_buf(idx++) = Yn(i, Ac); + send_buf(idx++) = Yn(equation_offset + i, Ac); } } } @@ -769,13 +775,19 @@ void CoupledBoundaryCondition::gather_global_mesh_state(ComMod& com_mod, const C const int root = cm_mod.master; const int nProcs = cm.np(); + // The coupled equation's velocity and displacement occupy solution rows + // starting at this offset. The gather must read them there rather than + // assuming the coupled equation is the first one (offset 0). This mirrors + // the offset used by compute_flowrates / compute_pressures. + const int equation_offset = com_mod.eq[com_mod.cplBC.equationIndex].s; + if (cm.seq()) { - gather_global_mesh_state_serial(com_mod, solutions, gather_Y, nsd, tnNo, cap_mesh_global_node_ids_, - cap_global_mesh_state_); + gather_global_mesh_state_serial(com_mod, solutions, gather_Y, equation_offset, nsd, tnNo, + cap_mesh_global_node_ids_, cap_global_mesh_state_); return; } - gather_global_mesh_state_parallel(com_mod, cm_mod, cm, solutions, gather_Y, nsd, tnNo, root, nProcs, + gather_global_mesh_state_parallel(com_mod, cm_mod, cm, solutions, gather_Y, equation_offset, nsd, tnNo, root, nProcs, cap_mesh_global_node_ids_, cap_g_to_cap_col_, cap_global_mesh_state_); } From f360c29ae6da183e4f5d3adafa26eddf9a3e8604 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 6 Aug 2026 09:19:48 -0500 Subject: [PATCH 7/9] Finite difference increment for evaluating the 0D solver tangent can be specified in the XML file --- Code/Source/solver/ComMod.h | 13 +++++++++++++ Code/Source/solver/Parameters.cpp | 8 ++++++++ Code/Source/solver/Parameters.h | 14 ++++++++++++++ Code/Source/solver/distribute.cpp | 4 +++- Code/Source/solver/read_files.cpp | 6 ++++++ Code/Source/solver/set_bc.cpp | 13 +++++-------- 6 files changed, 49 insertions(+), 9 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 1d378c2c6..0fecefc45 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -828,6 +828,19 @@ class cplBCType consts::CplBCType schm = consts::CplBCType::cplBC_NA; //int schm = cplBC_NA; + /// @brief Absolute floor on the flow-rate perturbation used to + /// finite-difference the coupled-BC tangent dP/dQ in + /// \c set_bc::calc_der_cpl_bc. This is a dimensional quantity, so it must + /// be set consistently with the unit system of the simulation. + double finite_difference_absolute_perturbation = 1.0e-7; + + /// @brief Flow-rate perturbation used to finite-difference the coupled-BC + /// tangent dP/dQ in \c set_bc::calc_der_cpl_bc, relative to the RMS coupled + /// flow rate. The perturbation actually applied is + /// \c max(rms(Q)*finite_difference_relative_perturbation, + /// \c finite_difference_absolute_perturbation). + double finite_difference_relative_perturbation = 1.0e-5; + /// @brief Path to the 0D code binary file std::string binPath; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 8f9eed3d6..baa8638ee 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1208,6 +1208,14 @@ svZeroDSolverInterfaceParameters::svZeroDSolverInterfaceParameters() { set_parameter("Initial_flows", 0.0, !required, initial_flows); set_parameter("Initial_pressures", 0.0, !required, initial_pressures); + // Finite-difference perturbation used for the coupled-BC tangent dP/dQ, + // diff = max(rms(Q) * Finite_difference_relative_perturbation, + // Finite_difference_absolute_perturbation). + set_parameter("Finite_difference_absolute_perturbation", 1.0e-7, !required, + finite_difference_absolute_perturbation); + set_parameter("Finite_difference_relative_perturbation", 1.0e-5, !required, + finite_difference_relative_perturbation); + set_parameter("Configuration_file", "", required, configuration_file); set_parameter("Shared_library", "", required, shared_library); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index b18bd4661..a1d18d437 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -657,6 +657,17 @@ class CoupleGenBCParameters : public ParameterLists //---------------------------------- // svZeroDSolverInterfaceParameters //---------------------------------- +/// @brief Parameters for coupling to the svZeroDSolver (0D lumped-parameter solver). +/// +/// XML element: \code {.xml} +/// +/// implicit +/// svzerod_3Dcoupling.json +/// /path/to/libsvzero_interface +/// 1.0e-7 +/// 1.0e-5 +/// +/// \endcode // class svZeroDSolverInterfaceParameters : public ParameterLists { @@ -674,6 +685,9 @@ class svZeroDSolverInterfaceParameters : public ParameterLists Parameter initial_flows; Parameter initial_pressures; + Parameter finite_difference_absolute_perturbation; + Parameter finite_difference_relative_perturbation; + Parameter shared_library; bool value_set = false; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 244dd1667..111148266 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -551,8 +551,10 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &cplBC.useGenBC); cm.bcast(cm_mod, &cplBC.useSvZeroD); cm.bcast(cm_mod, &cplBC.useSvOneD); + cm.bcast(cm_mod, &cplBC.finite_difference_absolute_perturbation); + cm.bcast(cm_mod, &cplBC.finite_difference_relative_perturbation); - if (cplBC.useGenBC) { + if (cplBC.useGenBC) { if (cm.slv(cm_mod)) { cplBC.nX = 0; cplBC.xo.resize(cplBC.nX); diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 7c4f770eb..08899995f 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1669,6 +1669,12 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) cplBC.useSvZeroD = true; cplbc_type_str = eq_params->svzerodsolver_interface_parameters.coupling_type.value(); cplBC.svzerod_solver_interface.set_data(eq_params->svzerodsolver_interface_parameters); + cplBC.finite_difference_absolute_perturbation = + eq_params->svzerodsolver_interface_parameters + .finite_difference_absolute_perturbation.value(); + cplBC.finite_difference_relative_perturbation = + eq_params->svzerodsolver_interface_parameters + .finite_difference_relative_perturbation.value(); } if (eq_params->svonedsolver_interface_parameters.defined()) { diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 4272ea598..30086dbb4 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -13,11 +13,12 @@ #include "lhsa.h" #include "mat_fun.h" #include "nn.h" +#include "svOneD_interface.h" +#include "svZeroD_interface.h" #include "ustruct.h" #include "utils.h" +#include #include -#include "svZeroD_interface.h" -#include "svOneD_interface.h" namespace set_bc { @@ -45,10 +46,8 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& #endif const int iEq = com_mod.cplBC.equationIndex; - // NOTE: For coupling with svZeroDPlus, absTol needs to be > 1e-8 to be compatible with the default convergence tolerance of svZeroDPlus (1e-8) - // If this is not true, the finite difference computation of bc.r below results in zero because the perturbation is below the svZeroDPlus tolerance - const double absTol = 1.0e-7; - const double relTol = 1.0e-5; + const double absTol = com_mod.cplBC.finite_difference_absolute_perturbation; + const double relTol = com_mod.cplBC.finite_difference_relative_perturbation; int nsd = com_mod.nsd; auto& eq = com_mod.eq[iEq]; @@ -759,8 +758,6 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) const auto& Yo = solutions.old.get_velocity(); const auto& Do = solutions.old.get_displacement(); - static double absTol = 1.E-8, relTol = 1.E-5; - using namespace consts; const int nsd = com_mod.nsd; From ac087eae5567b2140587a8fc1a6aa5840633762f Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 6 Aug 2026 09:43:10 -0500 Subject: [PATCH 8/9] Write svZeroD output to the same folder as the 3D simulation output --- Code/Source/solver/baf_ini.cpp | 2 +- Code/Source/solver/svZeroD_interface.cpp | 16 ++++++++++++---- Code/Source/solver/svZeroD_interface.h | 16 +++++++++++++++- 3 files changed, 28 insertions(+), 6 deletions(-) diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index 626c21761..4ac0e82f5 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -163,7 +163,7 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) } if (com_mod.cplBC.useSvZeroD) { - svZeroD::init_svZeroD(com_mod, cm_mod); + svZeroD::init_svZeroD(com_mod, cm_mod, simulation->get_chnl_mod().appPath); } if (com_mod.cplBC.useSvOneD) { diff --git a/Code/Source/solver/svZeroD_interface.cpp b/Code/Source/solver/svZeroD_interface.cpp index 2a3c8da8f..7a7b0caf6 100644 --- a/Code/Source/solver/svZeroD_interface.cpp +++ b/Code/Source/solver/svZeroD_interface.cpp @@ -58,6 +58,11 @@ static int numCoupledSrfs; static bool writeSvZeroD = true; static double svZeroDTime = 0.0; +/// Directory the svZeroD output files (svZeroD_data, Q_svZeroD, P_svZeroD) are +/// written to, set in init_svZeroD() to the simulation results directory +/// (chnl_mod.appPath). Includes the trailing separator. +static std::string svZeroD_output_dir = ""; + int num_output_steps; int system_size; int model_id; @@ -148,7 +153,7 @@ void write_svZeroD_solution(const double* lpn_time, std::vector& lpn_sol std::vector variable_names; variable_names = interface->variable_names_; std::ofstream out_file; - out_file.open("svZeroD_data", std::ios::out | std::ios::app); + out_file.open(svZeroD_output_dir + "svZeroD_data", std::ios::out | std::ios::app); out_file<(variable_names[i])<<" "; @@ -156,7 +161,7 @@ void write_svZeroD_solution(const double* lpn_time, std::vector& lpn_sol out_file<<'\n'; } else { std::ofstream out_file; - out_file.open("svZeroD_data", std::ios::out | std::ios::app); + out_file.open(svZeroD_output_dir + "svZeroD_data", std::ios::out | std::ios::app); out_file<<*lpn_time<<" "; for (int i = 0; i < system_size; i++) { out_file<& surfID, double Q[], double P[]) { int nParam = 2; - const char* fileNames[2] = {"Q_svZeroD", "P_svZeroD"}; + const std::string fileNames[2] = {svZeroD_output_dir + "Q_svZeroD", + svZeroD_output_dir + "P_svZeroD"}; std::vector> R(nParam, std::vector(*nSrfs)); if (*nSrfs == 0) return; @@ -226,7 +232,7 @@ void print_svZeroD(int* nSrfs, const std::vector& surfID, double Q[], doubl // init_svZeroD //-------------- // -void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) +void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod, const std::string& appPath) { using namespace consts; @@ -241,6 +247,8 @@ void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) auto& cm = com_mod.cm; double dt = com_mod.dt; + svZeroD_output_dir = appPath; + build_svzero_coupled_bc_idxs(com_mod); if (cplBC.nSvZeroD_coupled_bc == 0) { throw std::runtime_error( diff --git a/Code/Source/solver/svZeroD_interface.h b/Code/Source/solver/svZeroD_interface.h index 1de624f85..7187b7542 100644 --- a/Code/Source/solver/svZeroD_interface.h +++ b/Code/Source/solver/svZeroD_interface.h @@ -17,7 +17,21 @@ void get_coupled_QP(ComMod& com_mod, double QCoupled[], double QnCoupled[], doub void print_svZeroD(int* nSrfs, const std::vector& surfID, double Q[], double P[]); -void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod); +/** + * @brief Set up the svZeroD model and its coupled boundary conditions. + * + * Builds the list of svZeroD-coupled boundaries, creates the svZeroD model + * from the interface data in \c com_mod.cplBC, applies the initial flows and + * pressures, and writes the header of the svZeroD state output file. + * + * @param[in,out] com_mod Simulation data; the coupled boundary conditions and + * \c cplBC bookkeeping are initialized here. + * @param[in] cm_mod MPI communicator data. + * @param[in] appPath Directory the simulation results are written to. The + * svZeroD output files (svZeroD_data, Q_svZeroD, P_svZeroD) are written + * there as well. + */ +void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod, const std::string& appPath); void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag); From b76e4b059619e0e2c045b63c6b6c27bb42191d54 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 19 Aug 2026 13:48:37 -0500 Subject: [PATCH 9/9] svZeroDSolver output file is overridden by every simulation, instead of appended to --- Code/Source/solver/svZeroD_interface.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Source/solver/svZeroD_interface.cpp b/Code/Source/solver/svZeroD_interface.cpp index 7a7b0caf6..ad791f5c7 100644 --- a/Code/Source/solver/svZeroD_interface.cpp +++ b/Code/Source/solver/svZeroD_interface.cpp @@ -153,7 +153,7 @@ void write_svZeroD_solution(const double* lpn_time, std::vector& lpn_sol std::vector variable_names; variable_names = interface->variable_names_; std::ofstream out_file; - out_file.open(svZeroD_output_dir + "svZeroD_data", std::ios::out | std::ios::app); + out_file.open(svZeroD_output_dir + "svZeroD_data", std::ios::out); out_file<(variable_names[i])<<" ";