diff --git a/src/coupled/solvers/SegregatedSolver.cpp b/src/coupled/solvers/SegregatedSolver.cpp index 8782504c43546dd2a0b653d21e0f115398e20810..51ccceb2f018d462b1bb116975c5904dd36a2eb9 100644 --- a/src/coupled/solvers/SegregatedSolver.cpp +++ b/src/coupled/solvers/SegregatedSolver.cpp @@ -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); diff --git a/src/coupled/solvers/SegregatedSolver.hpp b/src/coupled/solvers/SegregatedSolver.hpp index 772206210e8b619422b8fc9e79a2487892a6de4c..ef0bbacb69a3683d18f8fc63a4a72015af7e2695 100644 --- a/src/coupled/solvers/SegregatedSolver.hpp +++ b/src/coupled/solvers/SegregatedSolver.hpp @@ -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); diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.cpp b/src/electrophysiology/solvers/LinearImplicitSolver.cpp index f8b9ddb166fbafff97cae55688de3b2fa0e80aa4..9c183da7ec5b5f2099a7fc116376d2439b67c36c 100644 --- a/src/electrophysiology/solvers/LinearImplicitSolver.cpp +++ b/src/electrophysiology/solvers/LinearImplicitSolver.cpp @@ -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()); diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.hpp b/src/electrophysiology/solvers/LinearImplicitSolver.hpp index 040ef9629877dcf236ad2f6a7760151138fd78cc..ef81917b06c9bc2ba889be4ff0dfd266b1536757 100644 --- a/src/electrophysiology/solvers/LinearImplicitSolver.hpp +++ b/src/electrophysiology/solvers/LinearImplicitSolver.hpp @@ -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);