diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 7db495dae..423b38316 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -3,17 +3,71 @@ #include "ActiveStress.h" +#include "mat_fun.h" + bool supports_active_stress(const consts::EquationType eq_type) { return eq_type == consts::EquationType::phys_struct || eq_type == consts::EquationType::phys_ustruct || eq_type == consts::EquationType::phys_FSI; } +void ActiveStress::Evaluator::update(const ActiveStress &active_stress, + const Vector &nodes) { + active_stress_ = &active_stress; + + const unsigned int n_states = active_stress.n_states; + + state_.resize(n_states, nodes.size()); + + // Vector::resize reallocates even when the size is unchanged. + if (interpolated_state.size() != n_states) + interpolated_state.resize(n_states); + + // Friend access to active_stress.states, so that gathering the state of an + // element does not need to go through an accessor. + for (int a = 0; a < nodes.size(); ++a) + for (unsigned int j = 0; j < n_states; ++j) + state_(j, a) = active_stress.states(j, nodes(a)); +} + +template +ActiveStress::ActiveTension ActiveStress::Evaluator::evaluate( + const Eigen::Ref &N, const mat_fun::Matrix &F, + const Eigen::Ref> &fN) + const { + if (active_stress_ == nullptr) + return {}; + + // The fiber stretch is only computed for the models that use it. Besides + // saving the work, this keeps the other models usable on a mesh with no fiber + // directions, where fN is zero and the stretch would come out zero too. + double fiber_stretch = 1.0; + + if (active_stress_->needs_fiber_stretch()) + fiber_stretch = (F * fN.col(0)).norm(); + + // Interpolate the nodal state to the quadrature point. + Eigen::Map(interpolated_state.data(), + interpolated_state.size()) = state_ * N; + + return active_stress_->compute_tension(interpolated_state, fiber_stretch); +} + +template ActiveStress::ActiveTension ActiveStress::Evaluator::evaluate<2>( + const Eigen::Ref &, const mat_fun::Matrix<2> &, + const Eigen::Ref> &) const; + +template ActiveStress::ActiveTension ActiveStress::Evaluator::evaluate<3>( + const Eigen::Ref &, const mat_fun::Matrix<3> &, + const Eigen::Ref> &) const; + void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { eta_f = params.get_eta_f(); eta_s = params.get_eta_s(); eta_n = params.get_eta_n(); + implicit_state_coupling_ = params.get_implicit_state_coupling(); + read_model_specific_parameters( params.get_parameters(params.get_model_name())); } @@ -24,6 +78,8 @@ void ActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &eta_s); cm.bcast(cm_mod, &eta_n); + cm.bcast(cm_mod, &implicit_state_coupling_); + distribute_model_specific_parameters(cm_mod, cm); } @@ -39,21 +95,23 @@ void ActiveStress::init(const unsigned int tnNo) { states(j, i) = state_loc(j); } - active_tension.resize(tnNo); + states_at_time_step_start.resize(n_states, tnNo); + states_at_time_step_start = states; } -void ActiveStress::advance_time_step(const double t, const double dt, - const Vector &calcium, - const Vector &fiber_stretch, - const Vector &fiber_stretch_rate) { +void ActiveStress::time_advance() { states_at_time_step_start = states; } + +void ActiveStress::update(const double t, const double dt, + const Vector &calcium, + const Vector &fiber_stretch, + const Vector &fiber_stretch_rate) { time = t; - for (unsigned int i = 0; i < states.ncols(); ++i) { - Vector state_loc = states.col(i); + // Advance the state from the beginning of the time step. + for (int i = 0; i < states.ncols(); ++i) { + Vector state_loc = states_at_time_step_start.col(i); advance_time_step_local(t, dt, calcium[i], fiber_stretch[i], fiber_stretch_rate[i], state_loc); states.set_col(i, state_loc); - - active_tension[i] = compute_active_tension_local(state_loc, fiber_stretch[i]); } } \ No newline at end of file diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index 29a419cf6..b90f0bc40 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -9,6 +9,7 @@ #include "Vector.h" #include "consts.h" #include "factory.h" +#include "mat_fun.h" #include "CmMod.h" @@ -23,6 +24,15 @@ bool supports_active_stress(const consts::EquationType eq_type); /** * @brief Abstract active stress class. * + * ## Table of contents + * + * - @ref activestress-overview + * - @ref activestress-directions + * - @ref activestress-implementing + * - @ref activestress-coupling + * + * ## Overview {#activestress-overview} + * * This class provides an interface for defining active stress models, i.e. * models that, in the context of structural mechanics of muscular tissue, * compute an active tension representing the contribution of muscular @@ -41,11 +51,11 @@ bool supports_active_stress(const consts::EquationType eq_type); * The expression assumed above implies that the active tension is a local * function of the variables it depends on, that is the active tension at a * given point only depends on the value of other variables at that same point. - * Accordingly, this class works nodally, by evaluating the active tension at - * every mesh node and storing it in a vector, whose values can be accessed - * through @ref ActiveStress::get_tension_fibers. + * Accordingly, this class works nodally, storing the state of contraction at + * every mesh node and evaluating the active tension from it, on demand, + * through @ref ActiveStress::compute_tension. * - * ### Directional distribution of active stress + * ## Directional distribution of active stress {#activestress-directions} * * In muscular mechanics models, active stress normally acts only along the * direction of fibers @f$\fiberdirection@f$, reflecting the fact that @@ -69,20 +79,20 @@ bool supports_active_stress(const consts::EquationType eq_type); * @f$\eta_f + \eta_s + \eta_n = 1@f$. * * This class stores the values of @f$\eta_f@f$, @f$\eta_s@f$ and @f$\eta_n@f$, - * and provides the functions @ref ActiveStress::get_tension_fibers, - * @ref ActiveStress::get_tension_sheets and @ref - * ActiveStress::get_tension_sheet_normals to access @f$\eta_f \Tact@f$, - * @f$\eta_s \Tact@f$ and @f$\eta_n \Tact@f$, respectively. + * and applies them in @ref ActiveStress::compute_tension, which returns + * @f$\eta_f \Tact@f$, @f$\eta_s \Tact@f$ and @f$\eta_n \Tact@f$, bundled in an + * @ref ActiveTension. * - * ### Implementing concrete active stress models + * ## Implementing concrete active stress models {#activestress-implementing} * * To implement a new active stress model, the following steps need to be taken: * * 1. Create a new class derived from @ref ActiveStress. - * 2. Override the methods @ref init_local, @ref advance_time_step_local and - * @ref compute_active_tension_local, defining the initial condition, - * time evolution and active tension computation, respectively, for a single - * node. + * 2. Override the methods @ref init_local, @ref advance_time_step_local, + * @ref compute_active_tension_local and + * @ref compute_active_tension_derivative_local, defining the initial + * condition, time evolution, active tension and its partial derivative with + * respect to the fiber stretch, respectively, for a single node. * 3. Create a new class derived from @ref ActiveStressModelParameters to store * the parameters specific to the new active stress model. * 4. Override the methods @ref get_parameters, @@ -96,19 +106,146 @@ bool supports_active_stress(const consts::EquationType eq_type); * Notice that if the model is expressed in terms of a system of ODEs, it can * be implemented by deriving from @ref ActiveStressODE, which already addresses * some of the points above. + * + * ## Coupling with the structural mechanics problem {#activestress-coupling} + * + * The active tension depends on the fiber stretch @f$\fiberstretch@f$ both + * directly, through the expression of @f$\Tact@f$, and indirectly, through the + * state @f$\astressstate@f$, which is itself driven by the fiber stretch. The + * mechanics problem, in turn, depends on @f$\Tact@f$. + * + * Explicit time discretization for the direct dependence was observed to lead + * to instabilities. Accordingly, the direct dependence is discretized + * implicitly, that is the active tension is recomputed within the nonlinear + * iterations for the structure problem that uses ActiveStress. To facilitate + * the convergence of nonlinear iterations, this class also allows to compute + * the derivative @f$\frac{\partial\Tact}{\partial\fiberstretch}@f$, which is + * used to assemble tangent terms associated to this in the structural system. + * + * Indirect dependence was not observed to give rise to instabilities. + * Accordingly, it is discretized explicitly by default, meaning that the state + * is updated once per time step evaluating the fiber stretch + * @f$\fiberstretch@f$ using the displacement from the previous time step. The + * user can change this behavior by setting the parameter @c + * Implicit_state_coupling to @c true in the XML file. This will make the state + * update every nonlinear iteration. No tangent terms are computed for this + * contribution, so nonlinear iterations can be expected to converge more slowly + * when this is enabled. */ class ActiveStress { public: + /** + * @brief Active tension information at a point. + * + * This struct bundles the active tension along the three principal + * directions (fibers @f$\fiberdirection@f$, sheets @f$\sheetdirection@f$ and + * sheet normals @f$\sheetnormaldirection@f$) and their partial derivatives + * with respect to the fiber stretch. + * + * It is a convenience data structure used to pass this information to + * functions that consume active tension information (e.g. the structural + * mechanics assembly functions). + */ + struct ActiveTension { + /// Tension along the fiber direction, @f$\eta_f \Tact@f$. + double fibers = 0.0; + + /// Tension along the sheet direction, @f$\eta_s \Tact@f$. + double sheets = 0.0; + + /// Tension along the sheet-normal direction, @f$\eta_n \Tact@f$. + double sheet_normals = 0.0; + + /// Derivative of @ref fibers with respect to the fiber stretch, at fixed + /// state, @f$\eta_f \pdv*{\Tact}{\fiberstretch}@f$. + double d_fibers = 0.0; + + /// Derivative of @ref sheets with respect to the fiber stretch, at fixed + /// state, @f$\eta_s \pdv*{\Tact}{\fiberstretch}@f$. + double d_sheets = 0.0; + + /// Derivative of @ref sheet_normals with respect to the fiber stretch, at + /// fixed state, @f$\eta_n \pdv*{\Tact}{\fiberstretch}@f$. + double d_sheet_normals = 0.0; + }; + + /** + * @brief Helper to evaluate active tension at quadrature points. + * + * The active stress model stores the state at the degrees of freedom, but the + * active tension needs to be evaluated at quadrature points. This class + * allows to do that by interpolating the state to quadrature points and then + * evaluating the active tension there. + */ + class Evaluator { + public: + /** + * @brief Update the evaluator for a given element. + * + * @param[in] active_stress Active stress model of the domain the element + * belongs to. + * @param[in] nodes Indices of the mesh nodes of the element. + */ + void update(const ActiveStress &active_stress, const Vector &nodes); + + /** + * @brief Reset this evaluator so that @ref evaluate returns zero tension, + * until the next call to @ref update. + * + * Used for elements whose domain has no active stress model. + */ + void clear() { active_stress_ = nullptr; } + + /** + * @brief Compute the active tension at a quadrature point. + * + * @param[in] N Shape functions at the quadrature point, of the same nodes + * the state was gathered at by @ref update. + * @param[in] F Deformation gradient at the quadrature point. + * @param[in] fN Fiber directions of the element, one per column, the first + * column being the fiber direction itself. + * + * @return Active tension along fibers, sheets and sheet normals, and their + * derivatives with respect to the fiber stretch, bundled in an object of + * type @ref ActiveTension. + * + * @tparam nsd Number of spatial dimensions. + */ + template + ActiveTension + evaluate(const Eigen::Ref &N, + const mat_fun::Matrix &F, + const Eigen::Ref> + &fN) const; + + private: + /// Active stress model of the domain the element belongs to, or null if + /// @ref clear was called last, or if this evaluator was never updated. + const ActiveStress *active_stress_ = nullptr; + + /// State variables at the element nodes, of size (n_states, element + /// nodes). + Eigen::MatrixXd state_; + + /// State variables interpolated at the quadrature point, of size + /// n_states. Scratch storage for @ref evaluate, mutable so that it can + /// stay const. + mutable Vector interpolated_state; + }; + + /// Grants @ref Evaluator direct access to @ref states. + friend class Evaluator; + /** * @brief Constructor. * * @param n_states_ Number of state variables for this model. * @param needs_fiber_stretch Whether this model uses the fiber stretch - * passed to @ref advance_time_step. This flag can be used to determine - * whether fiber stretch computation can be skipped for efficiency. + * passed to @ref update. This flag can be used to determine whether fiber + * stretch computation can be skipped for efficiency. * @param needs_fiber_stretch_rate Whether this model uses the fiber stretch - * rate passed to @ref advance_time_step. This flag can be used to determine - * whether fiber stretch rate computation can be skipped for efficiency. + * rate passed to @ref update. This flag can be used to determine whether + * fiber stretch rate computation can be skipped for efficiency. */ ActiveStress(const unsigned int n_states_, const bool needs_fiber_stretch, const bool needs_fiber_stretch_rate) @@ -137,25 +274,32 @@ class ActiveStress { void distribute_parameters(const CmMod &cm_mod, const cmType &cm); /** - * @brief Get the tension along fibers @f$\eta_f \Tact@f$ at a given point. - */ - double get_tension_fibers(const int idx) const { - return eta_f * active_tension[idx]; - } - - /** - * @brief Get the tension along sheets @f$\eta_s \Tact@f$ at a given point. + * @brief Get the state vector at a given mesh node. + * + * @param[in] Ac Index of the mesh node. */ - double get_tension_sheets(const int idx) const { - return eta_s * active_tension[idx]; + Vector get_state(const unsigned int Ac) const { + return states.col(Ac); } /** - * @brief Get the tension along sheet normals @f$\eta_n \Tact@f$ at a given - * point. + * @brief Compute the active tension given the state vector and fiber stretch. + * + * @param[in] state State vector at the point. + * @param[in] fiber_stretch Fiber stretch at the point. + * + * @return Active tension along fibers, sheets and sheet normals, and their + * derivatives with respect to the fiber stretch, bundled in an object of + * type @ref ActiveTension. */ - double get_tension_sheet_normals(const int idx) const { - return eta_n * active_tension[idx]; + ActiveTension compute_tension(const Vector &state, + const double fiber_stretch) const { + const double tension = compute_active_tension_local(state, fiber_stretch); + const double derivative = + compute_active_tension_derivative_local(state, fiber_stretch); + + return {eta_f * tension, eta_s * tension, eta_n * tension, + eta_f * derivative, eta_s * derivative, eta_n * derivative}; } /** @@ -169,7 +313,22 @@ class ActiveStress { virtual void init(const unsigned int tnNo); /** - * @brief Advance in time. + * @brief Begin a new time step. + * + * Stores the current state as the initial condition of the time step. Must be + * called once per time step, before any call to @ref update. + */ + virtual void time_advance(); + + /** + * @brief Update the state over the current time step. + * + * Advances the state stored by @ref time_advance over one time step, using + * the given calcium, fiber stretch and fiber stretch rate. + * + * This function may be called more than once per time step: every call + * restarts from the state stored by @ref time_advance, so the resulting state + * depends only on the arguments of the last call. * * @param[in] t Current time (i.e. the time instant being advanced to). * @param[in] dt Time step size. @@ -179,39 +338,36 @@ class ActiveStress { * @param[in] fiber_stretch_rate Fiber stretch rate at every node. This is * usually computed with post::fib_stretch_rate. */ - virtual void advance_time_step(const double t, const double dt, - const Vector &calcium, - const Vector &fiber_stretch, - const Vector &fiber_stretch_rate); - - /// Number of state variables for this model. - const unsigned int n_states; + virtual void update(const double t, const double dt, + const Vector &calcium, + const Vector &fiber_stretch, + const Vector &fiber_stretch_rate); /** - * @brief Whether this model uses the fiber stretch passed to - * @ref advance_time_step. This flag can be used to determine whether fiber - * stretch computation can be skipped for efficiency. + * @brief Whether this model uses the fiber stretch passed to @ref update. + * This flag can be used to determine whether fiber stretch computation can be + * skipped for efficiency. */ bool needs_fiber_stretch() const { return needs_fiber_stretch_; } /** * @brief Whether this model uses the fiber stretch rate passed to - * @ref advance_time_step. This flag can be used to determine whether fiber - * stretch rate computation can be skipped for efficiency. + * @ref update. This flag can be used to determine whether fiber stretch rate + * computation can be skipped for efficiency. */ bool needs_fiber_stretch_rate() const { return needs_fiber_stretch_rate_; } -protected: /** - * @brief Backing store for @ref needs_fiber_stretch. + * @brief Whether the state of this model is updated within the nonlinear + * iterations of the mechanics problem, i.e. whether the indirect dependence + * of the active tension on the fiber stretch is treated implicitly. */ - bool needs_fiber_stretch_; + bool implicit_state_coupling() const { return implicit_state_coupling_; } - /** - * @brief Backing store for @ref needs_fiber_stretch_rate. - */ - bool needs_fiber_stretch_rate_; + /// Number of state variables for this model. + const unsigned int n_states; +protected: /** * @brief Read model parameters from a parameter object. * @@ -266,14 +422,46 @@ class ActiveStress { compute_active_tension_local(const Vector &state, const double fiber_stretch) const = 0; - /// Current time. Updated whenever calling @ref advance_time_step. - double time; + /** + * @brief Compute the partial derivative of the active tension with respect to + * the fiber stretch, at fixed state, for a single node. + * + * @param[in] state State vector for a single node. + * @param[in] fiber_stretch Fiber stretch at the current node. + */ + virtual double + compute_active_tension_derivative_local(const Vector &state, + const double fiber_stretch) const = 0; + + /** + * @brief Backing store for @ref needs_fiber_stretch. + */ + bool needs_fiber_stretch_; + + /** + * @brief Backing store for @ref needs_fiber_stretch_rate. + */ + bool needs_fiber_stretch_rate_; + + /// Time instant being advanced to. Set by @ref update. + double time = 0.0; /// State variables for the model. Array states; - /// Active tension at every node. - Vector active_tension; + /** + * @brief State variables at the beginning of the current time step. + * + * Set by @ref time_advance and used by @ref update as the initial condition + * of every call within the time step. + */ + Array states_at_time_step_start; + + /** + * @brief Whether the state of this model is updated within the nonlinear + * iterations of the mechanics problem. + */ + bool implicit_state_coupling_; /// Active tension coefficient along the fiber direction. double eta_f; diff --git a/Code/Source/solver/ActiveStressNashPanfilov.h b/Code/Source/solver/ActiveStressNashPanfilov.h index 2f961560b..e00422c95 100644 --- a/Code/Source/solver/ActiveStressNashPanfilov.h +++ b/Code/Source/solver/ActiveStressNashPanfilov.h @@ -113,6 +113,18 @@ class ActiveStressNashPanfilov : public ActiveStressODE { compute_active_tension_local(const Vector &state, const double fiber_stretch) const override; + /** + * @brief Compute the partial derivative of the active tension with respect + * to the fiber stretch, at fixed state, for a single node. + * + * The active tension does not depend on the fiber stretch, so this is zero. + */ + virtual double compute_active_tension_derivative_local( + const Vector &state, + const double fiber_stretch) const override { + return 0.0; + } + /// @name Model parameters. /// @{ diff --git a/Code/Source/solver/ActiveStressRegazzoni.cpp b/Code/Source/solver/ActiveStressRegazzoni.cpp index 3c1dea8ed..8109d795b 100644 --- a/Code/Source/solver/ActiveStressRegazzoni.cpp +++ b/Code/Source/solver/ActiveStressRegazzoni.cpp @@ -152,6 +152,14 @@ double ActiveStressRegazzoni::compute_active_tension_local( fraction_single_overlap(sarcomere_length); } +double ActiveStressRegazzoni::compute_active_tension_derivative_local( + const Vector &state, const double fiber_stretch) const { + const double sarcomere_length = SL0 * fiber_stretch; + + return a_XB * (state[xb_index(1)] + state[xb_index(3)]) * SL0 * + fraction_single_overlap_derivative(sarcomere_length); +} + ActiveStressRegazzoni::RUArray ActiveStressRegazzoni::ru_transition_rates_tropomyosin() const { RUArray rates_T; @@ -322,4 +330,20 @@ ActiveStressRegazzoni::fraction_single_overlap(double sarcomere_length) const { return 0.0; } +double ActiveStressRegazzoni::fraction_single_overlap_derivative( + double sarcomere_length) const { + const double SL = sarcomere_length; + const double half_single_overlap = (LM - LB) * 0.5; + + if (SL > LA && SL <= LM) + return 1.0 / half_single_overlap; + if (SL > LM && SL <= 2.0 * LA - LB) + return 0.5 / half_single_overlap; + if (SL > 2.0 * LA - LB && SL <= 2.0 * LA + LB) + return 0.0; + if (SL > 2.0 * LA + LB && SL <= 2.0 * LA + LM) + return -0.5 / half_single_overlap; + return 0.0; +} + REGISTER_ACTIVE_STRESS_MODEL("Regazzoni", ActiveStressRegazzoni); diff --git a/Code/Source/solver/ActiveStressRegazzoni.h b/Code/Source/solver/ActiveStressRegazzoni.h index 331084533..e1127e782 100644 --- a/Code/Source/solver/ActiveStressRegazzoni.h +++ b/Code/Source/solver/ActiveStressRegazzoni.h @@ -45,8 +45,9 @@ * @c ActiveStress rather than @c ActiveStressODE because it requires a * customized time-stepping scheme to handle the stiffness of the model. * - * @todo Force-strain-rate feedback requires a stabilization strategy for robust - * use in coupled electromechanics. This will be addressed in a follow-up PR. + * @todo[michelebucelli] Force-strain-rate feedback requires a stabilization + * strategy for robust use in coupled electromechanics. This will be addressed + * in a follow-up PR. */ class ActiveStressRegazzoni : public ActiveStress { public: @@ -190,6 +191,22 @@ class ActiveStressRegazzoni : public ActiveStress { compute_active_tension_local(const Vector &state, const double fiber_stretch) const override; + /** + * @brief Compute the derivative of the active tension with respect to the + * fiber stretch, at fixed state, for a single node. + * + * The state enters @ref compute_active_tension_local as a factor, so the + * derivative is that same factor times the derivative of the single-overlap + * fraction, chained through @f$SL = SL_0 \, \fiberstretch@f$: + * @f[ + * \pdv{\Tact}{\fiberstretch} = a_\text{XB} + * \left(\mu_P^1 + \mu_N^1\right) SL_0 \, \phi'(SL)\;. + * @f] + */ + virtual double compute_active_tension_derivative_local( + const Vector &state, + const double fiber_stretch) const override; + private: /// Array indexed over the four binary RU configuration variables (TL, TC, TR, /// CC). @@ -270,6 +287,19 @@ class ActiveStressRegazzoni : public ActiveStress { */ double fraction_single_overlap(double sarcomere_length) const; + /** + * @brief Derivative of the single-overlap fraction with respect to the + * sarcomere length. + * + * @ref fraction_single_overlap is piecewise linear, so this is piecewise + * constant and jumps at the ends of its pieces. It is taken to be the + * derivative from the left there, matching the intervals of + * @ref fraction_single_overlap. + * + * @param[in] sarcomere_length Sarcomere length @f$SL@f$ [length]. + */ + double fraction_single_overlap_derivative(double sarcomere_length) const; + /// @} /// @name RU model parameters diff --git a/Code/Source/solver/ActiveStressUniformSteady.h b/Code/Source/solver/ActiveStressUniformSteady.h index a06b67855..513ab6fbf 100644 --- a/Code/Source/solver/ActiveStressUniformSteady.h +++ b/Code/Source/solver/ActiveStressUniformSteady.h @@ -86,6 +86,18 @@ class ActiveStressUniformSteady : public ActiveStress { return value; } + /** + * @brief Compute the partial derivative of the active tension with respect + * to the fiber stretch, at fixed state, for a single node. + * + * The active tension does not depend on the fiber stretch, so this is zero. + */ + virtual double compute_active_tension_derivative_local( + const Vector &state, + const double fiber_stretch) const override { + return 0.0; + } + /// Active tension value. double value; }; diff --git a/Code/Source/solver/ActiveStressUniformUnsteady.h b/Code/Source/solver/ActiveStressUniformUnsteady.h index 31733ee00..c9f74751b 100644 --- a/Code/Source/solver/ActiveStressUniformUnsteady.h +++ b/Code/Source/solver/ActiveStressUniformUnsteady.h @@ -55,6 +55,8 @@ class ActiveStressUniformUnsteady : public ActiveStress { * * Calls the parent class initialization method, and reads the Fourier * coefficient from file. + * + * @param[in] tnNo Total number of mesh nodes for the current rank. */ virtual void init(const unsigned int tnNo) override; @@ -96,6 +98,18 @@ class ActiveStressUniformUnsteady : public ActiveStress { compute_active_tension_local(const Vector &state, const double fiber_stretch) const override; + /** + * @brief Compute the partial derivative of the active tension with respect + * to the fiber stretch, at fixed state, for a single node. + * + * The active tension does not depend on the fiber stretch, so this is zero. + */ + virtual double compute_active_tension_derivative_local( + const Vector &state, + const double fiber_stretch) const override { + return 0.0; + } + /// Toggle between ramp or Fourier transform. bool ramp; diff --git a/Code/Source/solver/ArtificialNeuralNetMaterial.cpp b/Code/Source/solver/ArtificialNeuralNetMaterial.cpp index b9f3a8c8c..ba2756aea 100644 --- a/Code/Source/solver/ArtificialNeuralNetMaterial.cpp +++ b/Code/Source/solver/ArtificialNeuralNetMaterial.cpp @@ -111,93 +111,116 @@ void ArtificialNeuralNetMaterial::evaluate(const double aInv[9], double &psi, do } } -template +template void ArtificialNeuralNetMaterial::computeInvariantsAndDerivatives( -const Matrix& C, const Matrix& fl, int nfd, double J2d, double J4d, const Matrix& Ci, -const Matrix& Idm, const double Tfa, Matrix& N1, double& psi, double (&Inv)[9],std::array,9>& dInv, -std::array,9>& ddInv) const { - - // double Inv[9] = {0}; - Matrix C2 = C * C; - - Inv[0] = J2d * C.trace(); - Inv[1] = 0.50 * (Inv[0]*Inv[0] - J4d * (C*C).trace()); - Inv[2] = C.determinant(); - Inv[3] = J2d * (fl.col(0).dot(C * fl.col(0))); - Inv[4] = J4d * (fl.col(0).dot(C2 * fl.col(0))); - - Matrix dInv1 = -Inv[0]/3 * Ci + J2d * Idm; - Matrix dInv2 = (C2.trace()/3)*Ci + Inv[0]*dInv1 + J4d*C; - Matrix dInv3 = Inv[2]*Ci; - N1 = fl.col(0)*fl.col(0).transpose(); - Matrix dInv4 = -Inv[3]/3*Ci + J2d*N1; - Matrix dInv5 = J4d*(N1*C + C*N1) - Inv[4]/3*Ci; - - Matrix dInv6, dInv7, dInv8, dInv9; - Tensor ddInv6, ddInv7, ddInv8, ddInv9; - // initialize to 0 - dInv6.setZero(); - dInv7.setZero(); - dInv8.setZero(); - dInv9.setZero(); - ddInv6.setZero(); - ddInv7.setZero(); - ddInv8.setZero(); - ddInv9.setZero(); - - Tensor dCidC = -symmetric_dyadic_product(Ci, Ci); - Matrix dJ4ddC = -2.0/3.0 * J4d * Ci; - - Tensor ddInv1 = (-1.0/3.0)*(dyadic_product(dInv1,Ci) + Inv[0]*dCidC + J2d*dyadic_product(Ci,Idm)); - Tensor ddInv2 = dyadic_product(dInv1,dInv1) + Inv[0]*ddInv1 + (1.0/3.0)*C2.trace()*dCidC - + (1.0/3.0)*dyadic_product((C2.trace()*dJ4ddC + 2*J4d*C),Ci) - + dyadic_product(dJ4ddC,C) - J4d*fourth_order_identity(); - Tensor ddInv3 = dyadic_product(dInv3,Ci) + Inv[2]*dCidC; - Tensor ddInv4 = (-1.0/3.0)*(dyadic_product(dInv4,Ci) + J2d*dyadic_product(Ci,N1) + Inv[3]*dCidC); - Matrix sum1 = (N1*C + C*N1); - Tensor ddInv5 = (-1.0/3.0)*(dyadic_product(dInv5,Ci) + Inv[4]*dCidC + 2*J4d*dyadic_product(Ci,sum1)) - + J4d*(2*symmetric_dyadic_product(N1,Idm) - dyadic_product(N1,Idm) - + 2*symmetric_dyadic_product(Idm,N1) - dyadic_product(Idm,N1)); - - if (nfd == 2) { - Inv[5] = J2d * (fl.col(0).dot(C * fl.col(1))); - Inv[6] = J4d * (fl.col(0).dot(C2 * fl.col(1))); - Inv[7] = J2d * (fl.col(1).dot(C * fl.col(1))); - Inv[8] = J4d * (fl.col(1).dot(C2 * fl.col(1))); - - Matrix N2 = fl.col(1)*fl.col(1).transpose(); - Matrix N12 = 0.5*(fl.col(0)*fl.col(1).transpose() + fl.col(1)*fl.col(0).transpose()); - - - dInv6 = -Inv[5]/3*Ci + J2d*N12; - dInv7 = J4d*(N12*C + C*N12) - Inv[6]/3*Ci; - dInv8 = -Inv[7]/3*Ci + J2d*N2; - dInv9 = J4d*(N2*C + C*N2) - Inv[8]/3*Ci; - - ddInv6 = -1.0/3.0*(dyadic_product(dInv6,Ci) + J2d*dyadic_product(Ci,N12) + Inv[5]*dCidC); - Matrix sum12 = (N12*C + C*N12); - ddInv7 = -1.0/3.0*(dyadic_product(dInv7,Ci) + Inv[6]*dCidC + 2*J4d*dyadic_product(Ci,sum12)) - + J4d*(2*symmetric_dyadic_product(N12,Idm) - dyadic_product(N12,Idm) - + 2*symmetric_dyadic_product(Idm,N12) - dyadic_product(Idm,N12)); - ddInv8 = -1.0/3.0*(dyadic_product(dInv8,Ci) + J2d*dyadic_product(Ci,N2) + Inv[7]*dCidC); - Matrix sum2 = (N2*C + C*N2); - ddInv9 = -1.0/3.0*(dyadic_product(dInv9,Ci) + Inv[8]*dCidC + 2*J4d*dyadic_product(Ci,sum2)) - + J4d*(2*symmetric_dyadic_product(N2,Idm) - dyadic_product(N2,Idm) - + 2*symmetric_dyadic_product(Idm,N2) - dyadic_product(Idm,N2)); - } - - dInv = {dInv1, dInv2, dInv3, dInv4, dInv5, dInv6, dInv7, dInv8, dInv9}; - ddInv = {ddInv1, ddInv2, ddInv3, ddInv4, ddInv5, ddInv6, ddInv7, ddInv8, ddInv9}; + const Matrix &C, const Matrix &fl, int nfd, double J2d, + double J4d, const Matrix &Ci, const Matrix &Idm, Matrix &N1, + double &psi, double (&Inv)[9], std::array, 9> &dInv, + std::array, 9> &ddInv) const { + + // double Inv[9] = {0}; + Matrix C2 = C * C; + + Inv[0] = J2d * C.trace(); + Inv[1] = 0.50 * (Inv[0] * Inv[0] - J4d * (C * C).trace()); + Inv[2] = C.determinant(); + Inv[3] = J2d * (fl.col(0).dot(C * fl.col(0))); + Inv[4] = J4d * (fl.col(0).dot(C2 * fl.col(0))); + + Matrix dInv1 = -Inv[0] / 3 * Ci + J2d * Idm; + Matrix dInv2 = (C2.trace() / 3) * Ci + Inv[0] * dInv1 + J4d * C; + Matrix dInv3 = Inv[2] * Ci; + N1 = fl.col(0) * fl.col(0).transpose(); + Matrix dInv4 = -Inv[3] / 3 * Ci + J2d * N1; + Matrix dInv5 = J4d * (N1 * C + C * N1) - Inv[4] / 3 * Ci; + + Matrix dInv6, dInv7, dInv8, dInv9; + Tensor ddInv6, ddInv7, ddInv8, ddInv9; + // initialize to 0 + dInv6.setZero(); + dInv7.setZero(); + dInv8.setZero(); + dInv9.setZero(); + ddInv6.setZero(); + ddInv7.setZero(); + ddInv8.setZero(); + ddInv9.setZero(); + + Tensor dCidC = -symmetric_dyadic_product(Ci, Ci); + Matrix dJ4ddC = -2.0 / 3.0 * J4d * Ci; + + Tensor ddInv1 = + (-1.0 / 3.0) * (dyadic_product(dInv1, Ci) + Inv[0] * dCidC + + J2d * dyadic_product(Ci, Idm)); + Tensor ddInv2 = + dyadic_product(dInv1, dInv1) + Inv[0] * ddInv1 + + (1.0 / 3.0) * C2.trace() * dCidC + + (1.0 / 3.0) * + dyadic_product((C2.trace() * dJ4ddC + 2 * J4d * C), Ci) + + dyadic_product(dJ4ddC, C) - J4d * fourth_order_identity(); + Tensor ddInv3 = dyadic_product(dInv3, Ci) + Inv[2] * dCidC; + Tensor ddInv4 = + (-1.0 / 3.0) * (dyadic_product(dInv4, Ci) + + J2d * dyadic_product(Ci, N1) + Inv[3] * dCidC); + Matrix sum1 = (N1 * C + C * N1); + Tensor ddInv5 = + (-1.0 / 3.0) * (dyadic_product(dInv5, Ci) + Inv[4] * dCidC + + 2 * J4d * dyadic_product(Ci, sum1)) + + J4d * (2 * symmetric_dyadic_product(N1, Idm) - dyadic_product(N1, Idm) + + 2 * symmetric_dyadic_product(Idm, N1) - dyadic_product(Idm, N1)); + + if (nfd == 2) { + Inv[5] = J2d * (fl.col(0).dot(C * fl.col(1))); + Inv[6] = J4d * (fl.col(0).dot(C2 * fl.col(1))); + Inv[7] = J2d * (fl.col(1).dot(C * fl.col(1))); + Inv[8] = J4d * (fl.col(1).dot(C2 * fl.col(1))); + + Matrix N2 = fl.col(1) * fl.col(1).transpose(); + Matrix N12 = 0.5 * (fl.col(0) * fl.col(1).transpose() + + fl.col(1) * fl.col(0).transpose()); + + dInv6 = -Inv[5] / 3 * Ci + J2d * N12; + dInv7 = J4d * (N12 * C + C * N12) - Inv[6] / 3 * Ci; + dInv8 = -Inv[7] / 3 * Ci + J2d * N2; + dInv9 = J4d * (N2 * C + C * N2) - Inv[8] / 3 * Ci; + + ddInv6 = -1.0 / 3.0 * + (dyadic_product(dInv6, Ci) + J2d * dyadic_product(Ci, N12) + + Inv[5] * dCidC); + Matrix sum12 = (N12 * C + C * N12); + ddInv7 = + -1.0 / 3.0 * + (dyadic_product(dInv7, Ci) + Inv[6] * dCidC + + 2 * J4d * dyadic_product(Ci, sum12)) + + J4d * + (2 * symmetric_dyadic_product(N12, Idm) - dyadic_product(N12, Idm) + + 2 * symmetric_dyadic_product(Idm, N12) - dyadic_product(Idm, N12)); + ddInv8 = -1.0 / 3.0 * + (dyadic_product(dInv8, Ci) + J2d * dyadic_product(Ci, N2) + + Inv[7] * dCidC); + Matrix sum2 = (N2 * C + C * N2); + ddInv9 = + -1.0 / 3.0 * + (dyadic_product(dInv9, Ci) + Inv[8] * dCidC + + 2 * J4d * dyadic_product(Ci, sum2)) + + J4d * (2 * symmetric_dyadic_product(N2, Idm) - dyadic_product(N2, Idm) + + 2 * symmetric_dyadic_product(Idm, N2) - dyadic_product(Idm, N2)); + } + + dInv = {dInv1, dInv2, dInv3, dInv4, dInv5, dInv6, dInv7, dInv8, dInv9}; + ddInv = {ddInv1, ddInv2, ddInv3, ddInv4, ddInv5, + ddInv6, ddInv7, ddInv8, ddInv9}; } - // Template instantiation template void ArtificialNeuralNetMaterial::computeInvariantsAndDerivatives<2>( -const Matrix<2>& C, const Matrix<2>& fl, int nfd, double J2d, double J4d, const Matrix<2>& Ci, -const Matrix<2>& Idm, const double Tfa, Matrix<2>& N1, double& psi, double (&Inv)[9], std::array,9>& dInv, -std::array,9>& ddInv) const; + const Matrix<2> &C, const Matrix<2> &fl, int nfd, double J2d, double J4d, + const Matrix<2> &Ci, const Matrix<2> &Idm, Matrix<2> &N1, double &psi, + double (&Inv)[9], std::array, 9> &dInv, + std::array, 9> &ddInv) const; template void ArtificialNeuralNetMaterial::computeInvariantsAndDerivatives<3>( -const Matrix<3>& C, const Matrix<3>& fl, int nfd, double J2d, double J4d, const Matrix<3>& Ci, -const Matrix<3>& Idm, const double Tfa, Matrix<3>& N1, double& psi, double (&Inv)[9], std::array,9>& dInv, -std::array,9>& ddInv) const; + const Matrix<3> &C, const Matrix<3> &fl, int nfd, double J2d, double J4d, + const Matrix<3> &Ci, const Matrix<3> &Idm, Matrix<3> &N1, double &psi, + double (&Inv)[9], std::array, 9> &dInv, + std::array, 9> &ddInv) const; diff --git a/Code/Source/solver/ArtificialNeuralNetMaterial.h b/Code/Source/solver/ArtificialNeuralNetMaterial.h index 8500ed344..0b3d81be8 100644 --- a/Code/Source/solver/ArtificialNeuralNetMaterial.h +++ b/Code/Source/solver/ArtificialNeuralNetMaterial.h @@ -57,12 +57,13 @@ class ArtificialNeuralNetMaterial void evaluate(const double aInv[9], double &psi, double (&dpsi)[9], double (&ddpsi)[9]) const; // Helper for compute_pk2cc - template + template void computeInvariantsAndDerivatives( - const Matrix& C, const Matrix& fl, int nfd, double J2d, double J4d, const Matrix& Ci, - const Matrix& Idm, const double Tfa, Matrix& N1, double& psi, double (&Inv)[9], std::array,9>& dInv, - std::array,9>& ddInv) const; - + const Matrix &C, const Matrix &fl, int nfd, double J2d, + double J4d, const Matrix &Ci, const Matrix &Idm, + Matrix &N1, double &psi, double (&Inv)[9], + std::array, 9> &dInv, + std::array, 9> &ddInv) const; }; #endif // ArtificialNeuralNet_model_H \ No newline at end of file diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index 950b94896..0b78b0eef 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -186,37 +186,6 @@ class cepModelType std::shared_ptr ionic_model; }; -/// @brief Cardiac electromechanics model type -class cemModelType -{ - public: - /// @brief Whether electrophysiology and mechanics are coupled - bool cpld = false; - //bool cpld = .FALSE. - - /// @brief Whether active strain formulation is employed - bool aStrain = false; - //bool aStrain = .FALSE. - - /// @brief Activation along fibers. - /// - /// Corresponds to active tension along fibers if using active stress, and - /// to fiber stretch if using active strain. - Vector Ya_f; - - /// @brief Activation along sheets. - /// - /// Only used if using active stress, in which case it represents the active - /// tension along sheets. - Vector Ya_s; - - /// @brief Activation along sheet normals. - /// - /// Only used if using active stress, in which case it represents the active - /// tension along sheet normals. - Vector Ya_n; -}; - class CepMod { public: @@ -233,9 +202,6 @@ class CepMod /// @brief Calcium vector at all nodes. Vector calcium; - /// @brief Cardiac electromechanics type - cemModelType cem; - /// @brief ECG leads ecgLeadsType ecgleads; }; diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index dfdedc854..8affac66d 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -110,6 +110,19 @@ bool Integrator::step(bool save_results) { // Compute body forces set_body_forces(); + // Implicit state coupling of the active stress: re-advance the state of the + // active stress model from the displacement of the current nonlinear + // iterate. + if (supports_active_stress(eq.phys) && + has_implicit_active_stress_state_coupling()) { + Vector fiber_stretch; + Vector fiber_stretch_rate; + compute_fiber_stretch(fiber_stretch, fiber_stretch_rate); + + update_active_stress(eq, fiber_stretch, fiber_stretch_rate, + /* within_nonlinear_iterations = */ true); + } + // Assemble equations assemble_equations(); @@ -371,6 +384,140 @@ void Integrator::update_residual_arrays(eqType& eq) { } } +//------------------------ +// compute_fiber_stretch +//------------------------ +void Integrator::compute_fiber_stretch(Vector& fiber_stretch, Vector& fiber_stretch_rate) { + using namespace consts; + + auto& com_mod = simulation_->com_mod; + const auto& Dn = solutions_.current.get_displacement(); + + // Determine if we need to compute fiber stretch and stretch rate, by going + // through all domains of all equations until we find one for which active + // stress is enabled and the active stress model needs the stretch or stretch rate. + // + // have_active_stress is tracked separately from need_fiber_stretch because + // update() indexes both vectors for every node whether or not the + // model reads the values, so they must be allocated either way. + bool have_active_stress = false; + bool need_fiber_stretch = false; + bool need_fiber_stretch_rate = false; + int fiber_stretch_eq_index = -1; + for (int iEq = 0; iEq < com_mod.nEq; ++iEq) { + const auto &eq = com_mod.eq[iEq]; + + if (supports_active_stress(eq.phys)) { + fiber_stretch_eq_index = iEq; + + for (const auto &dmn : eq.dmn) { + if (dmn.active_stress != nullptr) { + have_active_stress = true; + need_fiber_stretch |= dmn.active_stress->needs_fiber_stretch(); + need_fiber_stretch_rate |= + dmn.active_stress->needs_fiber_stretch_rate(); + } + } + } + + else if (eq.phys == Equation_CEP) { + need_fiber_stretch = true; + } + } + + // If we need to compute fiber stretch, we iterate through all meshes, compute + // the stretch for each mesh, and then copy the mesh-local resulting vector + // into the global vector. + if (have_active_stress || need_fiber_stretch) { + fiber_stretch.resize(com_mod.tnNo); + + if (need_fiber_stretch && fiber_stretch_eq_index >= 0) { + for (const auto &mesh : com_mod.msh) { + Vector tmp(mesh.nNo); + + post::fib_stretch(com_mod, fiber_stretch_eq_index, mesh, Dn, tmp); + for (int a = 0; a < mesh.nNo; ++a) + fiber_stretch[mesh.gN[a]] = tmp[a]; + } + } else { + // No domain solves for the displacement, or no model reads the stretch: + // set the fiber stretch to 1, corresponding to no stretch. + fiber_stretch = 1.0; + } + } + + // Same for fiber stretch rate. + if (have_active_stress) { + fiber_stretch_rate.resize(com_mod.tnNo); + + if (need_fiber_stretch_rate && fiber_stretch_eq_index >= 0) { + for (const auto &mesh : com_mod.msh) { + Vector tmp(mesh.nNo); + + post::fib_stretch_rate(com_mod, fiber_stretch_eq_index, mesh, + solutions_, tmp); + for (int a = 0; a < mesh.nNo; ++a) + fiber_stretch_rate[mesh.gN[a]] = tmp[a]; + } + } else { + // No domain solves for displacement, or no model reads stretch rate: + // Set the fiber stretch rate to 0, corresponding to no movement. + fiber_stretch_rate = 0.0; + } + } +} + +//------------------------ +// has_implicit_active_stress_state_coupling +//------------------------ +bool Integrator::has_implicit_active_stress_state_coupling() const { + const auto& com_mod = simulation_->com_mod; + + for (const auto &eq : com_mod.eq) { + if (!supports_active_stress(eq.phys)) + continue; + + for (const auto &dmn : eq.dmn) { + if (dmn.active_stress != nullptr && + dmn.active_stress->implicit_state_coupling()) + return true; + } + } + + return false; +} + +void Integrator::time_advance_active_stress(eqType &eq) { + for (auto &dmn : eq.dmn) { + if (dmn.active_stress != nullptr) { + dmn.active_stress->time_advance(); + } + } +} + +//------------------------ +// update_active_stress +//------------------------ +void Integrator::update_active_stress( + eqType &eq, const Vector &fiber_stretch, + const Vector &fiber_stretch_rate, + const bool within_nonlinear_iterations) { + auto& com_mod = simulation_->com_mod; + auto& cep_mod = simulation_->get_cep_mod(); + + for (auto &dmn : eq.dmn) { + if (dmn.active_stress == nullptr) + continue; + + // Models with explicit coupling keep the state computed in the predictor. + if (within_nonlinear_iterations && + !dmn.active_stress->implicit_state_coupling()) + continue; + + dmn.active_stress->update(com_mod.time, com_mod.dt, cep_mod.calcium, + fiber_stretch, fiber_stretch_rate); + } +} // The code here replicates the Fortran code in PIC.f. // @@ -397,7 +544,6 @@ void Integrator::predictor() using namespace consts; auto& com_mod = simulation_->com_mod; - auto& cep_mod = simulation_->cep_mod; #define n_debug_picp #ifdef debug_picp @@ -464,79 +610,7 @@ void Integrator::predictor() Vector fiber_stretch; Vector fiber_stretch_rate; - - // Determine if we need to compute fiber stretch and stretch rate, by going - // through all domains of all equations until we find one for which active - // stress is enabled and the active stress model needs the stretch or stretch rate. - // - // have_active_stress is tracked separately from need_fiber_stretch because - // advance_time_step() indexes both vectors for every node whether or not the - // model reads the values, so they must be allocated either way. - bool have_active_stress = false; - bool need_fiber_stretch = false; - bool need_fiber_stretch_rate = false; - int fiber_stretch_eq_index = -1; - for (int iEq = 0; iEq < com_mod.nEq; ++iEq) { - const auto &eq = com_mod.eq[iEq]; - - if (supports_active_stress(eq.phys)) { - fiber_stretch_eq_index = iEq; - - for (const auto &dmn : eq.dmn) { - if (dmn.active_stress != nullptr) { - have_active_stress = true; - need_fiber_stretch |= dmn.active_stress->needs_fiber_stretch(); - need_fiber_stretch_rate |= - dmn.active_stress->needs_fiber_stretch_rate(); - } - } - } - - else if (eq.phys == Equation_CEP) { - need_fiber_stretch = true; - } - } - - // If we need to compute fiber stretch, we iterate through all meshes, compute - // the stretch for each mesh, and then copy the mesh-local resulting vector - // into the global vector. - if (have_active_stress || need_fiber_stretch) { - fiber_stretch.resize(com_mod.tnNo); - - if (need_fiber_stretch && fiber_stretch_eq_index >= 0) { - for (const auto &mesh : com_mod.msh) { - Vector tmp(mesh.nNo); - - post::fib_stretch(com_mod, fiber_stretch_eq_index, mesh, Dn, tmp); - for (int a = 0; a < mesh.nNo; ++a) - fiber_stretch[mesh.gN[a]] = tmp[a]; - } - } else { - // No domain solves for the displacement, or no model reads the stretch: - // set the fiber stretch to 1, corresponding to no stretch. - fiber_stretch = 1.0; - } - } - - // Same for fiber stretch rate. - if (have_active_stress) { - fiber_stretch_rate.resize(com_mod.tnNo); - - if (need_fiber_stretch_rate && fiber_stretch_eq_index >= 0) { - for (const auto &mesh : com_mod.msh) { - Vector tmp(mesh.nNo); - - post::fib_stretch_rate(com_mod, fiber_stretch_eq_index, mesh, - solutions_, tmp); - for (int a = 0; a < mesh.nNo; ++a) - fiber_stretch_rate[mesh.gN[a]] = tmp[a]; - } - } else { - // No domain solves for displacement, or no model reads stretch rate: - // Set the fiber stretch rate to 0, corresponding to no movement. - fiber_stretch_rate = 0.0; - } - } + compute_fiber_stretch(fiber_stretch, fiber_stretch_rate); for (int iEq = 0; iEq < com_mod.nEq; iEq++) { auto& eq = com_mod.eq[iEq]; @@ -571,54 +645,9 @@ void Integrator::predictor() // active stress if (supports_active_stress(eq.phys)) { - for (auto &dmn : eq.dmn) { - if (dmn.active_stress != nullptr) { - dmn.active_stress->advance_time_step(com_mod.time, com_mod.dt, - cep_mod.calcium, fiber_stretch, - fiber_stretch_rate); - } - } - - // Fill in the active tension vector. - // We go through all mesh nodes, find the domain they are associated with, - // and get the active stress from that domain. If a point is associated to - // multiple domains (which happens for points on domain interfaces), we - // average the active stresses from the domains. - for (int Ac = 0; Ac < com_mod.tnNo; Ac++) { - double Ta_f = 0.0; - double Ta_s = 0.0; - double Ta_n = 0.0; - unsigned int n_domains = 0; - - for (auto &dmn : eq.dmn) { - // Domains whose equations do not allow for active stress (e.g. fluid - // domains) do not contribute to the average, but domains that do - // allow for active stress (e.g. struct) for which active stress is - // not enabled contribute a zero value to the average. - if (!supports_active_stress(dmn.phys)) - continue; - - // Only domains that node Ac actually belongs to contribute to its - // average. Note that if there is only one domain dmnId may not be - // populated, so we only check domain membership if eq.nDmn > 1. - if (eq.nDmn > 1 && !utils::btest(com_mod.dmnId(Ac), dmn.Id)) - continue; - - if (dmn.active_stress != nullptr) { - Ta_f += dmn.active_stress->get_tension_fibers(Ac); - Ta_s += dmn.active_stress->get_tension_sheets(Ac); - Ta_n += dmn.active_stress->get_tension_sheet_normals(Ac); - } - - n_domains++; - } - - if (n_domains > 0) { - cep_mod.cem.Ya_f[Ac] = Ta_f / n_domains; - cep_mod.cem.Ya_s[Ac] = Ta_s / n_domains; - cep_mod.cem.Ya_n[Ac] = Ta_n / n_domains; - } - } + time_advance_active_stress(eq); + update_active_stress(eq, fiber_stretch, fiber_stretch_rate, + /* within_nonlinear_iterations = */ false); } // eqn 86 of Bazilevs 2007 diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h index 660a567bb..b0ea511d6 100644 --- a/Code/Source/solver/Integrator.h +++ b/Code/Source/solver/Integrator.h @@ -167,6 +167,51 @@ class Integrator { */ void update_residual_arrays(eqType& eq); + /** + * @brief Compute the fiber stretch and fiber stretch rate at every node. + * + * The stretch is computed from the current displacement, so that calling this + * function within the nonlinear iterations yields the stretch of the current + * displacement iterate. Both vectors are left empty if no equation needs + * them, and are filled with the neutral values (1 for the stretch, 0 for the + * stretch rate) if they are needed but no equation solves for the + * displacement. + * + * @param[out] fiber_stretch Fiber stretch at every node. + * @param[out] fiber_stretch_rate Fiber stretch rate at every node. + */ + void compute_fiber_stretch(Vector& fiber_stretch, + Vector& fiber_stretch_rate); + + /** + * @brief Whether any domain uses an active stress model with implicit state + * coupling, i.e. one whose state is updated within the nonlinear iterations. + * + * @return True if at least one such domain exists, false otherwise. + */ + bool has_implicit_active_stress_state_coupling() const; + + /** + * @brief Advance active stress models of an equation to the next time step. + */ + void time_advance_active_stress(eqType &eq); + + /** + * @brief Update the active stress models of an equation. + * + * @param[in,out] eq Equation whose domains carry the active stress models. + * @param[in] fiber_stretch Fiber stretch at every node. + * @param[in] fiber_stretch_rate Fiber stretch rate at every node. + * @param[in] within_nonlinear_iterations True when called within the + * nonlinear iterations, in which case only the models with implicit + * coupling are advanced again, from the state stored at the beginning of + * the time step. False when called once per time step from the predictor, + * in which case all models are advanced from that state. + */ + void update_active_stress(eqType &eq, const Vector &fiber_stretch, + const Vector &fiber_stretch_rate, + const bool within_nonlinear_iterations); + /** * @brief Initiator function for generalized-alpha method (initiator) * diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 9ff3eca0f..f6b56edb7 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1908,9 +1908,13 @@ void ActiveStressModelParameters::set_values( const std::string ActiveStressParameters::xml_element_name = "Active_stress"; ActiveStressParameters::ActiveStressParameters() { + constexpr bool required = true; + model_name = Parameter("Model", "", true); - set_parameter("Model", "", /* required = */ true, model_name); + set_parameter("Model", "", required, model_name); + set_parameter("Implicit_state_coupling", false, !required, + implicit_state_coupling); ActiveStressFactory::visit( [this](const std::string &name, const ActiveStress &model) { @@ -1984,6 +1988,10 @@ double ActiveStressParameters::get_eta_n() const { return directional_distribution.sheet_normal_direction.value(); } +bool ActiveStressParameters::get_implicit_state_coupling() const { + return implicit_state_coupling.value(); +} + const ActiveStressModelParameters & ActiveStressParameters::get_parameters(const std::string &model_name) const { if (active_stress_models.count(model_name) == 0) { diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index f97336688..547eef62e 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1558,6 +1558,10 @@ class ActiveStressParameters : public ParameterLists { /// Get the active tension coefficient along sheet normals. double get_eta_n() const; + /// Get whether the state of the active stress model is updated within the + /// nonlinear iterations of the mechanics problem. + bool get_implicit_state_coupling() const; + /// Get the parameters for a given active stress model. const ActiveStressModelParameters & get_parameters(const std::string &model_name) const; @@ -1569,6 +1573,10 @@ class ActiveStressParameters : public ParameterLists { /// Parameter for the model name. Parameter model_name; + /// Parameter selecting whether the state of the active stress model is + /// updated within the nonlinear iterations of the mechanics problem. + Parameter implicit_state_coupling; + /// Parameters for the directional distribution of active tension. DirectionalDistributionParameters directional_distribution; diff --git a/Code/Source/solver/cep.cpp b/Code/Source/solver/cep.cpp index 86416e86d..2e15ef101 100644 --- a/Code/Source/solver/cep.cpp +++ b/Code/Source/solver/cep.cpp @@ -105,7 +105,6 @@ void cep_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, con auto& cDmn = com_mod.cDmn; auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - const auto& cem = cep_mod.cem; Vector Dani(nFn), Vx(2), Ls(nFn), DVx(2); Array F(2,2), C(2,2), fl(2,nFn), D(2,2), DNx(2,eNoN); @@ -137,65 +136,17 @@ void cep_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, con Ls = 1.0; int i; - if (cem.cpld) { - for (int a = 0; a < com_mod.nEq; a++) { - if (com_mod.eq[a].phys == EquationType::phys_struct || - com_mod.eq[a].phys == EquationType::phys_ustruct) { - i = com_mod.eq[a].s; - break; - } - } - - // Compute deformation gradient tensor - // - F(0,0) = 1.0; - F(1,1) = 1.0; - - for (int a = 0; a < eNoN; a++) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(i+1,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(i+1,a); - } - - // Jacobian - double Jac = mat_fun::mat_det(F, 2); - - // Compute Cauchy-Green tensor and its inverse - C = mat_mul(transpose(F), F); - C = mat_inv(C, 2); + D = 0.0; + D(0, 0) = Diso; + D(1, 1) = Diso; + fl = fN; - // Compute fiber stretch - for (int i = 0; i < nFn; i++) { - Ls(i) = sqrt(fN.rcol(i) * mat_mul(C, fN.rcol(i))); - for (int j = 0; j < 2; j++) { - fl(j,i) = fN(j,i) / Ls(i); - } - } - - if (Ls(0) <= 1.0) { - Ls(0) = 1.0; - } - - // Diffusion tensor - spatial isotropy - // - Diso = Diso * Jac; - Dani = Dani * Jac; - D = Diso * C; - - } else { - D = 0.0; - D(0,0) = Diso; - D(1,1) = Diso; - fl = fN; - } - - for (int i = 0 ; i < nFn; i++) { - D(0,0) = D(0,0) + Dani(i)*fl(0,i)*fl(0,i); - D(0,1) = D(0,1) + Dani(i)*fl(0,i)*fl(1,i); + for (int i = 0; i < nFn; i++) { + D(0, 0) = D(0, 0) + Dani(i) * fl(0, i) * fl(0, i); + D(0, 1) = D(0, 1) + Dani(i) * fl(0, i) * fl(1, i); - D(1,0) = D(1,0) + Dani(i)*fl(1,i)*fl(0,i); - D(1,1) = D(1,1) + Dani(i)*fl(1,i)*fl(1,i); + D(1, 0) = D(1, 0) + Dani(i) * fl(1, i) * fl(0, i); + D(1, 1) = D(1, 1) + Dani(i) * fl(1, i) * fl(1, i); } i = eq.s; @@ -252,7 +203,6 @@ void cep_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, con auto& cDmn = com_mod.cDmn; auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - const auto& cem = cep_mod.cem; Vector Dani(nFn), Vx(3), Ls(nFn), DVx(3); Array F(3,3), C(3,3), fl(3,nFn), D(3,3), DNx(3,eNoN); @@ -283,65 +233,10 @@ void cep_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, con Ls = 1.0; int i = 0; - if (cem.cpld) { - // Get the displacement degrees of freedom - for (int a = 0; a < com_mod.nEq; a++) { - if (com_mod.eq[a].phys == EquationType::phys_struct || - com_mod.eq[a].phys == EquationType::phys_ustruct) { - i = com_mod.eq[a].s; - break; - } - } - - // Compute deformation gradient tensor - // - F(0,0) = 1.0; - F(1,1) = 1.0; - F(2,2) = 1.0; - - for (int a = 0; a < eNoN; a++) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(i+1,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(i+1,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(i+1,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(i+2,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(i+2,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(i+2,a); - } - - // Jacobian - double Jac = mat_fun::mat_det(F, 3); - - // Compute Cauchy-Green tensor and its inverse - C = mat_mul(transpose(F), F); - C = mat_inv(C, 3); - - // Compute fiber stretch - for (int i = 0; i < nFn; i++) { - Ls(i) = sqrt(fN.rcol(i) * mat_mul(C, fN.rcol(i))); - for (int j = 0; j < 3; j++) { - fl(j,i) = fN(j,i) / Ls(i); - } - } - - if (Ls(0) <= 1.0) { - Ls(0) = 1.0; - } - - // Diffusion tensor - spatial isotropy - // - Diso = Diso * Jac; - Dani = Dani * Jac; - D = Diso * C; - - } else { - D(0,0) = Diso; - D(1,1) = Diso; - D(2,2) = Diso; - fl= fN; - } + D(0, 0) = Diso; + D(1, 1) = Diso; + D(2, 2) = Diso; + fl = fN; // Compute anisotropic components of diffusion tensor // diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 81a9a7493..f01b86e68 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -140,14 +140,12 @@ void cep_integ(Simulation *simulation, const int iEq, const int iDof, double time = com_mod.time; auto &cep_mod = simulation->cep_mod; - auto &cem = cep_mod.cem; auto &eq = com_mod.eq[iEq]; auto &Xion = cep_mod.Xion; int nXion = cep_mod.nXion; #ifdef debug_cep_integ - dmsg << "cem.cpld: " << cem.cpld; dmsg << "time: " << time; #endif @@ -223,10 +221,6 @@ void cep_integ(Simulation *simulation, const int iEq, const int iDof, all_fun::commu(com_mod, sF); all_fun::commu(com_mod, cep_mod.calcium); - if (cem.cpld) { - all_fun::commu(com_mod, sY); - } - for (int Ac = 0; Ac < tnNo; Ac++) { if (!utils::is_zero(sA(Ac))) { Xion.set_col(Ac, sF.col(Ac) / sA(Ac)); diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index acc49e3e3..88ce08630 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1599,14 +1599,6 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: } } - // Distribute cardiac electromechanics parameters - // - cm.bcast(cm_mod, &cep_mod.cem.cpld); - - if (cep_mod.cem.cpld) { - cm.bcast(cm_mod, &cep_mod.cem.aStrain); - } - if (com_mod.ibFlag) { if (cm.slv(cm_mod)) { lEq.dmnIB.resize(lEq.nDmnIB); diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index bfd54c254..e867f9b67 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -47,7 +47,6 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So } // l = 3, if nsd==2 ; else 6; - auto& cem = cep_mod.cem; const int l = com_mod.nsymd; const int nsd = com_mod.nsd; const int tDof = com_mod.tDof; @@ -69,7 +68,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So Array3 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN); Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), lR(dof,eNoN); - Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN); + Vector pSl(nsymd); + ActiveStress::Evaluator active_stress_evaluator; std::array fs_1; fs::get_thood_fs(com_mod, fs_1, lM, vmsStab, 1); @@ -98,9 +98,6 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So // Create local copies fN = 0.0; pS0l = 0.0; - ya_l_f = 0.0; - ya_l_s = 0.0; - ya_l_n = 0.0; if (lM.fN.size() != 0) { for (int iFn = 0; iFn < nFn; iFn++) { @@ -128,11 +125,12 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So pS0l.set_col(a, pS0.col(Ac)); } - if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l_f(a) = cep_mod.cem.Ya_f[Ac]; - ya_l_s(a) = cep_mod.cem.Ya_s[Ac]; - ya_l_n(a) = cep_mod.cem.Ya_n[Ac]; - } + } + + if (eq.dmn[cDmn].active_stress != nullptr) { + active_stress_evaluator.update(*eq.dmn[cDmn].active_stress, ptr); + } else { + active_stress_evaluator.clear(); } // For FSI, fluid domain should be in the current configuration @@ -223,8 +221,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So case Equation_struct: { auto N0 = fs_1[0].N.col(g); struct_ns::struct_3d(com_mod, cep_mod, fs_1[0].eNoN, nFn, w, N0, - Nwx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l_f, - ya_l_s, ya_l_n, lR, lK, recompute_visc); + Nwx, al, yl, dl, bfl, fN, pS0l, pSl, + active_stress_evaluator, lR, lK, recompute_visc); } break; case Equation_lElas: throw std::runtime_error("[construct_fsi] LELAS3D not implemented"); @@ -236,7 +234,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So auto N1 = fs_1[1].N.col(g); ustruct::ustruct_3d_m(com_mod, cep_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, - yl, dl, bfl, fN, ya_l_f, ya_l_s, ya_l_n, lR, + yl, dl, bfl, fN, active_stress_evaluator, lR, lK, lKd, recompute_visc); break; } @@ -259,8 +257,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So case Equation_struct: { auto N0 = fs_1[0].N.col(g); struct_ns::struct_2d(com_mod, cep_mod, fs_1[0].eNoN, nFn, w, N0, - Nwx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l_f, - ya_l_s, ya_l_n, lR, lK, recompute_visc); + Nwx, al, yl, dl, bfl, fN, pS0l, pSl, + active_stress_evaluator, lR, lK, recompute_visc); } break; case Equation_ustruct: diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 89dad332e..99479aa3a 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -69,7 +69,6 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< auto const recLn = com_mod.recLn; auto& cm_mod = simulation->cm_mod; auto& cep_mod = simulation->cep_mod; - auto& cem = cep_mod.cem; bool ibFlag = com_mod.ibFlag; bool dFlag = com_mod.dFlag; @@ -128,13 +127,10 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< } else if (cepEq) { bin_file.read((char*)Ad.data(), Ad.msize()); bin_file.read((char*)Xion.data(), Xion.msize()); - bin_file.read((char*)cem.Ya_f.data(), cem.Ya_f.msize()); - bin_file.read((char*)cem.Ya_s.data(), cem.Ya_s.msize()); - bin_file.read((char*)cem.Ya_n.data(), cem.Ya_n.msize()); } else if (risFlag) { bin_file.read((char*)Ad.data(), Ad.msize()); - init_ris_data(com_mod, bin_file); + init_ris_data(com_mod, bin_file); } else if (urisFlag) { bin_file.read((char*)Ad.data(), Ad.msize()); @@ -151,12 +147,9 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< } else if (cepEq) { bin_file.read((char*)Xion.data(), Xion.msize()); - bin_file.read((char*)cem.Ya_f.data(), cem.Ya_f.msize()); - bin_file.read((char*)cem.Ya_s.data(), cem.Ya_s.msize()); - bin_file.read((char*)cem.Ya_n.data(), cem.Ya_n.msize()); } else if (risFlag) { - init_ris_data(com_mod, bin_file); + init_ris_data(com_mod, bin_file); } else if (urisFlag) { init_uris_data(com_mod, bin_file); @@ -543,7 +536,6 @@ void initialize(Simulation* simulation, Vector& timeP) if (com_mod.sstEq) i = i + nsd; if (cep_mod.cepEq) { i = i + cep_mod.nXion; - if (cep_mod.cem.cpld) i = i + 1; } if (com_mod.risFlag) { i = i + com_mod.ris.nbrRIS; @@ -707,9 +699,6 @@ void initialize(Simulation* simulation, Vector& timeP) // initialized. { cep_mod.calcium.resize(tnNo); - cep_mod.cem.Ya_f.resize(tnNo); - cep_mod.cem.Ya_s.resize(tnNo); - cep_mod.cem.Ya_n.resize(tnNo); } // Setup the initial conditions for the active stress models. diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index 2afdcff2f..cff2462bd 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -95,56 +95,42 @@ void cc_to_voigt(const int nsd, const Tensor4& CC, Array& Dm) } } +/** + * @brief Write a 4th order elasticity tensor in Voigt notation. + * + * Every entry is read from the tensor independently, so that an elasticity + * tensor that is not major symmetric is written faithfully. Active stress + * produces such a tensor, since it does not derive from a strain energy. + */ template void cc_to_voigt_eigen(const Tensor& CC, Matrix<3*(nsd-1)>& Dm) { - if (nsd == 3) { - Dm(0,0) = CC(0,0,0,0); - Dm(0,1) = CC(0,0,1,1); - Dm(0,2) = CC(0,0,2,2); - Dm(0,3) = CC(0,0,0,1); - Dm(0,4) = CC(0,0,1,2); - Dm(0,5) = CC(0,0,2,0); - - Dm(1,1) = CC(1,1,1,1); - Dm(1,2) = CC(1,1,2,2); - Dm(1,3) = CC(1,1,0,1); - Dm(1,4) = CC(1,1,1,2); - Dm(1,5) = CC(1,1,2,0); - - Dm(2,2) = CC(2,2,2,2); - Dm(2,3) = CC(2,2,0,1); - Dm(2,4) = CC(2,2,1,2); - Dm(2,5) = CC(2,2,2,0); - - Dm(3,3) = CC(0,1,0,1); - Dm(3,4) = CC(0,1,1,2); - Dm(3,5) = CC(0,1,2,0); - - Dm(4,4) = CC(1,2,1,2); - Dm(4,5) = CC(1,2,2,0); - - Dm(5,5) = CC(2,0,2,0); - - for (int i = 1; i < 6; i++) { - for (int j = 0; j <= i-1; j++) { - Dm(i,j) = Dm(j,i); - } + // Index pairs of the tensor corresponding to each index in Voigt notation. + constexpr int n_voigt = 3 * (nsd - 1); + + constexpr std::array voigt_row = []() { + if constexpr (nsd == 3) + return std::array{0, 1, 2, 0, 1, 2}; + else + return std::array{0, 1, 0}; + }(); + + constexpr std::array voigt_col = []() { + if constexpr (nsd == 3) + return std::array{0, 1, 2, 1, 2, 0}; + else + return std::array{0, 1, 1}; + }(); + + for (int i = 0; i < n_voigt; i++) { + for (int j = 0; j < n_voigt; j++) { + const int i_row = voigt_row[i]; + const int i_col = voigt_col[i]; + const int j_row = voigt_row[j]; + const int j_col = voigt_col[j]; + + Dm(i,j) = CC(i_row, i_col, j_row, j_col); } - - } else if (nsd == 2) { - Dm(0,0) = CC(0,0,0,0); - Dm(0,1) = CC(0,0,1,1); - Dm(0,2) = CC(0,0,0,1); - - Dm(1,1) = CC(1,1,1,1); - Dm(1,2) = CC(1,1,0,1); - - Dm(2,2) = CC(0,1,0,1); - - Dm(1,0) = Dm(0,1); - Dm(2,0) = Dm(0,2); - Dm(2,1) = Dm(1,2); } } @@ -285,7 +271,7 @@ template void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, const dmnType &lDmn, const Matrix &F, const int nfd, const FiberRef &fl, - const double ya_f, const double ya_s, const double ya_n, + const ActiveStress::ActiveTension &active_tension, Matrix &S, Matrix<3 * (nsd - 1)> &Dm, double &Ja) { using namespace consts; using namespace mat_fun; @@ -315,26 +301,13 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, // Active stress from active stress models, already distributed among the // fiber, sheet and sheet-normal directions by the active stress model. - double Tfa = ya_f; // Fiber direction - double Tsa = ya_s; // Sheet direction - double Tna = ya_n; // Sheet-normal direction - - // Validate directional distribution is supported for this constitutive model - // Only Guccione, HO, and HO-ma models support sheet and sheet-normal stress contributions - bool supports_directional_distribution = (stM.isoType == ConstitutiveModelType::stIso_Gucci || - stM.isoType == ConstitutiveModelType::stIso_HO || - stM.isoType == ConstitutiveModelType::stIso_HO_ma); - - if (!supports_directional_distribution && (ya_s > 0.0 || ya_n > 0.0)) { - throw std::runtime_error("Directional distribution of active stress (eta_s > 0 or eta_n > 0) " - "is only supported for Guccione, Holzapfel-Ogden (HO), and Holzapfel-Ogden Modified Anisotropy (HO-ma) models. " - "Current model does not support sheet or sheet-normal stress contributions. " - "Set Fiber_direction=1.0, Sheet_direction=0.0, Sheet_normal_direction=0.0."); - } + double Tfa = active_tension.fibers; // Fiber direction + double Tsa = active_tension.sheets; // Sheet direction + double Tna = active_tension.sheet_normals; // Sheet-normal direction // Aliases for fiber directions const auto& fib_dir1 = fl.col(0); - + // fib_dir2 only exists when nfd >= 2 Eigen::Matrix fib_dir2; if (nfd >= 2) { @@ -357,16 +330,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, Matrix Fa = Matrix::Identity(); Matrix Fai = Fa; - // This commented block implements the active strain formulation, taken from svFSI - // It is commented out because the active strain formulation is not used in the - // current implementation. However, it is left here for reference when we decide to - // implement it. - // if (cep_mod.cem.aStrain) { - // actv_strain(com_mod, cep_mod, ya, nfd, fl, Fa); - // Fai = Fa.inverse(); - // Fe = F * Fai; - // } - Ja = Fa.determinant(); double J = Fe.determinant(); double J2d = pow(J, (-2.0/nd)); @@ -403,7 +366,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, case ConstitutiveModelType::stIso_lin: { double g1 = stM.C10; // mu S += g1*Idm; - return; } break; // St.Venant-Kirchhoff @@ -430,12 +392,9 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, // Compute fictious stress and elasticity tensor Matrix S_bar = 2.0 * stM.C10 * Idm; - Tensor CC_bar; + Tensor CC_bar; CC_bar.setZero(); - // Add fiber reinforcement/active stress - S_bar += Tfa * Hff; - // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -451,9 +410,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, Tensor CC_bar = 4.0 * J4d * stM.C01 * (dyadic_product(Idm, Idm) - fourth_order_identity()); - // Add fiber reinforcement/active stress - S_bar += Tfa * Hff; - // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -492,10 +448,7 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g1 = 4.0*J4d*g1; g2 = 4.0*J4d*g2; Tensor CC_bar = g1 * dyadic_product(Hff_disp, Hff_disp) + g2 * dyadic_product(Hss_disp, Hss_disp); - - // Add fiber reinforcement/active stress - S_bar += Tfa * Hff; - + // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -560,14 +513,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, dyadic_product(RmRm_20, RmRm_20)); CC_bar = r2 * CC_bar; - // Add fiber reinforcement/active stress in all three orthogonal directions - S_bar += Tfa * Hff; // Fiber direction - S_bar += Tsa * Hss; // Sheet direction - if (Tna > 0.0) { - auto Hnn = fib_dir3 * fib_dir3.transpose(); - S_bar += Tna * Hnn; // Sheet-normal direction - } - // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -580,9 +525,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, throw std::runtime_error("[compute_pk2cc] Min fiber directions not defined for Holzapfel material model."); } - // Compute sheet-normal direction - auto fib_dir3 = compute_sheet_normal(fl); - // Compute cross fiber-sheet structure tensor Matrix Hfs = 0.5 * (fib_dir1 * fib_dir2.transpose() + fib_dir2 * fib_dir1.transpose()); @@ -623,11 +565,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g2 = 4.0*J4d*stM.afs*(1.0 + 2.0*stM.bfs*Efs*Efs)* exp(stM.bfs*Efs*Efs); Tensor CC_bar = g1 * dyadic_product(Idm, Idm) + g2 * dyadic_product(Hfs, Hfs); - // 2.S) Add fiber-fiber interaction stress + additional fiber reinforcement/active stress (Tfa) + // 2.S) Add fiber-fiber interaction stress double rexp = exp(stM.bff*Eff*Eff); g1 = c4f * Eff * rexp; g1 = g1 + (0.5*dc4f/stM.bff) * (rexp - 1.0); - g1 = 2.0 * stM.aff * g1 + Tfa; + g1 = 2.0 * stM.aff * g1; S_bar += g1*Hff; // 2.CC) Add fiber-fiber interaction stiffness @@ -637,11 +579,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g1 = 4.0 * J4d * stM.aff * g1; CC_bar += g1*dyadic_product(Hff, Hff); - // 3.S) Add sheet-sheet interaction stress + additional cross-fiber active stress (Tsa) + // 3.S) Add sheet-sheet interaction stress rexp = exp(stM.bss*Ess*Ess); g2 = c4s * Ess * rexp; g2 = g2 + (0.5*dc4s/stM.bss) * (rexp - 1.0); - g2 = 2.0 * stM.ass * g2 + Tsa; + g2 = 2.0 * stM.ass * g2; S_bar += g2 * Hss; // 3.CC) Add sheet-sheet interaction stiffness @@ -651,24 +593,10 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g2 = 4.0 * J4d * stM.ass * g2; CC_bar += g2*dyadic_product(Hss, Hss); - // 4.S) Add sheet-normal active stress (Tna) - if (Tna > 0.0) { - auto Hnn = fib_dir3 * fib_dir3.transpose(); - S_bar += Tna * Hnn; // Sheet-normal direction (fib_dir3 already normalized) - } - // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; CC += CC_iso; - - // Modify S and CC if using active strain - if (cep_mod.cem.aStrain) { - S = Fa * S * Fai.transpose(); - CC_bar = dyadic_product(Fai, Fai); // Reusing CC_bar - CC = double_dot_product(CC, {2,3}, CC_bar, {1,3}); - CC = double_dot_product(CC_bar, {1,3}, CC, {0,1}); - } } break; // HO (Holzapfel-Ogden)-MA model for myocardium with full invariants for the anisotropy terms (modified-anisotropy) @@ -677,9 +605,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, //err = "Min fiber directions not defined for Holzapfel material model (2)" } - // Compute sheet-normal direction - auto fib_dir3 = compute_sheet_normal(fl); - // Compute cross fiber-sheet structure tensor auto Hfs = 0.5 * (fib_dir1 * fib_dir2.transpose() + fib_dir2 * fib_dir1.transpose()); @@ -732,11 +657,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g1 = g1 * 2.0*(1.0 + 2.0*stM.bfs*Efs*Efs); CC += g1*dyadic_product(Hfs, Hfs); - // 2.S) Add fiber-fiber interaction stress + additional reinforcement/active stress (Tfa) + // 2.S) Add fiber-fiber interaction stress double rexp = exp(stM.bff * Eff * Eff); g1 = c4f*Eff*rexp; g1 = g1 + (0.5*dc4f/stM.bff)*(rexp - 1.0); - g1 = (2.0*stM.aff*g1) + Tfa; + g1 = 2.0*stM.aff*g1; S += g1*Hff; // 2.CC) Add fiber-fiber interaction stiffness @@ -746,11 +671,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g1 = 4.0*stM.aff*g1; CC += g1*dyadic_product(Hff, Hff); - // 3.S) Add sheet-sheet interaction stress + additional cross-fiber active stress (Tsa) + // 3.S) Add sheet-sheet interaction stress rexp = exp(stM.bss * Ess * Ess); double g2 = c4s*Ess*rexp; g2 = g2 + (0.5*dc4s/stM.bss)*(rexp - 1.0); - g2 = 2.0*stM.ass*g2 + Tsa; + g2 = 2.0*stM.ass*g2; S += g2*Hss; // 3.CC) Add sheet-sheet interaction stiffness @@ -759,12 +684,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g2 = g2 + (0.5*ddc4s/stM.bss)*(rexp - 1.0); g2 = 4.0*stM.ass*g2; CC += g2*dyadic_product(Hss, Hss); - - // 4.S) Add sheet-normal active stress (Tna) - if (Tna > 0.0) { - auto Hnn = fib_dir3 * fib_dir3.transpose(); - S += Tna * Hnn; // Sheet-normal direction (fib_dir3 already normalized) - } } break; // Universal Material Subroutine - CANN Model @@ -780,8 +699,10 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, std::array,9> ddInv; Matrix N1; - // Compute and store invariants and derivatives wrt C in array of matrices/tensors - CANNModel.computeInvariantsAndDerivatives(C, fl, nfd, J2d, J4d, Ci, Idm, Tfa, N1, psi, Inv, dInv, ddInv); + // Compute and store invariants and derivatives wrt C in array of + // matrices/tensors + CANNModel.computeInvariantsAndDerivatives( + C, fl, nfd, J2d, J4d, Ci, Idm, N1, psi, Inv, dInv, ddInv); // Strain energy function and derivatives CANNModel.evaluate(Inv, psi, dpsi, ddpsi); @@ -790,9 +711,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, S += 2*dInv[i]*dpsi[i]; } - // Fiber reinforcement/active stress - S += Tfa*N1; - // Stiffness Tensor for(int x = 0; x < 9; x++){ CC += 4*dpsi[x]*ddInv[x]; @@ -804,7 +722,72 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, default: throw std::runtime_error("Undefined material constitutive model."); - } + } + + // Active stress. + // + // The sheet and sheet-normal components need a sheet direction to be defined, + // and the sheet-normal one is only defined in 3D (compute_sheet_normal raises + // in 2D). + svmp::check( + nfd >= 1, + "At least one fiber direction must be defined for active stress."); + + const bool has_tangent_f = !utils::is_zero(active_tension.d_fibers); + const bool has_tangent_s = !utils::is_zero(active_tension.d_sheets); + const bool has_tangent_n = !utils::is_zero(active_tension.d_sheet_normals); + + Matrix dS_act = Matrix::Zero(); + + S += Tfa * Hff; + if (has_tangent_f) + dS_act += active_tension.d_fibers * Hff; + + if (!utils::is_zero(Tsa) || !utils::is_zero(active_tension.d_sheets)) { + svmp::check( + nfd >= 2, "Applying active stress along sheets (eta_s > 0) requires a " + "sheet direction, but only " + + std::to_string(nfd) + " fiber directions are defined."); + + S += Tsa * Hss; + if (has_tangent_s) + dS_act += active_tension.d_sheets * Hss; + } + + if (!utils::is_zero(Tna) || !utils::is_zero(active_tension.d_sheet_normals)) { + svmp::check( + nfd >= 2, + "Applying active stress along normals (eta_n > 0) requires both a " + "fiber and a sheet direction, but only " + + std::to_string(nfd) + " fiber directions are defined."); + + auto fib_dir3 = compute_sheet_normal(fl); + const Matrix Hnn = fib_dir3 * fib_dir3.transpose(); + S += Tna * Hnn; + + if (has_tangent_n) + dS_act += active_tension.d_sheet_normals * Hnn; + } + + // Tangent of the active stress. + // + // The active stress depends on the deformation through the fiber stretch + // @f$\lambda = |F f_0| = \sqrt{C : H_{ff}}@f$, so that + // @f$\partial\lambda/\partial C = H_{ff} / (2\lambda)@f$ and + // @f[ + // CC_\text{act} = 2 \frac{\partial S_\text{act}}{\partial C} + // = \frac{1}{\lambda} \frac{\partial S_\text{act}}{\partial\lambda} + // \otimes H_{ff} \;. + // @f] + // + // Only the direct dependence of the active stress on the fiber stretch is + // differentiated here. The active stress also depends on it through the state + // of the active stress model, but differentiating that would mean + // differentiating through the ODE solver of the model. + if (has_tangent_f || has_tangent_s || has_tangent_n) { + const double fiber_stretch = sqrt(fib_dir1.dot(C * fib_dir1)); + CC += (1.0 / fiber_stretch) * dyadic_product(dS_act, Hff); + } // Convert to Voigt Notation cc_to_voigt_eigen(CC, Dm); @@ -813,11 +796,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, // Explicitly instantiate compute_pk2cc for 2D and 3D. template void compute_pk2cc<2>(const ComMod&, const CepMod&, const dmnType&, const Matrix<2>&, const int, const FiberRef<2>&, - const double, const double, const double, Matrix<2>&, Matrix<3>&, double&); + const ActiveStress::ActiveTension &, Matrix<2>&, Matrix<3>&, double&); template void compute_pk2cc<3>(const ComMod&, const CepMod&, const dmnType&, const Matrix<3>&, const int, const FiberRef<3>&, - const double, const double, const double, Matrix<3>&, Matrix<6>&, double&); + const ActiveStress::ActiveTension &, Matrix<3>&, Matrix<6>&, double&); /// @brief Compute 2nd Piola-Kirchhoff stress and material stiffness tensors diff --git a/Code/Source/solver/mat_models.h b/Code/Source/solver/mat_models.h index e5477d610..a85b6a1b3 100644 --- a/Code/Source/solver/mat_models.h +++ b/Code/Source/solver/mat_models.h @@ -4,6 +4,7 @@ #ifndef MAT_MODELS_H #define MAT_MODELS_H +#include "ActiveStress.h" #include "Array.h" #include "CepMod.h" #include "ComMod.h" @@ -45,9 +46,11 @@ using FiberRef = Eigen::Ref>; * @param[in] F Deformation gradient tensor. * @param[in] nfd Number of fiber directions. * @param[in] fl Fiber directions. - * @param[in] ya_f Active tension along the fiber direction. - * @param[in] ya_s Active tension along the sheet direction. - * @param[in] ya_n Active tension along the sheet-normal direction. + * @param[in] active_tension Active tension along the fiber, sheet and + * sheet-normal directions, and its derivatives with respect to the fiber + * stretch. The derivatives contribute the tangent of the active stress to + * @p Dm, which is therefore not symmetric unless the active stress acts along + * the fiber direction alone. * @param[out] S 2nd Piola-Kirchhoff stress tensor (modified in place). * @param[out] Dm Material stiffness tensor (modified in place). * @param[out] Ja Jacobian for active strain @@ -59,9 +62,9 @@ using FiberRef = Eigen::Ref>; template void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, const dmnType &lDmn, const Matrix &F, const int nfd, - const FiberRef &fl, const double ya_f, - const double ya_s, const double ya_n, Matrix &S, - Matrix<3 * (nsd - 1)> &Dm, double &Ja); + const FiberRef &fl, + const ActiveStress::ActiveTension &active_tension, + Matrix &S, Matrix<3 * (nsd - 1)> &Dm, double &Ja); void compute_pk2cc_shlc(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array& fNa0, const Array& gg_0, const Array& gg_x, double& g33, Vector& Sml, Array& Dml); diff --git a/Code/Source/solver/mesh.cpp b/Code/Source/solver/mesh.cpp index b7b8cdfb9..6e90dad6a 100644 --- a/Code/Source/solver/mesh.cpp +++ b/Code/Source/solver/mesh.cpp @@ -32,7 +32,6 @@ void construct_mesh(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const S using namespace consts; - auto& cem = cep_mod.cem; const int nsd = com_mod.nsd; const int tDof = com_mod.tDof; const int dof = com_mod.dof; diff --git a/Code/Source/solver/output.cpp b/Code/Source/solver/output.cpp index 9a321b8b7..0bd0d7e1d 100644 --- a/Code/Source/solver/output.cpp +++ b/Code/Source/solver/output.cpp @@ -233,7 +233,6 @@ void write_restart(Simulation* simulation, std::array& timeP, const So auto& Ad = com_mod.Ad; auto& pS0 = com_mod.pS0; auto& Xion = cep_mod.Xion; - auto& cem = cep_mod.cem; #ifdef debug_write_restart dmsg << "stFileName: " << stFileName; @@ -311,9 +310,6 @@ void write_restart(Simulation* simulation, std::array& timeP, const So } else if (cepEq) { restart_file.write((char*)Ad.data(), Ad.msize()); restart_file.write((char*)Xion.data(), Xion.msize()); - restart_file.write((char*)cem.Ya_f.data(), cem.Ya_f.msize()); - restart_file.write((char*)cem.Ya_s.data(), cem.Ya_s.msize()); - restart_file.write((char*)cem.Ya_n.data(), cem.Ya_n.msize()); } else if (risFlag) { restart_file.write((char*)Ad.data(), Ad.msize()); @@ -336,9 +332,6 @@ void write_restart(Simulation* simulation, std::array& timeP, const So } else if (cepEq) { restart_file.write((char*)Xion.data(), Xion.msize()); - restart_file.write((char*)cem.Ya_f.data(), cem.Ya_f.msize()); - restart_file.write((char*)cem.Ya_s.data(), cem.Ya_s.msize()); - restart_file.write((char*)cem.Ya_n.data(), cem.Ya_n.msize()); } else if (risFlag) { write_ris_data(com_mod, restart_file); diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index ad0a88a0c..015c7d9c2 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -3,6 +3,7 @@ #include "post.h" +#include "ActiveStress.h" #include "Core/Exception.h" #include "FE/Common/FEException.h" #include "FE/Math/DenseLinearAlgebra.h" @@ -802,6 +803,65 @@ void fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, } } +/// @brief Compute active tension along fibers, sheets and sheet normals at +/// every mesh node. +// +void active_tension(const ComMod& com_mod, const int iEq, const mshType& lM, + const Array& lD, Vector& res_f, Vector& res_s, + Vector& res_n) +{ + auto& eq = com_mod.eq[iEq]; + + Vector fiber_stretch(lM.nNo); + if (lM.nFn != 0) { + fib_stretch(com_mod, iEq, lM, lD, fiber_stretch); + } + + res_f = 0.0; + res_s = 0.0; + res_n = 0.0; + + for (int a = 0; a < lM.nNo; a++) { + int Ac = lM.gN(a); + + double Ta_f = 0.0; + double Ta_s = 0.0; + double Ta_n = 0.0; + unsigned int n_domains = 0; + + for (auto& dmn : eq.dmn) { + // Domains whose equations do not allow for active stress (e.g. fluid + // domains) do not contribute to the average, but domains that do + // allow for active stress (e.g. struct) for which active stress is + // not enabled contribute a zero value to the average. + if (!supports_active_stress(dmn.phys)) + continue; + + // Only domains that node Ac actually belongs to contribute to its + // average. Note that if there is only one domain dmnId may not be + // populated, so we only check domain membership if eq.nDmn > 1. + if (eq.nDmn > 1 && !utils::btest(com_mod.dmnId(Ac), dmn.Id)) + continue; + + if (dmn.active_stress != nullptr) { + auto tension = dmn.active_stress->compute_tension( + dmn.active_stress->get_state(Ac), fiber_stretch[a]); + Ta_f += tension.fibers; + Ta_s += tension.sheets; + Ta_n += tension.sheet_normals; + } + + n_domains++; + } + + if (n_domains > 0) { + res_f[a] = Ta_f / n_domains; + res_s[a] = Ta_s / n_domains; + res_n[a] = Ta_n / n_domains; + } + } +} + /// @brief Compute fiber stretch rate dλ/dt via backward finite difference. // void fib_stretch_rate(const ComMod& com_mod, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector& res) @@ -1795,6 +1855,8 @@ void tensor_post_impl(Simulation* simulation, const mshType& lM, const int m, Ar Vector resl(m); Array Nx(nsd,fs.eNoN); Vector N(fs.eNoN); + Vector element_nodes(fs.eNoN); + ActiveStress::Evaluator active_stress_evaluator; int insd = nsd; if (lM.lFib) { @@ -1844,6 +1906,7 @@ void tensor_post_impl(Simulation* simulation, const mshType& lM, const int m, Ar for (int a = 0; a < fs.eNoN; a++) { int Ac = lM.IEN(a,e); + element_nodes(a) = Ac; for (int i = 0; i < nsd; i++) { xl(i,a) = com_mod.x(i,Ac); } @@ -1853,6 +1916,12 @@ void tensor_post_impl(Simulation* simulation, const mshType& lM, const int m, Ar } } + if (eq.dmn[cDmn].active_stress != nullptr) { + active_stress_evaluator.update(*eq.dmn[cDmn].active_stress, element_nodes); + } else { + active_stress_evaluator.clear(); + } + Je = 0.0; double Jac = 0.0; @@ -1950,20 +2019,10 @@ void tensor_post_impl(Simulation* simulation, const mshType& lM, const int m, Ar Matrix sigma = Matrix::Zero(); Matrix S = Matrix::Zero(); - // Interpolate the active stress from active stress models to the - // current Gauss point so that the active contribution is included in - // the reported stress, consistently with the residual assembly. - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - if (eq.dmn[cDmn].active_stress != nullptr) { - for (int a = 0; a < fs.eNoN; a++) { - int Ac = lM.IEN(a,e); - ya_g_f = ya_g_f + N(a)*cep_mod.cem.Ya_f[Ac]; - ya_g_s = ya_g_s + N(a)*cep_mod.cem.Ya_s[Ac]; - ya_g_n = ya_g_n + N(a)*cep_mod.cem.Ya_n[Ac]; - } - } + // Evaluate the active stress at the current Gauss point, the same + // way the residual assembly does, so that the active contribution to + // the reported stress matches the one the solver used. + const auto Ta = active_stress_evaluator.evaluate(eigen_view(N), F, eigen_view(fN)); if (cPhys == EquationType::phys_lElas) { if (nsd == 3) { @@ -1997,7 +2056,7 @@ void tensor_post_impl(Simulation* simulation, const mshType& lM, const int m, Ar Matrix<3*(nsd-1)> Dm; double Ja; mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, - eigen_view(fN), ya_g_f, ya_g_s, ya_g_n, S, Dm, Ja); + eigen_view(fN), Ta, S, Dm, Ja); // TODO: Add viscous stress @@ -2015,7 +2074,7 @@ void tensor_post_impl(Simulation* simulation, const mshType& lM, const int m, Ar Matrix<3*(nsd-1)> Dm; double Ja; mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, - eigen_view(fN), ya_g_f, ya_g_s, ya_g_n, S, Dm, Ja); + eigen_view(fN), Ta, S, Dm, Ja); // TODO: Add viscous stress diff --git a/Code/Source/solver/post.h b/Code/Source/solver/post.h index a147b842a..48cb97d5e 100644 --- a/Code/Source/solver/post.h +++ b/Code/Source/solver/post.h @@ -13,6 +13,9 @@ namespace post { void all_post(Simulation* simulation, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp, const int iEq); +void active_tension(const ComMod& com_mod, const int iEq, const mshType& lM, const Array& lD, + Vector& res_f, Vector& res_s, Vector& res_n); + void bpost(Simulation* simulation, const mshType& lM, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp); diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 6c8c02070..744aefc16 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1340,7 +1340,7 @@ void read_cep_domain(Simulation* simulation, EquationParameters* eq_params, Doma lDmn.cep.odes.relTol = domain_params->relative_tolerance.value(); } - if (domain_params->feedback_parameter_for_stretch_activated_currents.defined() && cep_mod.cem.cpld) { + if (domain_params->feedback_parameter_for_stretch_activated_currents.defined()) { lDmn.cep.Ksac = domain_params->feedback_parameter_for_stretch_activated_currents.value(); } else { lDmn.cep.Ksac = 0.0; @@ -2060,51 +2060,30 @@ void read_files(Simulation* simulation, const std::string& file_name) read_cep_equation(&cep_mod, simulation, eq_params); } } - } - if (cep_mod.cem.cpld) { - if (nEq == 1) { - throw std::runtime_error("Min equations (2) not solved for electro-mechanics coupling"); - } + // Stretch-activated currents require the fiber stretch, which is computed + // only when an equation solving for the displacement is present. + bool has_sac = false; + bool has_displacement_eq = false; - int i = 0; for (int iEq = 0; iEq < nEq; iEq++) { auto& eq = com_mod.eq[iEq]; - if ((eq.phys == EquationType::phys_CEP) || (eq.phys == EquationType::phys_struct) || - (eq.phys == EquationType::phys_ustruct)) { - i = i + 1; - } - } - - if (i != 2) { - throw std::runtime_error("Both electrophysiology and struct have to be solved for electro-mechanics"); - } - if (cep_mod.cem.aStrain) { - if (com_mod.nsd != 3) { - throw std::runtime_error("Active strain coupling is allowed only for 3D bodies"); + if (supports_active_stress(eq.phys)) { + has_displacement_eq = true; } - for (int iEq = 0; iEq < nEq; iEq++) { - auto& eq = com_mod.eq[iEq]; - for (int i = 0; i < eq.nDmn; i++) { - auto& dmn = eq.dmn[i]; - - if ((dmn.phys != EquationType::phys_ustruct) && (dmn.phys != EquationType::phys_struct)) { - continue; - } - - if (dmn.active_stress != nullptr) { - svmp::raise( - "Active strain and active stress cannot be used together."); - } - - if ((dmn.stM.isoType != ConstitutiveModelType::stIso_HO)) { - throw std::runtime_error("Active strain is allowed with Holzapfel-Ogden passive constitutive model only"); - } + for (int iDmn = 0; iDmn < eq.nDmn; iDmn++) { + auto& dmn = eq.dmn[iDmn]; + if ((dmn.phys == EquationType::phys_CEP) && (dmn.cep.Ksac != 0.0)) { + has_sac = true; } } } + + svmp::throw_if(has_sac && !has_displacement_eq, + "Feedback_parameter_for_stretch_activated_currents requires an " + "equation solving for the displacement (struct, ustruct or FSI)."); } // [NOTE] what's going on here? diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index b6dfd718a..2c99d54ce 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -199,7 +199,6 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const dmsg.banner(); #endif - auto& cem = cep_mod.cem; const int nsd = com_mod.nsd; const int tDof = com_mod.tDof; const int dof = com_mod.dof; @@ -229,7 +228,8 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // STRUCT: dof = nsd Vector ptr(eNoN); - Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN), N(eNoN); + Vector pSl(nsymd), N(eNoN); + ActiveStress::Evaluator active_stress_evaluator; Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), Nx(nsd,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN); @@ -252,9 +252,6 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // Create local copies fN = 0.0; pS0l = 0.0; - ya_l_f = 0.0; - ya_l_s = 0.0; - ya_l_n = 0.0; if (lM.fN.size() != 0) { for (int iFn = 0; iFn < nFn; iFn++) { @@ -282,12 +279,12 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const if (pS0.size() != 0) { pS0l.set_col(a, pS0.col(Ac)); } + } - if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l_f(a) = cep_mod.cem.Ya_f[Ac]; - ya_l_s(a) = cep_mod.cem.Ya_s[Ac]; - ya_l_n(a) = cep_mod.cem.Ya_n[Ac]; - } + if (eq.dmn[cDmn].active_stress != nullptr) { + active_stress_evaluator.update(*eq.dmn[cDmn].active_stress, ptr); + } else { + active_stress_evaluator.clear(); } // Gauss integration @@ -316,7 +313,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const if (nsd == 3) { struct_3d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - pS0l, pSl, ya_l_f, ya_l_s, ya_l_n, lR, lK, recompute_visc); + pS0l, pSl, active_stress_evaluator, lR, lK, recompute_visc); #if 0 if (e == 0 && g == 0) { @@ -330,7 +327,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } else if (nsd == 2) { struct_2d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - pS0l, pSl, ya_l_f, ya_l_s, ya_l_n, lR, lK, recompute_visc); + pS0l, pSl, active_stress_evaluator, lR, lK, recompute_visc); } // Prestress @@ -356,8 +353,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, - Vector &pSl, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + Vector &pSl, const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, const bool recompute_visc) { using namespace consts; using namespace mat_fun; @@ -406,36 +402,33 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // Inertia, damping and body force: the term the residual weights with N const Eigen::Vector2d ud = (rho*(acc - bfm) + dmp*vel) * Nm - rho * fb; - // Active stress activation along fiber, sheet and sheet-normal - const double ya_g_f = eigen_view(ya_l_f).dot(Nm); - const double ya_g_s = eigen_view(ya_l_s).dot(Nm); - const double ya_g_n = eigen_view(ya_l_n).dot(Nm); - // Prestress at this Gauss point, in Voigt order [11, 22, 12] const Eigen::Vector pS0g = eigen_view<3>(pS0l) * Nm; Matrix<2> S0; S0 << pS0g(0), pS0g(2), pS0g(2), pS0g(1); - - #ifdef debug_struct_2d - dmsg << "ud: " << ud(0) << " " << ud(1); - dmsg << "F: " << F(0,0); - dmsg << "ya_g_f: " << ya_g_f; - dmsg << "ya_g_s: " << ya_g_s; - dmsg << "ya_g_n: " << ya_g_n; -#endif // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) const Matrix<2> vx = vel * Nxm.transpose(); const Matrix<2> F = Matrix<2>::Identity() + disp * Nxm.transpose(); + // Active tension, evaluated here from the fiber stretch of F. + const auto Ta = active_stress_evaluator.evaluate<2>(Nm, F, eigen_view<2>(fN)); + + #ifdef debug_struct_2d + dmsg << "ud: " << ud(0) << " " << ud(1); + dmsg << "F: " << F(0,0); + dmsg << "Ta.fibers: " << Ta.fibers; + dmsg << "Ta.sheets: " << Ta.sheets; + dmsg << "Ta.sheet_normals: " << Ta.sheet_normals; + #endif + // 2nd Piola-Kirchhoff stress (S) and material stiffness tensor in Voight notation (Dm) Matrix<2> S; Matrix<3> Dm; double Ja; - mat_models::compute_pk2cc<2>(com_mod, cep_mod, dmn, F, nFn, eigen_view<2>(fN), ya_g_f, ya_g_s, - ya_g_n, S, Dm, Ja); + mat_models::compute_pk2cc<2>(com_mod, cep_mod, dmn, F, nFn, eigen_view<2>(fN), Ta, S, Dm, Ja); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. // Reuse from the previous Gauss point when shape function gradients @@ -520,8 +513,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, - Vector &pSl, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + Vector &pSl, const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, const bool recompute_visc) { using namespace consts; using namespace mat_fun; @@ -575,11 +567,6 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // Inertia, damping and body force. const Eigen::Vector3d ud = (rho*(acc - bfm) + dmp*vel) * Nm - rho * fb; - // Active stress activation along fiber, sheet and sheet-normal - const double ya_g_f = eigen_view(ya_l_f).dot(Nm); - const double ya_g_s = eigen_view(ya_l_s).dot(Nm); - const double ya_g_n = eigen_view(ya_l_n).dot(Nm); - // Prestress at this Gauss point, in Voigt order [11, 22, 33, 12, 23, 31] const Eigen::Vector pS0g = eigen_view<6>(pS0l) * Nm; @@ -592,14 +579,16 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Matrix<3> vx = vel * Nxm.transpose(); const Matrix<3> F = Matrix<3>::Identity() + disp * Nxm.transpose(); + // Active tension, evaluated here from the fiber stretch of F. + const auto Ta = active_stress_evaluator.evaluate<3>(Nm, F, eigen_view<3>(fN)); + // 2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in // Voigt notation (Dm) // Matrix<3> S; Matrix<6> Dm; double Ja; - mat_models::compute_pk2cc<3>(com_mod, cep_mod, dmn, F, nFn, eigen_view<3>(fN), ya_g_f, ya_g_s, - ya_g_n, S, Dm, Ja); + mat_models::compute_pk2cc<3>(com_mod, cep_mod, dmn, F, nFn, eigen_view<3>(fN), Ta, S, Dm, Ja); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. // Reuse from the previous Gauss point when shape function gradients diff --git a/Code/Source/solver/sv_struct.h b/Code/Source/solver/sv_struct.h index 89d082d71..eb76db1da 100644 --- a/Code/Source/solver/sv_struct.h +++ b/Code/Source/solver/sv_struct.h @@ -4,6 +4,7 @@ #ifndef STRUCT_H #define STRUCT_H +#include "ActiveStress.h" #include "ComMod.h" #include "SolutionStates.h" @@ -24,8 +25,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, - Vector &pSl, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + Vector &pSl, const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, const bool recompute_visc); void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, @@ -33,8 +33,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, - Vector &pSl, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + Vector &pSl, const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, const bool recompute_visc); }; diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 673c77d7e..54be7cd89 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -32,8 +32,8 @@ namespace ustruct { -void b_ustruct_2d(const ComMod& com_mod, const int eNoN, const double w, const Vector& N, - const Array& Nx, const Array& dl, const Vector& hl, const Vector& nV, +void b_ustruct_2d(const ComMod& com_mod, const int eNoN, const double w, const Vector& N, + const Array& Nx, const Array& dl, const Vector& hl, const Vector& nV, Array& lR, Array3& lK, Array3& lKd) { int cEq = com_mod.cEq; @@ -45,10 +45,10 @@ void b_ustruct_2d(const ComMod& com_mod, const int eNoN, const double w, const V int i = eq.s; int j = i + 1; - Vector nFi(2); + Vector nFi(2); Array NxFi(2,eNoN); - Array F(2,2); + Array F(2,2); F(0,0) = 1.0; F(1,1) = 1.0; @@ -90,8 +90,8 @@ void b_ustruct_2d(const ComMod& com_mod, const int eNoN, const double w, const V } /// @brief Add follower pressure load contributions to the local residual and stiffness matrix. -/// @param com_mod -/// @param eNoN +/// @param com_mod +/// @param eNoN /// @param w Gauss point weight times reference configuration area /// @param N Shape function values at the Gauss point /// @param Nx Shape function derivatives at the Gauss point @@ -101,8 +101,8 @@ void b_ustruct_2d(const ComMod& com_mod, const int eNoN, const double w, const V /// @param lR Local residual /// @param lK Local stiffness matrix /// @param lKd Local stiffness matrix (displacement) -void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const Vector& N, - const Array& Nx, const Array& dl, const Vector& hl, const Vector& nV, +void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const Vector& N, + const Array& Nx, const Array& dl, const Vector& hl, const Vector& nV, Array& lR, Array3& lK, Array3& lKd) { int cEq = com_mod.cEq; @@ -115,10 +115,10 @@ void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const V int j = i + 1; int k = j + 1; - Vector nFi(3); + Vector nFi(3); Array NxFi(3,eNoN); - Array F(3,3); + Array F(3,3); F(0,0) = 1.0; F(1,1) = 1.0; F(2,2) = 1.0; @@ -186,17 +186,17 @@ void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const V /// @brief Check is a constitutive model is valid for the ustruct equation. // bool constitutive_model_is_valid(consts::ConstitutiveModelType model) -{ +{ using namespace consts; - - static std::set unsupported_models { + + static std::set unsupported_models { ConstitutiveModelType::stIso_lin, ConstitutiveModelType::stIso_StVK, ConstitutiveModelType::stIso_mStVK }; return unsupported_models.count(model) == 0; -} +} /// @brief Reproduces Fortran CONSTRUCT_uSOLID. // @@ -215,7 +215,6 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const dmsg << "lM.nFs: " << lM.nFs; #endif - auto& cem = cep_mod.cem; const int nsd = com_mod.nsd; const int tDof = com_mod.tDof; const int dof = com_mod.dof; @@ -249,10 +248,11 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // USTRUCT: dof = nsd+1 Vector ptr(eNoN); - Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN), N(eNoN); + Vector pSl(nsymd), N(eNoN); Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), Nx(nsd,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN); + ActiveStress::Evaluator active_stress_evaluator; for (int e = 0; e < lM.nEl; e++) { // Change the current domain which will be used in later function calls. @@ -264,9 +264,6 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // Create local copies fN = 0.0; - ya_l_f = 0.0; - ya_l_s = 0.0; - ya_l_n = 0.0; if (lM.fN.size() != 0) { for (int iFn = 0; iFn < nFn; iFn++) { @@ -290,12 +287,12 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const dl(i,a) = Dg(i,Ac); yl(i,a) = Yg(i,Ac); } + } - if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l_f(a) = cep_mod.cem.Ya_f[Ac]; - ya_l_s(a) = cep_mod.cem.Ya_s[Ac]; - ya_l_n(a) = cep_mod.cem.Ya_n[Ac]; - } + if (eq.dmn[cDmn].active_stress != nullptr) { + active_stress_evaluator.update(*eq.dmn[cDmn].active_stress, ptr); + } else { + active_stress_evaluator.clear(); } // Initialize residual and tangents @@ -345,15 +342,15 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); ustruct_3d_m(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, nFn, w, - Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, ya_l_f, ya_l_s, - ya_l_n, lR, lK, lKd, recompute_visc); + Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, + active_stress_evaluator, lR, lK, lKd, recompute_visc); } else if (nsd == 2) { auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); ustruct_2d_m(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, nFn, w, - Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, ya_l_f, ya_l_s, - ya_l_n, lR, lK, lKd, recompute_visc); + Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, + active_stress_evaluator, lR, lK, lKd, recompute_visc); } } // for g = 0 to fs[0].nG @@ -385,13 +382,13 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const if (nsd == 3) { auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); - ustruct_3d_c(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, Jac, N0, N1, Nwx, + ustruct_3d_c(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, Jac, N0, N1, Nwx, Nqx, al, yl, dl, bfl, lR, lK, lKd); } else if (nsd == 2) { auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); - ustruct_2d_c(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, Jac, N0, N1, Nwx, + ustruct_2d_c(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, Jac, N0, N1, Nwx, Nqx, al, yl, dl, bfl, lR, lK, lKd); } @@ -411,25 +408,25 @@ int get_col_ptr(ComMod& com_mod, const int rowN, const int colN) int left = rowPtr(rowN); int right = rowPtr(rowN+1); int ptr = (right + left) / 2; - + while (colN != colPtr(ptr)) { if (colN > colPtr(ptr)) { left = ptr; - } else { + } else { right = ptr; } ptr = (right + left) / 2; } return ptr; -} +} /// @brief Reproduces Fortran USTRUCT2D_C. // void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, const double w, const double Je, const Vector& Nw, const Vector& Nq, - const Array& Nwx, const Array& Nqx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, Array& lR, Array3& lK, + const Array& Nwx, const Array& Nqx, const Array& al, const Array& yl, + const Array& dl, const Array& bfl, Array& lR, Array3& lK, Array3& lKd) { using namespace consts; @@ -489,7 +486,7 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in double Jac = F.determinant(); const Matrix<2> Fi = F.inverse(); - // Pressure and its gradients + // Pressure and its gradients // double p = 0.0; double pd = 0.0; @@ -556,7 +553,7 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in for (int a = 0; a < eNoNq; a++) { NxNx = NqxFi.col(a).dot(NwxFi.col(b)); - // dC/dV_1 + af/am *dC/dU_1 + // dC/dV_1 + af/am *dC/dU_1 // T0 = Nq(a)*(rC*NwxFi(0,b) - VxNwx(0,b)); T1 = tauM*(rMNqx(a)*NwxFi(0,b) - rMNwx(b)*NqxFi(0,a)); @@ -567,7 +564,7 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = (am*tauM*rho)*NqxFi(0,a)*Nw(b) + af*Nq(a)*NwxFi(0,b); lK(6,a,b) += w*Jac*T1 + afm*Ku; - // dC/dV_2 + af/am *dC/dU_2 + // dC/dV_2 + af/am *dC/dU_2 T0 = Nq(a)*(rC*NwxFi(1,b) - VxNwx(1,b)); T1 = tauM*(rMNqx(a)*NwxFi(1,b) - rMNwx(b)*NqxFi(1,a)); T2 = -tauM*NxNx*PxFi(1); @@ -595,8 +592,8 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in // void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, const double w, const double Je, const Vector& Nw, const Vector& Nq, - const Array& Nwx, const Array& Nqx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, Array& lR, Array3& lK, + const Array& Nwx, const Array& Nqx, const Array& al, const Array& yl, + const Array& dl, const Array& bfl, Array& lR, Array3& lK, Array3& lKd) { using namespace consts; @@ -657,7 +654,7 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in double Jac = F.determinant(); const Matrix<3> Fi = F.inverse(); - // Pressure and its gradients + // Pressure and its gradients // double p = 0.0; double pd = 0.0; @@ -726,7 +723,7 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in for (int a = 0; a < eNoNq; a++) { NxNx = NqxFi.col(a).dot(NwxFi.col(b)); - // dC/dV_1 + af/am *dC/dU_1 + // dC/dV_1 + af/am *dC/dU_1 // T0 = Nq(a)*(rC*NwxFi(0,b) - VxNwx(0,b)); T1 = tauM*(rMNqx(a)*NwxFi(0,b) - rMNwx(b)*NqxFi(0,a)); @@ -737,7 +734,7 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = (am*tauM*rho)*NqxFi(0,a)*Nw(b) + af*Nq(a)*NwxFi(0,b); lK(12,a,b) += w*Jac*T1 + afm*Ku; - // dC/dV_2 + af/am *dC/dU_2 + // dC/dV_2 + af/am *dC/dU_2 T0 = Nq(a)*(rC*NwxFi(1,b) - VxNwx(1,b)); T1 = tauM*(rMNqx(a)*NwxFi(1,b) - rMNwx(b)*NqxFi(1,a)); T2 = -tauM*NxNx*PxFi(1); @@ -747,7 +744,7 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = (am*tauM*rho)*NqxFi(1,a)*Nw(b) + af*Nq(a)*NwxFi(1,b); lK(13,a,b) += w*Jac*T1 + afm*Ku; - // dC/dV_3 + af/am *dC/dU_3 + // dC/dV_3 + af/am *dC/dU_3 // T0 = Nq(a)*(rC*NwxFi(2,b) - VxNwx(2,b)); T1 = tauM*(rMNqx(a)*NwxFi(2,b) - rMNwx(b)*NqxFi(2,a)); @@ -780,8 +777,8 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Vector &Nq, const Array &Nwx, const Array &al, const Array &yl, const Array &dl, const Array &bfl, - const Array &fN, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + const Array &fN, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, Array3 &lKd, const bool recompute_visc) { using namespace consts; @@ -835,11 +832,6 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Eigen::Vector2d v = vel * Nwm; const Eigen::Vector2d vd = (acc - bfm) * Nwm - fb; - // Active stress activation along fiber, sheet and sheet-normal - const double ya_g_f = eigen_view(ya_l_f).dot(Nwm); - const double ya_g_s = eigen_view(ya_l_s).dot(Nwm); - const double ya_g_n = eigen_view(ya_l_n).dot(Nwm); - // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) const Matrix<2> vx = vel * Nwxm.transpose(); const Matrix<2> F = Matrix<2>::Identity() + disp * Nwxm.transpose(); @@ -847,6 +839,9 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, double Jac = F.determinant(); const Matrix<2> Fi = F.inverse(); + // Active tension, evaluated here from the fiber stretch of F. + const auto Ta = active_stress_evaluator.evaluate<2>(Nwm, F, eigen_view<2>(fN)); + // Pressure and its time derivative // double p = 0.0; @@ -862,8 +857,8 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, Matrix<2> Siso; Matrix<3> Dm; double Ja = 0; - mat_models::compute_pk2cc<2>(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, eigen_view<2>(fN), ya_g_f, - ya_g_s, ya_g_n, Siso, Dm, Ja); + mat_models::compute_pk2cc<2>(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, eigen_view<2>(fN), + Ta, Siso, Dm, Ja); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. Reuse the // previous Gauss point's when shape function gradients are constant within an @@ -930,7 +925,7 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // Tangent (stiffness) matrices // - double NxSNx{0.0}, BtDB{0.0}, + double NxSNx{0.0}, BtDB{0.0}, T1{0.0}, T2{0.0}, T3{0.0}, Tv{0.0}, Ku{0.0}; @@ -1026,8 +1021,8 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Vector &Nq, const Array &Nwx, const Array &al, const Array &yl, const Array &dl, const Array &bfl, - const Array &fN, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + const Array &fN, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, Array3 &lKd, const bool recompute_visc) { using namespace consts; @@ -1083,11 +1078,6 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Eigen::Vector3d v = vel * Nwm; const Eigen::Vector3d vd = (acc - bfm) * Nwm - fb; - // Active stress activation along fiber, sheet and sheet-normal - const double ya_g_f = eigen_view(ya_l_f).dot(Nwm); - const double ya_g_s = eigen_view(ya_l_s).dot(Nwm); - const double ya_g_n = eigen_view(ya_l_n).dot(Nwm); - // Velocity and deformation gradients: Grad(v) and F = I + Grad(u) const Matrix<3> vx = vel * Nwxm.transpose(); const Matrix<3> F = Matrix<3>::Identity() + disp * Nwxm.transpose(); @@ -1095,6 +1085,9 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, double Jac = F.determinant(); const Matrix<3> Fi = F.inverse(); + // Active tension, evaluated here from the fiber stretch of F. + const auto Ta = active_stress_evaluator.evaluate<3>(Nwm, F, eigen_view<3>(fN)); + // Pressure and its time derivative // double p = 0.0; @@ -1111,8 +1104,8 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, Matrix<3> Siso; Matrix<6> Dm; double Ja = 0; - mat_models::compute_pk2cc<3>(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, eigen_view<3>(fN), ya_g_f, - ya_g_s, ya_g_n, Siso, Dm, Ja); + mat_models::compute_pk2cc<3>(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, eigen_view<3>(fN), + Ta, Siso, Dm, Ja); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. Reuse the // previous Gauss point's when shape function gradients are constant within an @@ -1134,7 +1127,7 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, if (vmsFlag) { mat_models::compute_tau(com_mod, eq.dmn[cDmn], Jac, Je, tauM, tauC); - } else { + } else { tauM = 0.0; tauC = 0.0; } @@ -1198,10 +1191,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, BtDB = Bm[a].col(0).dot(DBm.col(0)); T1 = Jac*rho*vd(0)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(0,b); - + Ku = w*af*(T1 + T2 + BtDB + NxSNx + visc.du(0,a,b)); lKd(0,a,b) += Ku; - + T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(0,a)*NxFi(0,b); Tv = af*visc.dv(0,a,b); @@ -1212,10 +1205,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T1 = Jac*rho*vd(0)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(1,b); T3 = Jac*rCl*(NxFi(0,a)*NxFi(1,b) - NxFi(1,a)*NxFi(0,b)); - + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(1,a,b)); lKd(1,a,b) += Ku; - + T2 = af*Jac*tauC*rho*NxFi(0,a)*NxFi(1,b); Tv = af*visc.dv(1,a,b); lK(1,a,b) += w*(T2 + Tv) + afm*Ku; @@ -1226,10 +1219,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T1 = Jac*rho*vd(0)*Nw(a)*NxFi(2,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(2,b); T3 = Jac*rCl*(NxFi(0,a)*NxFi(2,b) - NxFi(2,a)*NxFi(0,b)); - + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(2,a,b)); lKd(2,a,b) += Ku; - + T2 = af*Jac*tauC*rho*NxFi(0,a)*NxFi(2,b); Tv = af*visc.dv(2,a,b); lK(2,a,b) += w*(T2 + Tv) + afm*Ku; @@ -1241,10 +1234,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T1 = Jac*rho*vd(1)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(0,b); T3 = Jac*rCl*(NxFi(1,a)*NxFi(0,b) - NxFi(0,a)*NxFi(1,b)); - + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(3,a,b)); lKd(3,a,b) += Ku; - + T2 = af*Jac*tauC*rho*NxFi(1,a)*NxFi(0,b); Tv = af*visc.dv(3,a,b); @@ -1258,10 +1251,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T2 = -tauC*Jac*NxFi(1,a)*VxNx(1,b); - + Ku = w*af*(T1 + T2 + BtDB + NxSNx + visc.du(4,a,b)); lKd(4,a,b) += Ku; - + T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(1,a)*NxFi(1,b); Tv = af*visc.dv(4,a,b); @@ -1275,10 +1268,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T2 = -tauC*Jac*NxFi(1,a)*VxNx(2,b); T3 = Jac*rCl*(NxFi(1,a)*NxFi(2,b) - NxFi(2,a)*NxFi(1,b)); - + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(5,a,b)); lKd(5,a,b) += Ku; - + T2 = af*Jac*tauC*rho*NxFi(1,a)*NxFi(2,b); Tv = af*visc.dv(5,a,b); lK(6,a,b) += w*(T2 + Tv) + afm*Ku; @@ -1290,10 +1283,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T1 = Jac*rho*vd(2)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(0,b); T3 = Jac*rCl*(NxFi(2,a)*NxFi(0,b) - NxFi(0,a)*NxFi(2,b)); - + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(6,a,b)); lKd(6,a,b) += Ku; - + T2 = af*Jac*tauC*rho*NxFi(2,a)*NxFi(0,b); Tv = af*visc.dv(6,a,b); lK(8,a,b) += w*(T2 + Tv) + afm*Ku; @@ -1305,10 +1298,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T1 = Jac*rho*vd(2)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(1,b); T3 = Jac*rCl*(NxFi(2,a)*NxFi(1,b) - NxFi(1,a)*NxFi(2,b)); - + Ku = w*af*(T1 + T2 + T3 + BtDB + visc.du(7,a,b)); lKd(7,a,b) += Ku; - + T2 = af*Jac*tauC*rho*NxFi(2,a)*NxFi(1,b); Tv = af*visc.dv(7,a,b); @@ -1320,10 +1313,10 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, T1 = Jac*rho*vd(2)*Nw(a)*NxFi(2,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(2,b); - + Ku = w*af*(T1 + T2 + BtDB + NxSNx + visc.du(8,a,b)); lKd(8,a,b) += Ku; - + T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(2,a)*NxFi(2,b); Tv = af*visc.dv(8,a,b); @@ -1356,7 +1349,7 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // /// @brief Replicates 'SUBROUTINE USTRUCT_DOASSEM(d, eqN, lKd, lK, lR)' // -void ustruct_do_assem(ComMod& com_mod, const int d, const Vector& eqN, const Array3& lKd, +void ustruct_do_assem(ComMod& com_mod, const int d, const Vector& eqN, const Array3& lKd, const Array3& lK, const Array& lR) { const int nsd = com_mod.nsd; @@ -1438,7 +1431,7 @@ void ustruct_do_assem(ComMod& com_mod, const int d, const Vector& eqN, cons } } - } else { + } else { // Stiffness matrix (A) is assembled using mapped rows and columns // Gradient matrix (B) is assembled using mapped row but unmapped @@ -1508,8 +1501,8 @@ void ustruct_r(ComMod& com_mod, const SolutionStates& solutions) const auto& Yg = solutions.intermediate.get_velocity(); using namespace consts; - #define n_debug_ustruct_r - #ifdef debug_ustruct_r + #define n_debug_ustruct_r + #ifdef debug_ustruct_r DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); #endif @@ -1518,7 +1511,7 @@ void ustruct_r(ComMod& com_mod, const SolutionStates& solutions) const auto& eq = com_mod.eq[cEq]; if (eq.phys != EquationType::phys_ustruct && eq.phys != EquationType::phys_FSI) { - return; + return; } const int tnNo = com_mod.tnNo; @@ -1535,7 +1528,7 @@ void ustruct_r(ComMod& com_mod, const SolutionStates& solutions) int s = eq.s; double amg = (eq.gam - eq.am) / (eq.gam - 1.0); double ami = 1.0 / eq.am; - #ifdef debug_ustruct_r + #ifdef debug_ustruct_r dmsg << "nsd: " << nsd; dmsg << "s: " << s; dmsg << "eq.itr: " << eq.itr; @@ -1546,7 +1539,7 @@ void ustruct_r(ComMod& com_mod, const SolutionStates& solutions) if (eq.itr > 1) { Rd = 0.0; } else { - for (int a = 0; a < tnNo; a++) { + for (int a = 0; a < tnNo; a++) { if (!all_fun::is_domain(com_mod, eq, a, EquationType::phys_ustruct)) { continue; } @@ -1558,7 +1551,7 @@ void ustruct_r(ComMod& com_mod, const SolutionStates& solutions) if (nsd == 3) { Array KU(4,tnNo); - for (int a = 0; a < tnNo; a++) { + for (int a = 0; a < tnNo; a++) { if (!all_fun::is_domain(com_mod, eq, a, EquationType::phys_ustruct)) { continue; } @@ -1575,7 +1568,7 @@ void ustruct_r(ComMod& com_mod, const SolutionStates& solutions) all_fun::commu(com_mod, KU); - for (int a = 0; a < tnNo; a++) { + for (int a = 0; a < tnNo; a++) { R(0,a) = R(0,a) - ami*KU(0,a); R(1,a) = R(1,a) - ami*KU(1,a); R(2,a) = R(2,a) - ami*KU(2,a); @@ -1584,7 +1577,7 @@ void ustruct_r(ComMod& com_mod, const SolutionStates& solutions) } else { Array KU(3,tnNo); - for (int a = 0; a < tnNo; a++) { + for (int a = 0; a < tnNo; a++) { if (!all_fun::is_domain(com_mod, eq, a, EquationType::phys_ustruct)) { continue; } @@ -1599,14 +1592,14 @@ void ustruct_r(ComMod& com_mod, const SolutionStates& solutions) all_fun::commu(com_mod, KU); - for (int a = 0; a < tnNo; a++) { + for (int a = 0; a < tnNo; a++) { R(0,a) = R(0,a) - ami*KU(0,a); R(1,a) = R(1,a) - ami*KU(1,a); R(2,a) = R(2,a) - ami*KU(2,a); } } - } -} + } +} }; diff --git a/Code/Source/solver/ustruct.h b/Code/Source/solver/ustruct.h index c89e99c61..a6838c263 100644 --- a/Code/Source/solver/ustruct.h +++ b/Code/Source/solver/ustruct.h @@ -4,6 +4,7 @@ #ifndef USTRUCT_H #define USTRUCT_H +#include "ActiveStress.h" #include "ComMod.h" #include "SolutionStates.h" @@ -35,8 +36,7 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Vector &Nq, const Array &Nwx, const Array &al, const Array &yl, const Array &dl, const Array &bfl, - const Array &fN, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + const Array &fN, const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, Array3 &lKd, const bool recompute_visc); @@ -52,8 +52,7 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Vector &Nq, const Array &Nwx, const Array &al, const Array &yl, const Array &dl, const Array &bfl, - const Array &fN, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + const Array &fN, const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, Array3 &lKd, const bool recompute_visc); diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index ab89e1456..264fff0a6 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1066,6 +1066,23 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b for (int iEq = 0; iEq < nEq; iEq++) { auto& eq = eqs[iEq]; + // Active tension along fibers, sheets and sheet normals, computed + // lazily the first time any of the three is requested below, and + // reused for the others. + Vector active_tension_f, active_tension_s, active_tension_n; + bool active_tension_computed = false; + auto compute_active_tension = [&]() { + if (active_tension_computed) { + return; + } + active_tension_f.resize(msh.nNo); + active_tension_s.resize(msh.nNo); + active_tension_n.resize(msh.nNo); + post::active_tension(simulation->com_mod, iEq, msh, solutions.current.get_displacement(), + active_tension_f, active_tension_s, active_tension_n); + active_tension_computed = true; + }; + for (int iOut = 0; iOut < eq.nOutput; iOut++) { if (!eq.output[iOut].options.spatial) { continue; @@ -1335,23 +1352,23 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b } break; case OutputNameType::outGrp_activeTensionFibers: { + compute_active_tension(); for (int a = 0; a < msh.nNo; a++) { - int Ac = msh.gN(a); - d[iM].x(is, a) = simulation->cep_mod.cem.Ya_f[Ac]; + d[iM].x(is, a) = active_tension_f[a]; } } break; case OutputNameType::outGrp_activeTensionSheets: { + compute_active_tension(); for (int a = 0; a < msh.nNo; a++) { - int Ac = msh.gN(a); - d[iM].x(is, a) = simulation->cep_mod.cem.Ya_s[Ac]; + d[iM].x(is, a) = active_tension_s[a]; } } break; case OutputNameType::outGrp_activeTensionNormal: { + compute_active_tension(); for (int a = 0; a < msh.nNo; a++) { - int Ac = msh.gN(a); - d[iM].x(is, a) = simulation->cep_mod.cem.Ya_n[Ac]; + d[iM].x(is, a) = active_tension_n[a]; } } break; diff --git a/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu index c423af8b6..e327f6d18 100644 --- a/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu +++ b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:d09218c7505190adf4963e2c080d4aeab5712d7f22376bfd53abc13388d7e7a4 -size 1441592 +oid sha256:7fc5f9fa247860346783ef5ac147f1adc0072b744281d4aa85e87251d4f98867 +size 1476088 diff --git a/tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu b/tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu index 3ab600fb8..e850a8c89 100644 --- a/tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu +++ b/tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:906c4b394fe4e0901f5b0fbd7ebc8b55b96bf39faf5b4a75d81efa40ea61f36a -size 76838 +oid sha256:dc56180243bd9c7df9e0cced187be1dc91163fb01787ee12097d007e1add4f3e +size 76532 diff --git a/tests/cases/ustruct/LV_Guccione_active/result_001.vtu b/tests/cases/ustruct/LV_Guccione_active/result_001.vtu index 1d7dee31a..e1adb7b16 100644 --- a/tests/cases/ustruct/LV_Guccione_active/result_001.vtu +++ b/tests/cases/ustruct/LV_Guccione_active/result_001.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3c16cafb77f58a0a3d54f42d5d84bbc18ba9fad98d826533f80f7ea5de4eea06 -size 73482 +oid sha256:3df354be4fcf0d0543b20760f1b1048644b3d2e2a7c38f2f14efc50effca7441 +size 73036 diff --git a/tests/unitTests/active_stress_tests/active_stress_test_helpers.h b/tests/unitTests/active_stress_tests/active_stress_test_helpers.h index e83fe583f..00c29247e 100644 --- a/tests/unitTests/active_stress_tests/active_stress_test_helpers.h +++ b/tests/unitTests/active_stress_tests/active_stress_test_helpers.h @@ -46,6 +46,13 @@ struct ActiveStressTrajectoryConfiguration { /// Base tolerance in the pointwise comparison formula. double tolerance = 1.0e-10; + + /// Fiber-stretch increment of the central finite difference used to check + /// the active tension derivative. + double derivative_step = 1.0e-6; + + /// Base tolerance in the active tension derivative comparison formula. + double derivative_tolerance = 1.0e-8; }; /** @@ -152,6 +159,18 @@ struct ActiveStressTrajectoryConfiguration { * Checkpoint @c N is compared after the update starting at * @f$ t = N \Delta t @f$, which advances the interval * @f$[N\Delta t,(N+1)\Delta t]@f$. + * + * ### Active tension derivative + * + * After every update, the partial derivative of the active tension with + * respect to the fiber stretch, at fixed state, is compared with the central + * finite difference + * @f[ + * \frac{T_a(\mathbf{s}, \lambda + h) - T_a(\mathbf{s}, \lambda - h)}{2h}, + * @f] + * where @f$\mathbf{s}@f$ is the updated state and @f$h@f$ is + * @c derivative_step. This check needs no reference data, and is performed at + * every simulation step rather than only at the checkpoints. */ template class ActiveStressTrajectoryTest { @@ -164,6 +183,7 @@ class ActiveStressTrajectoryTest { public: using ConcreteModel::advance_time_step_local; using ConcreteModel::compute_active_tension_local; + using ConcreteModel::compute_active_tension_derivative_local; using ConcreteModel::init_local; using ConcreteModel::read_model_specific_parameters; }; @@ -218,6 +238,19 @@ class ActiveStressTrajectoryTest { const double active_tension = model_->compute_active_tension_local(state, fiber_stretch); + const double h = configuration_.derivative_step; + const double active_tension_derivative = + model_->compute_active_tension_derivative_local(state, fiber_stretch); + const double finite_difference = + (model_->compute_active_tension_local(state, fiber_stretch + h) - + model_->compute_active_tension_local(state, fiber_stretch - h)) / + (2.0 * h); + EXPECT_NEAR(active_tension_derivative, finite_difference, + configuration_.derivative_tolerance * + (1.0 + std::fabs(finite_difference))) + << "trajectory update index " << step + << ", active tension derivative with respect to fiber stretch"; + compare_checkpoint_if_present(step, state, active_tension, checkpoint_index); } @@ -255,6 +288,14 @@ class ActiveStressTrajectoryTest { !(configuration_.tolerance > 0.0) || !std::isfinite(configuration_.tolerance), "ActiveStress test tolerance must be finite and positive"); + svmp::throw_if( + !(configuration_.derivative_step > 0.0) || + !std::isfinite(configuration_.derivative_step), + "ActiveStress test derivative step must be finite and positive"); + svmp::throw_if( + !(configuration_.derivative_tolerance > 0.0) || + !std::isfinite(configuration_.derivative_tolerance), + "ActiveStress test derivative tolerance must be finite and positive"); svmp::throw_if( configuration_.reference_csv_filename.empty(), "ActiveStress trajectory reference CSV filename must not be empty"); diff --git a/tests/unitTests/material_model_tests/test_material_common.h b/tests/unitTests/material_model_tests/test_material_common.h index 1ef07b1da..b9f449d61 100644 --- a/tests/unitTests/material_model_tests/test_material_common.h +++ b/tests/unitTests/material_model_tests/test_material_common.h @@ -300,9 +300,7 @@ class TestMaterialModel : public TestBase { public: int nFn; Array fN; - double ya_g_f; - double ya_g_s; - double ya_g_n; + ActiveStress::ActiveTension active_tension; bool ustruct; TestMaterialModel(const consts::ConstitutiveModelType matType, const consts::ConstitutiveModelType penType) { @@ -317,9 +315,7 @@ class TestMaterialModel : public TestBase { // Initialize fibers and other material parameters nFn = 2; // Number of fiber directions fN = Array(nsd, nFn); // Fiber directions array (initialized to zeros) - ya_g_f = 0.0; // Active tension along fibers. - ya_g_s = 0.0; // Active tension along sheets. - ya_g_n = 0.0; // Active tension along sheet normals. + active_tension = ActiveStress::ActiveTension{}; // No active tension. // Flag to use struct or ustruct material models // If struct, calls compute_pk2cc() and uses strain energy composed of isochoric and volumetric parts @@ -365,7 +361,7 @@ class TestMaterialModel : public TestBase { mat_models::compute_pk2cc<3>(com_mod, cep_mod, dmn, convert_to_eigen_matrix>(F), nFn, - eigen_view<3>(fN), ya_g_f, ya_g_s, ya_g_n, S_e, Dm_e, J); + eigen_view<3>(fN), active_tension, S_e, Dm_e, J); convert_to_array(S_e, S); copy_Dm(Dm_e, Dm); @@ -375,7 +371,7 @@ class TestMaterialModel : public TestBase { mat_models::compute_pk2cc<2>(com_mod, cep_mod, dmn, convert_to_eigen_matrix>(F), nFn, - eigen_view<2>(fN), ya_g_f, ya_g_s, ya_g_n, S_e, Dm_e, J); + eigen_view<2>(fN), active_tension, S_e, Dm_e, J); convert_to_array(S_e, S); copy_Dm(Dm_e, Dm);