Skip to content
Snippets Groups Projects
Commit e3e4c877 authored by Laura Lindner's avatar Laura Lindner
Browse files

first version of usng cellevaluations in coupled

parent 6b03866d
No related branches found
No related tags found
1 merge request!261Resolve "Fix evaluation on cells"
Pipeline #136664 passed
......@@ -128,21 +128,6 @@ void SegregatedSolver::printFullEvaluation() {
mout << "===================================================" << endl;
}
void SegregatedSolver::Initialize(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble,
Vector &potential,
Vector &displacement,
Vector& iota_c,
Vector& gamma_f_c) {
elphySolver->Initialize(elphyAssemble, potential, gamma_f_c);
mechSolver->Initialize(mechAssemble, displacement);
initEvaluation(elphyAssemble, mechAssemble);
elphySolver->SetVtk(-1);
cardiacTransfer = GetCardiacTransfer(displacement, potential);
}
void SegregatedSolver::Initialize(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble,
Vector &potential,
Vector &displacement) {
......@@ -159,14 +144,20 @@ bool SegregatedSolver::Method(IElphyAssemble &elphyAssemble, IElasticity &mechAs
Vector &potential,
Vector &displacement) {
const Meshes &M = displacement.GetMeshes();
//const Meshes &M = potential.GetMeshes();
LagrangeDiscretization disc0(M, 0, 1);
Vector gamma_f_c(disc0);
gamma_f_c = 0;
const Meshes &Mpot = potential.GetMeshes();
LagrangeDiscretization disc0pot(Mpot, 0, 1);
Vector gamma_f_c_pot(disc0pot);
gamma_f_c_pot = 0;
Vector iota_c(disc0);
iota_c = 1;
Vector iota_c_pot(disc0pot);
iota_c_pot = 1;
Initialize(elphyAssemble, mechAssemble, potential, displacement, iota_c, gamma_f_c);
Initialize(elphyAssemble, mechAssemble, potential, displacement);//, iota_c, gamma_f_c_pot);
Vectors vNew(3, potential, false);
v_potential(vNew) = potential;
......@@ -184,7 +175,7 @@ bool SegregatedSolver::Method(IElphyAssemble &elphyAssemble, IElasticity &mechAs
PlotIteration(elphyAssemble, mechAssemble, vNew, uNew);
while (!mechAssemble.IsFinished()) {
bool converged = Step(elphyAssemble, mechAssemble, vNew, uNew, iota_c, gamma_f_c);
bool converged = Step(elphyAssemble, mechAssemble, vNew, uNew, iota_c, iota_c_pot, gamma_f_c,gamma_f_c_pot);
if (converged) {
// mpp::plot("gamma_f_c" + std::to_string(mechAssemble.Time())) << gamma_f_c << mpp::endp;
// mpp::plot("iota_c" + std::to_string(mechAssemble.Time())) << iota_c << mpp::endp;
......@@ -212,7 +203,9 @@ bool SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAsse
Vectors &elphyValues,
Vectors &mechValues,
Vector& iota_c,
Vector& gamma_f_c) {
Vector& iota_c_pot,
Vector& gamma_f_c,
Vector& gamma_f_c_pot) {
/* TODO: Implement deformation of conductivity tensor.
* Is Sigma always diagonal?
* Is it then possible to assign VectorFields or Scalars instead of using a deformation gradient
......@@ -241,17 +234,13 @@ bool SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAsse
<< elphyAssemble.LastTStep()
<< "] - step " << elphyAssemble.Step() + 1 << " at time " << elphyAssemble.Time() << "."
<< endl;
elphySolver->Step(elphyAssemble, elphyValues, iota_c, gamma_f_c);
elphySolver->Step(elphyAssemble, elphyValues, iota_c, gamma_f_c);//hier eigentlich iota_c_pot und gamma_f_c_pot reinstecken
}
elphyAssemble.GetElphyProblem().PrintMinMaxJacobian(v_potential(elphyValues));
// Projects stretch from elphy to mech
cardiacTransfer->Project(u_stretch(mechValues), v_stretch(elphyValues));
//Es fehlt ne Projektion von gamma_f_c auf die Größe von gamma_c
//Für untrschiedliche Gittergrößen in displacement und potential brauchen wir noch mapping von fine to corse auf Zellebene!!!
//cardiacTransfer->Restrict(gamma_f_c,gamma_f_c_pot);
mechAssemble.UpdateStretch(u_stretch(mechValues),gamma_f_c);
// Solves one mech step
......@@ -276,8 +265,8 @@ bool SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAsse
// Interpolates iota_4 from mech to elphy
mechAssemble.GetInvariant(u_displacement(mechValues), u_iota(mechValues), iota_c);
cardiacTransfer->Prolongate(u_iota(mechValues), v_iota(elphyValues));
//iota_c->iota_c_pot ,mappen!
//cardiacTransfer->Prolongate(u_iota(mechValues), v_iota(elphyValues));
mechAssemble.PrintPointEvaluation(u_iota(mechValues), "iota4 ", evaluationPoints);
mechAssemble.PointEvaluationGamma( "gamma ", evaluationPoints);
......
......@@ -72,13 +72,13 @@ public:
void Initialize(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble, Vector &potential, Vector &displacement);
void Initialize(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble, Vector &potential, Vector &displacement, Vector& iota, Vector& gamma_f);
//void Initialize(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble, Vector &potential, Vector &displacement, Vector& iota);//, Vector& gamma_f);
/// Solves coupled equation from mechAssemble.firstTimeStep to mechAssemble.lastTimeStep
bool Method(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble, Vector &potential, Vector &displacement);
/// Calculates next timestep for the mechanics as well as for the electrophysiology.
bool Step(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble, Vectors &elphyValues, Vectors &mechValues, Vector& iota, Vector& gamma_f);
bool Step(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble, Vectors &elphyValues, Vectors &mechValues, Vector& iota_c,Vector& iota_c_pot, Vector& gamma_f_c,Vector& gamma_f_c_pot);
bool Step(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble, Vectors &elphyValues, Vectors &mechValues);
......
......@@ -238,7 +238,7 @@ void LinearImplicitSolver::SolveConcentration(double t, double dt, const Vectors
}
}
void LinearImplicitSolver::SolveConcentrationOnCells(IElphyAssemble &A, const Vectors &values, Vector &iota_c, Vector &gamma_f_c) {
void LinearImplicitSolver::SolveConcentrationOnCells(IElphyAssemble &A, const Vectors &values, Vector &iota_c_pot) {
double t = A.Time();
double dt = A.StepSize();
......@@ -257,7 +257,7 @@ void LinearImplicitSolver::SolveConcentrationOnCells(IElphyAssemble &A, const Ve
vcw[i] = (*gating_c)[i](r, j);
}
if (M.Type() == TENSION) {
vcw[M.Size()] = iota_c(r, j);
vcw[M.Size()] = iota_c_pot(r, j);
}
SolveConcentrationStep(M, t, M.TimeScale() * dt, vcw, G);
for (int i: ConcentrationIndices(M)) {
......@@ -309,9 +309,9 @@ void LinearImplicitSolver::selectOrder(std::string o) {
}
};
void LinearImplicitSolver::StaggeredStepOnCells(IElphyAssemble &A, const Vectors &values,Vector &iota_c, Vector &gamma_f_c){
void LinearImplicitSolver::StaggeredStepOnCells(IElphyAssemble &A, const Vectors &values,Vector &iota_c){
SolveGatingOnCells(A.Time(), A.StepSize());
SolveConcentrationOnCells(A,values, iota_c, gamma_f_c);
SolveConcentrationOnCells(A,values, iota_c);
A.updateIionVecs(*gating_c);
SolvePDEOnCells(A);
if (plotCa) A.PlotCa((*gating_c)[1]);
......@@ -387,18 +387,12 @@ void SemiImplicitSolverOnCells::updateValues(Vectors &values) {
}
}
void SemiImplicitSolverOnCells::updateValues(Vectors &values, Vector &gamma_f_c) {
void SemiImplicitSolverOnCells::updateValues(Vectors &values, Vector &gamma_f_c_pot) {
values[0] = *potential;
if (M.Type() == TENSION){
//ProjectFromfinerToCoarseCells()
/* for (row r = gamma_f_c.rows(); r != gamma_f_c.rows_end(); ++r) {
if ((*gating_c).find_row(r()) == (*gating_c).rows_end()) continue;
for (int j = 0; j < r.n(); ++j) {
gamma_f_c(r(), j) = (*gating_c)[M.ElphySize()](r(), j);
}*/
for (row r = (*gating_c)[0].rows(); r != (*gating_c)[0].rows_end(); ++r) {
for (int j = 0; j < r.n(); ++j) {
gamma_f_c(r(),j) = (*gating_c)[M.ElphySize()](r(),j);
gamma_f_c_pot(r(),j) = (*gating_c)[M.ElphySize()](r(),j);
}
}
......@@ -431,7 +425,7 @@ void SemiImplicitSolverOnCells::Step(IElphyAssemble &A, Vectors &values) {
A.Initialize(values[0]);
A.updateExternalCurrent(values[0]);
SolveGatingOnCells(A.Time(), A.StepSize());
SolveConcentrationOnCells(A,values, values[0], values[0]);
SolveConcentrationOnCells(A,values, values[0]);
A.updateIionVecs(*gating_c);
SolvePDEOnCells(A);
if (plotCa) A.PlotCa((*gating_c)[1]);
......@@ -439,12 +433,12 @@ void SemiImplicitSolverOnCells::Step(IElphyAssemble &A, Vectors &values) {
A.NextTimeStep();
A.updateActivationTime(values[0], A.Time());
}
void SemiImplicitSolverOnCells::Step(IElphyAssemble &A, Vectors &values, Vector &iota_c, Vector &gamma_f_c) {
void SemiImplicitSolverOnCells::Step(IElphyAssemble &A, Vectors &values, Vector &iota_c, Vector &gamma_f_c_pot) {
A.Initialize(values[0]);
A.updateExternalCurrent(values[0]);
StaggeredStepOnCells(A, values, iota_c, gamma_f_c);
StaggeredStepOnCells(A, values, iota_c);
updateValues(values,gamma_f_c);
updateValues(values,gamma_f_c_pot);
A.NextTimeStep();
A.updateActivationTime(values[0], A.Time());
......
......@@ -95,7 +95,7 @@ public:
void SolveConcentration(double t, double dt, const Vectors &values);
void SolveConcentrationOnCells(IElphyAssemble &A, const Vectors &values, Vector &iota_c, Vector &gamma_f_c);
void SolveConcentrationOnCells(IElphyAssemble &A, const Vectors &values, Vector &iota_c);
void SolveGating(double t, double dt);
void SolveGatingOnCells(double t,double dt);
......@@ -108,7 +108,7 @@ public:
void Step(IElphyAssemble &A, Vectors &values) override;
std::function<void(IElphyAssemble &A, const Vectors &values)> StaggeredStep;
void StaggeredStepOnCells(IElphyAssemble &A, const Vectors &values,Vector &iota_c, Vector &gamma_f_c);
void StaggeredStepOnCells(IElphyAssemble &A, const Vectors &values,Vector &iota_c);
void selectOrder(std::string o);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment