diff --git a/src/coupled/solvers/SegregatedSolver.cpp b/src/coupled/solvers/SegregatedSolver.cpp index 03ff16a025667d53dc2b35d76e6027d994443b2a..7a48f2c541bf1a0c165fe356b375cfe442a1b6b5 100644 --- a/src/coupled/solvers/SegregatedSolver.cpp +++ b/src/coupled/solvers/SegregatedSolver.cpp @@ -144,15 +144,14 @@ void SegregatedSolver::Method(IElphyAssemble &elphyAssemble, IElasticity &mechAs cardiacTransfer->Project(uNew[3], vNew[0]); - elphyAssemble.SetInitialValue(vNew[0]); - mechAssemble.SetInitialValue(uNew[0]); + elphyAssemble.SetInitialValue(vNew); + mechAssemble.SetInitialValue(uNew); EvaluateIteration(elphyAssemble, mechAssemble, vNew, uNew); PlotIteration(elphyAssemble, mechAssemble, vNew, uNew); while (!mechAssemble.IsFinished()) { - mechAssemble.NextTimeStep(); Step(elphyAssemble, mechAssemble, vNew, uNew); EvaluateIteration(elphyAssemble, mechAssemble, vNew, uNew); PlotIteration(elphyAssemble, mechAssemble, vNew, uNew); @@ -172,8 +171,9 @@ void SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAsse * Is it then possible to assign VectorFields or Scalars instead of using a deformation gradient * inside Elphy? */ - - elphyAssemble.ResetTime(elphyAssemble.Time(), mechAssemble.Time(), elphySteps); + vout(2) << "Solving MechStep from " << mechAssemble.Time() << " to " + << mechAssemble.NextTimeStep(false) << endl; + elphyAssemble.ResetTime(mechAssemble.Time(), mechAssemble.NextTimeStep(false), elphySteps); //elphyDeformation->Update(uNew[0]); while (!elphyAssemble.IsFinished()) { ElphyStep(elphyAssemble, elphyValues); diff --git a/src/elasticity/assemble/IElasticity.cpp b/src/elasticity/assemble/IElasticity.cpp index 68221d606a326b588bf24f12a501b677cbb6f0cc..701456c95d574a2f9c536c57b7b0693f89a7da99 100644 --- a/src/elasticity/assemble/IElasticity.cpp +++ b/src/elasticity/assemble/IElasticity.cpp @@ -36,6 +36,12 @@ void IElasticity::SetInitialValue(Vector &u) { } +void IElasticity::SetInitialValue(Vectors &values) { + SetInitialValue(values[0]); + values[1] = *gamma; + GetInvariant(values[0], values[2]); +} + void IElasticity::InitializePrestress(const Vector &u) { Prestress = std::make_unique<Vector>(0.0, u); @@ -133,10 +139,11 @@ void IElasticity::PlotIteration(const Vectors &mechValues, int step) const { PlotDeformed(mechValues[0], { {"Displacement", mechValues[0]}, - {"Stretch", mechValues[1]}, - {"Potential", mechValues[3]}, - {"Iota", mechValues[2]}}, - step, "U");} + {"Stretch", mechValues[1]}, + {"Potential", mechValues[3]}, + {"Iota", mechValues[2]}}, + step, "U"); +} std::unique_ptr<IElasticity> diff --git a/src/elasticity/assemble/IElasticity.hpp b/src/elasticity/assemble/IElasticity.hpp index 569e9280e25619b5e4cb0e9852e6a85089453c23..37b29722f90275d32850722baf2b98163c7a4625 100644 --- a/src/elasticity/assemble/IElasticity.hpp +++ b/src/elasticity/assemble/IElasticity.hpp @@ -212,6 +212,8 @@ public: void SetInitialValue(Vector &u) override; + void SetInitialValue(Vectors &values); + void PlotIteration(const Vector &u) const override; void PlotIteration(const Vectors &u, int step = -1) const; diff --git a/src/electrophysiology/assemble/IElphyAssemble.hpp b/src/electrophysiology/assemble/IElphyAssemble.hpp index f0e985e8e3f3f928d22b53e60d85a9b7232a5055..b96c57678606e0406b49821345118441f91c02c7 100644 --- a/src/electrophysiology/assemble/IElphyAssemble.hpp +++ b/src/electrophysiology/assemble/IElphyAssemble.hpp @@ -46,10 +46,12 @@ public: return std::make_unique<MultiPartDiscretization>(disc, size); } - virtual void Initialize (Vector& u) const {}; + virtual void Initialize(Vector &u) const {}; void SetInitialValue(Vector &u) override; + virtual void SetInitialValue(Vectors &values) { SetInitialValue(values[0]); }; + void FinishTimeStep(const Vector &u) override; virtual void RHS(const Vector &u, Vector &A) const = 0; diff --git a/src/electrophysiology/assemble/Monodomain.cpp b/src/electrophysiology/assemble/Monodomain.cpp index 6d1a1f5a45d0fa67c52f7440d39cd9b7d8536f68..7c96d46ae99cc4af094cc4e7b4e8ecd01bd01c49 100644 --- a/src/electrophysiology/assemble/Monodomain.cpp +++ b/src/electrophysiology/assemble/Monodomain.cpp @@ -25,11 +25,17 @@ void Monodomain::InitializeMonodomainParameters() { void Monodomain::SetInitialValue(Vector &u) { IElphyAssemble::SetInitialValue(u); - u = elphyModel.InitialValue(); + u = elphyModel.InitialValue(ELPHY); updateExternalCurrent(u); } +void Monodomain::SetInitialValue(Vectors &values) { + SetInitialValue(values[0]); + values[1] = elphyModel.InitialValue(TENSION); + +} + void Monodomain::updateExternalCurrent(const Vector &V) { if (!isNormIextZero) { double tNew = Time() + StepSize(); diff --git a/src/electrophysiology/assemble/Monodomain.hpp b/src/electrophysiology/assemble/Monodomain.hpp index 44f8b11819c16a2fbac023f92bd5a4efaea170cd..e0d099a22c688ac236f525f1aee32497d105eb22 100644 --- a/src/electrophysiology/assemble/Monodomain.hpp +++ b/src/electrophysiology/assemble/Monodomain.hpp @@ -26,6 +26,8 @@ public: void SetInitialValue(Vector &u) override; + void SetInitialValue(Vectors &values) override; + void updateExternalCurrent(const Vector &V) override; virtual void updateIionVecs(const Vectors &VCW) override {}; @@ -88,6 +90,7 @@ public: virtual void SystemMatrix(const cell &c, const Vector &u, Matrix &A) const { Exit("Not implemented for this assemble class!") }; + };