Skip to content
Snippets Groups Projects
Commit ac9759d2 authored by Laura Stengel's avatar Laura Stengel
Browse files

Merge remote-tracking branch 'origin/212-document-code' into 212-document-code

parents b08d3bd8 c37c0533
No related branches found
No related tags found
1 merge request!231Resolve "Document Code"
Pipeline #125731 passed
Stage: build
Stage: test
Hide whitespace changes
Inline Side-by-side
conf/coupled/gamm-paper/coarse-biventricle.conf
+ 1
− 0
View file @ ac9759d2
......@@ -2,6 +2,7 @@ Model = Coupled
CoupledModel = Segregated
CoupledProblem = BiventricleCoarse
#CoupledProblem = BiventricleDirichlet
ElphyModelClass = MElphyModel
ElphyModelName = TenTusscher
......
src/coupled/solvers/SegregatedSolver.cpp
+ 53
− 36
View file @ ac9759d2
......@@ -8,6 +8,15 @@
#include "solver/Newton.hpp"
#include "SegregatedSolver.hpp"
#define v_potential(v) v[0]
#define v_stretch(v) v[1]
#define v_iota(v) v[2]
#define u_displacement(u) u[0]
#define u_stretch(u) u[1]
#define u_iota(u) u[2]
#define u_pot(u) u[3]
std::vector<std::string> metrics = {
"Time",
......@@ -37,33 +46,35 @@ void
SegregatedSolver::EvaluateIteration(const IElphyAssemble &elphyAssemble, const IElasticity &mechAssemble,
Vectors elphyValues, Vectors mechValues) {
auto volumes = mechAssemble.ChamberVolume(mechValues[0]);
auto volumes = mechAssemble.ChamberVolume(u_displacement(mechValues));
evaluationMetrics[metrics[0]].emplace_back(mechAssemble.Time());
evaluationMetrics[metrics[1]].emplace_back(mechAssemble.Step());
evaluationMetrics[metrics[2]].emplace_back(mechAssemble.DeformedVolume(mechValues[0]));
evaluationMetrics[metrics[2]].emplace_back(
mechAssemble.DeformedVolume(u_displacement(mechValues)));
evaluationMetrics[metrics[3]].emplace_back(volumes[0]);
evaluationMetrics[metrics[4]].emplace_back(volumes[1]);
evaluationMetrics[metrics[5]].emplace_back(volumes[2]);
evaluationMetrics[metrics[6]].emplace_back(volumes[3]);
evaluationMetrics[metrics[7]].emplace_back(norm(mechValues[0]));
evaluationMetrics[metrics[8]].emplace_back(mechAssemble.L2Norm(mechValues[0]));
evaluationMetrics[metrics[9]].emplace_back(mechAssemble.L2AvgNorm(mechValues[0]));
evaluationMetrics[metrics[10]].emplace_back(mechAssemble.H1Norm(mechValues[0]));
evaluationMetrics[metrics[11]].emplace_back(mechAssemble.EnergyError(mechValues[0]));
evaluationMetrics[metrics[12]].emplace_back(norm(elphyValues[0]));
evaluationMetrics[metrics[13]].emplace_back(elphyAssemble.L2Norm(elphyValues[0]));
evaluationMetrics[metrics[14]].emplace_back(elphyAssemble.H1Norm(elphyValues[0]));
evaluationMetrics[metrics[15]].emplace_back(norm(elphyValues[1]));
evaluationMetrics[metrics[16]].emplace_back(elphyAssemble.L2Norm(elphyValues[1]));
evaluationMetrics[metrics[17]].emplace_back(elphyAssemble.H1Norm(elphyValues[1]));
evaluationMetrics[metrics[18]].emplace_back(norm(mechValues[1]));
evaluationMetrics[metrics[19]].emplace_back(elphyAssemble.L2Norm(mechValues[1]));
evaluationMetrics[metrics[20]].emplace_back(elphyAssemble.H1Norm(mechValues[1]));
evaluationMetrics[metrics[7]].emplace_back(norm(u_displacement(mechValues)));
evaluationMetrics[metrics[8]].emplace_back(mechAssemble.L2Norm(u_displacement(mechValues)));
evaluationMetrics[metrics[9]].emplace_back(mechAssemble.L2AvgNorm(u_displacement(mechValues)));
evaluationMetrics[metrics[10]].emplace_back(mechAssemble.H1Norm(u_displacement(mechValues)));
evaluationMetrics[metrics[11]].emplace_back(mechAssemble.EnergyError(u_displacement(mechValues)));
evaluationMetrics[metrics[12]].emplace_back(norm(v_potential(elphyValues)));
evaluationMetrics[metrics[13]].emplace_back(elphyAssemble.L2Norm(v_potential(elphyValues)));
evaluationMetrics[metrics[14]].emplace_back(elphyAssemble.H1Norm(v_potential(elphyValues)));
evaluationMetrics[metrics[15]].emplace_back(norm(v_stretch(elphyValues)));
evaluationMetrics[metrics[16]].emplace_back(elphyAssemble.L2Norm(v_stretch(elphyValues)));
evaluationMetrics[metrics[17]].emplace_back(elphyAssemble.H1Norm(v_stretch(elphyValues)));
evaluationMetrics[metrics[18]].emplace_back(norm(u_stretch(mechValues)));
evaluationMetrics[metrics[19]].emplace_back(elphyAssemble.L2Norm(u_stretch(mechValues)));
evaluationMetrics[metrics[20]].emplace_back(elphyAssemble.H1Norm(u_stretch(mechValues)));
// Get some metrices from src/elasticity/problem/ElasticityProblem.hpp (for the mechanics part)
auto problemMetrics = mechAssemble.GetElasticityProblem().EvaluationMetrics();
auto problemEntries = mechAssemble.GetElasticityProblem().EvaluationResults(mechValues[0]);
auto problemEntries = mechAssemble.GetElasticityProblem().EvaluationResults(
u_displacement(mechValues));
for (int i = 0; i < problemMetrics.size(); ++i) {
evaluationMetrics[problemMetrics[i]].emplace_back(problemEntries[i]);
}
......@@ -126,16 +137,16 @@ void SegregatedSolver::Method(IElphyAssemble &elphyAssemble, IElasticity &mechAs
Initialize(elphyAssemble, mechAssemble, potential, displacement);
Vectors vNew(3, potential, false);
vNew[0] = potential;
vNew[1] = 0.0; // Stretch
vNew[2] = 1.0; // Iota4
v_potential(vNew) = potential;
v_stretch(vNew) = 0.0;
v_iota(vNew) = 1.0;
Vectors uNew(4, displacement, false);
uNew[0] = displacement;
uNew[1] = 0.0; // Stretch
uNew[2] = 1.0; // Iota4
uNew[3] = 0.0; // Potential for plotting
cardiacTransfer->Project(uNew[3], vNew[0]);
u_displacement(uNew) = displacement;
u_stretch(uNew) = 0.0; // Stretch
u_iota(uNew) = 1.0; // Iota4
u_pot(uNew) = 0.0; // Potential for plotting
cardiacTransfer->Project(u_pot(uNew), vNew[0]);
EvaluateIteration(elphyAssemble, mechAssemble, vNew, uNew);
PlotIteration(elphyAssemble, mechAssemble, vNew, uNew);
......@@ -146,8 +157,8 @@ void SegregatedSolver::Method(IElphyAssemble &elphyAssemble, IElasticity &mechAs
PlotIteration(elphyAssemble, mechAssemble, vNew, uNew);
}
potential = vNew[0];
displacement = uNew[0];
potential = v_potential(vNew);
displacement = u_displacement(uNew);
printFullEvaluation();
}
......@@ -162,8 +173,10 @@ void SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAsse
*/
vout(2) << "Solving MechStep from " << mechAssemble.Time() << " to "
<< mechAssemble.NextTimeStep(false) << endl;
elphyAssemble.ResetTime(mechAssemble.Time(), mechAssemble.NextTimeStep(false), elphySteps);
elphyAssemble.UpdateDeformation(mechValues[0], elphyValues[0]);
elphyAssemble.UpdateDeformation(u_displacement(mechValues), v_potential(elphyValues));
while (!elphyAssemble.IsFinished()) {
vout(3) << "Solving ElphyProblem in [" << elphyAssemble.FirstTStep() << ", "
<< elphyAssemble.LastTStep()
......@@ -172,19 +185,23 @@ void SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAsse
elphySolver->SolveStep(elphyAssemble, elphyValues);
}
cardiacTransfer->Project(mechValues[1], elphyValues[1]);
mechAssemble.UpdateStretch(mechValues[1]);
// Projects stretch from elphy to mech
cardiacTransfer->Project(u_stretch(mechValues), v_stretch(elphyValues));
mechAssemble.UpdateStretch(u_stretch(mechValues));
mechSolver->Step(mechAssemble, mechValues[0]);
// Solves one mech step
mechSolver->Step(mechAssemble, u_displacement(mechValues));
mechAssemble.GetInvariant(mechValues[0], mechValues[2]);
cardiacTransfer->Prolongate(mechValues[2], elphyValues[2]);
cardiacTransfer->Project(mechValues[3], elphyValues[0]);
// Interpolates iota_4 from mech to elphy
mechAssemble.GetInvariant(u_displacement(mechValues), u_iota(mechValues));
cardiacTransfer->Prolongate(u_iota(mechValues), v_iota(elphyValues));
}
void SegregatedSolver::PlotIteration(const IElphyAssemble &elphyAssemble,
const IElasticity &mechAssemble, const Vectors &elphyValues,
const Vectors &mechValues) {
const IElasticity &mechAssemble, Vectors &elphyValues,
Vectors &mechValues) const {
cardiacTransfer->Project(u_pot(mechValues), v_potential(elphyValues));
if (vtk <= 0) return;
if (elphyVtk > 0) elphyAssemble.PlotIteration(elphyValues, mechAssemble.Step());
......
src/coupled/solvers/SegregatedSolver.hpp
+ 13
− 6
View file @ ac9759d2
......@@ -12,6 +12,15 @@
#include "solvers/ElphySolver.hpp"
#include "solvers/DynamicMechanicsSolver.hpp"
/**
* Solves coupled electro-elastodynamical problem with the following scheme:
* 0. Initialize ELPHY and MECH
* 1. Solve ELPHY from t_n to t_{n+1} in m steps
* 2. Project values from ELPHY to MECH
* 3. Solve MECH from t_n to t_{n+1} in 1 step
* 4. Interpolate values from MECH to ELPHY
* 5. If t_{n+1} = t_N STOP, else return to 1.
*/
class SegregatedSolver : public CardiacSolver {
int elphyVtk{0};
int mechVtk{0};
......@@ -25,14 +34,12 @@ class SegregatedSolver : public CardiacSolver {
/// Creates some output for each step.
void printStep();
void printFullEvaluation();
void initEvaluation(const IElphyAssemble &elphyAssemble,
const IElasticity &elasticityAssemble);
public:
/*explicit SegregatedSolver() : SegregatedSolver(GetElphySolver(), GetDynamicSolver()) {}
explicit SegregatedSolver(const string &elphySolverName, const string& mechSolverName)
: SegregatedSolver(GetElphySolver(elphySolverName), GetDynamicSolver(mechSolverName)) {}*/
explicit SegregatedSolver(
std::unique_ptr<ElphySolver> &&eSolver,
std::unique_ptr<ElastodynamicTimeIntegrator> &&mSolver,
......@@ -56,9 +63,9 @@ public:
void EvaluateIteration(const IElphyAssemble &elphyAssemble, const IElasticity &mechAssemble,
Vectors elphyValues, Vectors mechValues);
void PlotIteration(const IElphyAssemble &elphyAssemble, const IElasticity &mechAssemble,
const Vectors &elphyValues, const Vectors &mechValues);
void PlotIteration(const IElphyAssemble &elphyAssemble, const IElasticity &mechAssemble,
Vectors &elphyValues, Vectors &mechValues) const;
void Initialize(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble, Vector &potential, Vector &displacement);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Help | Imprint | Privacy policy | Accessibility | Contact