diff --git a/conf/coupled/gamm-paper/coarse-biventricle.conf b/conf/coupled/gamm-paper/coarse-biventricle.conf
index cfbed51b6d2febb0321fb46b1d20a883635b3d0c..3cdf22a703858f85432a0460caba8f3a4432df14 100644
--- a/conf/coupled/gamm-paper/coarse-biventricle.conf
+++ b/conf/coupled/gamm-paper/coarse-biventricle.conf
@@ -1,3 +1,6 @@
+WithPrestress = false
+CoupledProblem = BiventricleCoarse
+
#VolumetricPenalty = 10000
#PrestressSteps = 20
#CoupledProblem = BiventricleCoarse
@@ -7,8 +10,8 @@ ActiveMaterial = Bonet
PassiveMaterial = Bonet
GuccioneMat_C = 0.29
GuccioneMat_bf = 42.0
-GuccioneMat_bs = 16.8
-GuccioneMat_bfs = 29.4
+GuccioneMat_bs = 16.8
+GuccioneMat_bfs = 29.4
# = Removes all files in the data folder upon start ======
ClearData = true
@@ -159,7 +162,8 @@ VolumetricPenalty = 1000 # ?
PermeabilityPenalty = 0 #1000 # ?
-ElphyModel = SemiImplicit
+#ElphyModel = SemiImplicit
+ElphyModel = SemiImplicitOnCells
OrderStaggeredScheme = wcV
ReassembleRHSonce = 0
@@ -262,13 +266,13 @@ ElphyVerbose = 1
MechVerbose = 1
PrintVSteps = 10
# = Plot Verbose =
-ElphyVTK = 1 # 0 Disables all plots from Electrophysiology
-MechVTK = 1 # 0 Disables all plots from Elasticity
+ElphyVTK = -1 # 0 Disables all plots from Electrophysiology
+MechVTK = -1 # 0 Disables all plots from Elasticity
# = Plot Verbose =
PressureSolverVTK = 0
-DynamicSolverVTK = 10
-ElphySolverVTK = 10
-CoupledSolverVTK = 10
+DynamicSolverVTK = -10
+ElphySolverVTK = -10
+CoupledSolverVTK = -10
PlotCalcium=0
-PlotVTK = 1
\ No newline at end of file
+PlotVTK = -1
\ No newline at end of file
diff --git a/conf/electrophysiology/BiVentricleTest/test.conf b/conf/electrophysiology/BiVentricleTest/test.conf
index 650df2d9c07951e0d872283d9557a0d212210ce2..39872646c9ea2b482f6d4dbc1f847bcdfe704f53 100644
--- a/conf/electrophysiology/BiVentricleTest/test.conf
+++ b/conf/electrophysiology/BiVentricleTest/test.conf
@@ -2,7 +2,10 @@ ClearData = 1 # Removes existing plots from data folder.
Distribution=RCB
Model = Monodomain;
ElphyModel = SemiImplicit
+#mit activation times
ElphyProblem =BiVentricleWAProblem
+# ohne activation times
+#ElphyProblem= BiVentricleProblem
# Spielt für diesen Test keine Rolle, nur beim Splitting (Iext splitted in PDE =1, Iext in ODE =0)
ElphySplittingMethod =Godunov
IextInPDE=1
diff --git a/conf/electrophysiology/cuboidTests.conf b/conf/electrophysiology/cuboidTests.conf
index 30f7c1d5577f8294845c74464dce48ee279b9625..7b10e1f10c83376f129973285953e35d25957c1c 100644
--- a/conf/electrophysiology/cuboidTests.conf
+++ b/conf/electrophysiology/cuboidTests.conf
@@ -7,10 +7,12 @@ Model = Monodomain;
#ElphyModel = Diffusion
#ElphyModel = Splitting
+#ElphySplittingMethod =Godunov
#ElphyModel = ImplictEuler
-ElphyModel = SemiImplicit
-#ElphyModel = SemiImplicitNodes
-#ElphyModel = LinearImplicit
+#ElphyModel = SemiImplicit
+ElphyModel = SemiImplicitOnCells
+ElphyModel = SemiImplicitNodes
+ElphyModel = LinearImplicit
#ElphyModel = LinearImplicitNodes
#ElphyModel = LinearImplicitQuadrature
OrderStaggeredScheme=wcV
@@ -20,11 +22,13 @@ ReassembleRHSonce=1
ElphyProblem = ElphyBenchmarkNiederer
-#ElphyProblem=VariableMeshProblem
+
ElphyModelClass = MElphyModel
ElphyModelName = TenTusscher
+#ElphyModelName=BeelerReuter
+#IonScheme = ExplicitEuler
CellModelScheme = ExponentialIntegrator
IonScheme = Quadratic
ThresholdVForActivation=0.0
@@ -33,7 +37,7 @@ ThresholdVForActivation=0.0
IextInPDE=1
#Glaettung von Anregung in Zeit: Stepfunction,Arctan,Tanh
-ExternalCurrentSmoothingInTime = Arctan;
+ExternalCurrentSmoothingInTime = Stepfunction;
#Glättung Anregung im Raum: discrete, linear
ExternalCurrentSmoothingInSpace =linear;
@@ -42,13 +46,14 @@ ExternalCurrentSmoothingInSpace =linear;
######################
# Meshes #
######################
+#Mesh = TestBiventricle
#Mesh=NiedererBM05_3DQuad
ElphyPolynomialDegree = 1
#HexaQuadExactupTo =6
StartTime = 0.0
-EndTime=0.0001
-EndTime = 0.06
+EndTime=0.0002
+EndTime = 0.008
EndTime = 0.4
DeltaTime = 0.0001
@@ -113,7 +118,7 @@ GatingReduction = 1e-12
# === Elphy Solver ================================================
ElphySolver = gmres
ElphyPreconditioner=Jacobi
-ElphyEpsilon = 1e-8
+ElphyEpsilon = 1e-12
ElphyReduction = 1e-12
diff --git a/conf/electrophysiology/m++.conf b/conf/electrophysiology/m++.conf
index 2dfcc29dfac146a098a42df2a3f74af2f6a6dd06..92557e858c3fc95b1590a4ecc5112cb6a2d14018 100644
--- a/conf/electrophysiology/m++.conf
+++ b/conf/electrophysiology/m++.conf
@@ -4,16 +4,17 @@ Executable = Elphy
ClearData = true
#loadconf = electrophysiology/niedererbenchmark.conf
loadconf=electrophysiology/cuboidTests.conf
+#loadconf=electrophysiology/BiVentricleTest/biventricelTest.conf
#loadconf = electrophysiology/biventricelTest.conf
# ======================================
# = Choose Configuration for Problems =
# ======================================
# === Electrophysiology ================
-loadconf = electrophysiology/elphy.conf
+#loadconf = electrophysiology/elphy.conf
#loadconf = electrophysiology/monodomain.conf;
# ======================================
# = Default Solver Configuration Files =
-loadconf = electrophysiology/solvers.conf;
-loadconf = electrophysiology/cellmodels.conf
-loadconf = electrophysiology/verbose.conf
+#loadconf = electrophysiology/solvers.conf;
+#loadconf = electrophysiology/cellmodels.conf
+#loadconf = electrophysiology/verbose.conf
# ======================================
diff --git a/conf/m++.conf b/conf/m++.conf
index 334d494eb17d220ec5d1c298265b0bd1c02c20a4..2e658adb4cadf059cc08d97e48942b12d40a5600 100644
--- a/conf/m++.conf
+++ b/conf/m++.conf
@@ -1,5 +1,5 @@
-#loadconf = electrophysiology/m++.conf;
-loadconf = elasticity/m++.conf;
+loadconf = electrophysiology/m++.conf;
+#loadconf = elasticity/m++.conf;
#loadconf = coupled/m++.conf;
#loadconf = coupled/gamm-paper/coarse-biventricle.conf
\ No newline at end of file
diff --git a/mpp b/mpp
index 6f97c4f4ce3e198686af72f98524349285e21fc7..20c09dfd1f52ee0593c136ced501b1d54c78d524 160000
--- a/mpp
+++ b/mpp
@@ -1 +1 @@
-Subproject commit 6f97c4f4ce3e198686af72f98524349285e21fc7
+Subproject commit 20c09dfd1f52ee0593c136ced501b1d54c78d524
diff --git a/src/coupled/solvers/SegregatedSolver.cpp b/src/coupled/solvers/SegregatedSolver.cpp
index a158c4612dab71573180eec671cb08d44e15b493..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) {
@@ -162,10 +147,17 @@ bool SegregatedSolver::Method(IElphyAssemble &elphyAssemble, IElasticity &mechAs
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;
@@ -183,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;
@@ -211,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
@@ -228,22 +222,25 @@ bool SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAsse
elphyAssemble.UpdateDeformation(u_displacement(mechValues), dyn[0], v_potential(elphyValues));
elphyAssemble.PrintV(v_potential(elphyValues));
- Vector ca(elphyValues[0]);
- elphySolver->GatingVectorAti(ca,1);
- elphyAssemble.PrintVariable(ca,"Ca");
+ if(plotCa){
+ Vector ca(elphyValues[0]);
+ elphySolver->GatingVectorAti(ca,1);
+ elphyAssemble.PrintVariable(ca,"Ca");
+ }
+
while (!elphyAssemble.IsFinished()) {
vout(3) << "Solving ElphyProblem in [" << elphyAssemble.FirstTStep() << ", "
<< 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));
+ //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
@@ -268,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);
@@ -295,9 +292,11 @@ bool SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAsse
elphyAssemble.UpdateDeformation(u_displacement(mechValues), dyn[0], v_potential(elphyValues));
elphyAssemble.PrintV(v_potential(elphyValues));
- Vector ca(elphyValues[0]);
- elphySolver->GatingVectorAti(ca, 1);
- elphyAssemble.PrintVariable(ca, "Ca");
+ if(plotCa){
+ Vector ca(elphyValues[0]);
+ elphySolver->GatingVectorAti(ca,1);
+ elphyAssemble.PrintVariable(ca,"Ca");
+ }
while (!elphyAssemble.IsFinished()) {
vout(3) << "Solving ElphyProblem in [" << elphyAssemble.FirstTStep() << ", "
diff --git a/src/coupled/solvers/SegregatedSolver.hpp b/src/coupled/solvers/SegregatedSolver.hpp
index 2f445bb2a3385947978d5c3ee995fcc53ee9b133..ef0bbacb69a3683d18f8fc63a4a72015af7e2695 100644
--- a/src/coupled/solvers/SegregatedSolver.hpp
+++ b/src/coupled/solvers/SegregatedSolver.hpp
@@ -25,6 +25,7 @@ class SegregatedSolver : public CardiacSolver {
int elphyVtk{0};
int mechVtk{0};
int elphySteps;
+ bool plotCa{false};
std::unordered_map<std::string, std::vector<double>> evaluationMetrics{};
@@ -71,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/elasticity/problem/EllipsoidPoints.cpp b/src/elasticity/problem/EllipsoidPoints.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..30c353349abe2f8b5b231bfb0df6ce39adb7d016
--- /dev/null
+++ b/src/elasticity/problem/EllipsoidPoints.cpp
@@ -0,0 +1,843 @@
+
+#include "EllipsoidPoints.hpp"
+
+
+const std::vector<Point> ellipsoidEndocard = {
+ Point(3.764178, 5.744528, 3.286529), Point(4.171387, 5.230753, 5.0),
+ Point(5.066103, 4.599616, 3.584207), Point(-4.171387, -5.230753, 5.0),
+ Point(-5.066103, -4.599616, 3.584207), Point(-3.764178, -5.744528, 3.286529),
+ Point(-4.918419, 4.45966, -5.387152), Point(-5.778401, 3.471484, -4.581496),
+ Point(-5.031684, 4.598732, -3.865793), Point(5.480983, 2.575521, -8.52616),
+ Point(5.208109, 3.488358, -7.566456), Point(4.785261, 3.580168, -8.851359),
+ Point(0.0, -6.832549, -3.696116), Point(1.298003, -6.650078, -4.270177),
+ Point(1.360458, -6.757807, -2.955677), Point(1.362385, 6.43437, -5.819858),
+ Point(2.021953, 6.081706, -6.836571), Point(2.669063, 5.991028, -5.940665),
+ Point(-6.672435, 0.0, -5.139517), Point(-6.614786, 1.288465, -4.597935),
+ Point(-6.429871, 1.085941, -6.180935), Point(1.086903, -5.012064, -11.570264),
+ Point(0.0, -5.405972, -10.799763), Point(1.19845, -5.473064, -10.191118),
+ Point(-5.937444, 2.315824, -7.031846), Point(-6.183089, 1.113262, -7.497194),
+ Point(0.0, 6.946014, 2.107258), Point(0.0, 6.845144, 3.556014),
+ Point(-1.228985, 6.803006, 2.670005), Point(-6.505261, 2.35533, -2.586167),
+ Point(-6.864095, 1.121411, -1.922462), Point(-6.783435, 1.16804, -3.091733),
+ Point(0.0, -6.994789, 0.655811), Point(0.0, -6.946014, 2.107258),
+ Point(1.264143, -6.869007, 1.135684), Point(-3.597701, 5.061102, -7.847628),
+ Point(-4.38254, 4.693776, -6.766007), Point(-3.000892, 5.682453, -6.740929),
+ Point(-1.292619, 6.863583, 1.140043), Point(0.0, -6.992313, -0.796456),
+ Point(0.0, -6.938543, -2.247727), Point(1.396412, -6.829523, -1.550571),
+ Point(-6.992313, 0.0, -0.796456), Point(-6.994789, 0.0, 0.655811),
+ Point(-6.915463, -1.084514, -0.033912), Point(-1.1619, -6.902212, -0.236159),
+ Point(-5.639774, -3.690821, -4.589154), Point(-5.549712, -4.056619, -3.207687),
+ Point(-6.276254, -2.735845, -3.539206), Point(-3.006887, 5.906251, -5.470939),
+ Point(-3.885439, 5.523773, -4.472327), Point(-3.897481, 5.306937, -5.770636),
+ Point(-3.994127, -4.385334, -9.026905), Point(-4.928345, -3.496317, -8.581893),
+ Point(-4.191071, -4.615092, -7.731806), Point(-6.907732, 1.123651, -0.348976),
+ Point(-6.544469, 2.435699, -1.183025), Point(-6.224226, 2.303312, -5.405285),
+ Point(-6.393261, 2.320781, -4.020138), Point(1.273326, -6.88318, -0.052855),
+ Point(-1.022907, 6.7379, 3.881439), Point(0.0, -6.672435, -5.139517),
+ Point(1.294508, -6.469965, -5.676732), Point(2.689765, -6.286326, -3.640649),
+ Point(6.915463, 1.084514, -0.033912), Point(6.992313, 0.0, -0.796456),
+ Point(6.906038, 1.00325, -1.330468), Point(-4.590256, 5.157856, 2.796621),
+ Point(-3.475798, 5.990087, 2.473905), Point(-3.562444, 5.815385, 3.832447),
+ Point(-6.907251, 1.071535, 0.914166), Point(2.580489, -6.49884, -0.791249),
+ Point(2.444718, -6.556273, 0.477359), Point(-6.257566, -1.156146, -7.083023),
+ Point(-6.4731, -1.353275, -5.573885), Point(-6.038944, -2.468604, -6.16166),
+ Point(-1.141034, -6.881578, 1.420093), Point(5.49659, 1.241092, -10.085792),
+ Point(5.142728, 1.102389, -11.218077), Point(4.775967, 2.255788, -11.15599),
+ Point(3.631444, 5.554046, -5.411514), Point(2.608928, 6.21096, -4.619056),
+ Point(-6.633346, 2.229998, 0.390246), Point(-3.684563, 4.632347, -9.075634),
+ Point(-4.533108, 4.171412, -8.07304), Point(-4.634236, 3.543856, -9.394819),
+ Point(6.994789, 0.0, 0.655811), Point(0.0, 6.672435, -5.139517),
+ Point(-1.628852, 6.369241, -5.838482), Point(-1.168711, 6.638295, -4.586818),
+ Point(-2.413235, 6.51683, 2.042338), Point(-2.613759, 6.490791, 0.472807),
+ Point(-5.119003, 4.675925, -2.343952), Point(-5.906287, 3.509931, -3.25479),
+ Point(5.587209, -3.969419, 3.457458), Point(4.69061, -4.920564, 4.053772),
+ Point(5.230753, -4.171387, 5.0), Point(4.779495, 5.096302, 1.042119),
+ Point(3.800155, 5.851847, 1.362401), Point(4.432044, 5.418168, 0.050855),
+ Point(-4.432044, -5.418168, 0.050854), Point(-3.800155, -5.851847, 1.362401),
+ Point(-4.779495, -5.096302, 1.042118), Point(5.789065, 2.4836, -7.413491),
+ Point(3.244772, -6.175742, 1.398625), Point(3.145492, -6.166969, 2.517258),
+ Point(2.256435, -6.59013, 1.680252), Point(6.176336, 0.0, -8.000614),
+ Point(6.257569, 1.15614, -7.083009), Point(5.855623, 1.368803, -8.70163),
+ Point(5.94349, 3.559017, 2.438656), Point(5.785114, 3.596067, 3.916523),
+ Point(-1.488751, 6.522643, 5.0), Point(-2.346662, 6.445606, 3.389004),
+ Point(5.305606, 4.56608, -0.092688), Point(4.517225, 5.330002, -1.0467),
+ Point(-5.305606, -4.56608, -0.092689), Point(-4.517225, -5.330002, -1.046701),
+ Point(-6.740215, -1.482316, -2.844893), Point(-6.625694, -1.402576, -4.298609),
+ Point(5.468692, 4.333318, 1.364678), Point(4.881653, 4.922069, 2.35797),
+ Point(-5.468692, -4.333318, 1.364677), Point(-5.94349, -3.559017, 2.438656),
+ Point(-4.881653, -4.922069, 2.35797), Point(6.400101, 2.417633, 3.597039),
+ Point(6.519726, 2.413664, 1.984015), Point(6.803821, 1.182255, 2.779933),
+ Point(-6.519726, -2.413664, 1.984015), Point(-6.803821, -1.182255, 2.779933),
+ Point(-6.400101, -2.417633, 3.597039), Point(6.276254, 2.735845, -3.539206),
+ Point(6.740215, 1.482316, -2.844893), Point(6.62569, 1.4026, -4.298602),
+ Point(5.179294, 2.473864, -9.730943), Point(6.007858, 2.494804, -6.27756),
+ Point(6.4731, 1.353275, -5.573885), Point(1.042871, 6.227708, -7.337371),
+ Point(-6.938543, 0.0, -2.247727), Point(-6.906038, -1.00325, -1.330468),
+ Point(6.908147, -1.089212, -0.7331), Point(-2.60984, 5.548196, -8.202048),
+ Point(-1.225997, 5.90257, -8.639912), Point(-1.428351, 6.186775, -7.156512),
+ Point(-6.832549, 0.0, -3.696116), Point(-4.515363, -3.39208, -10.044218),
+ Point(-3.301864, -4.60484, -9.981738), Point(-5.871594, 1.435111, -8.574213),
+ Point(-6.176336, 0.0, -8.000614), Point(-5.82959, 0.0, -9.410796),
+ Point(-5.405972, 0.0, -10.799763), Point(-5.581368, 1.119722, -9.893131),
+ Point(3.550644, -5.817544, 3.877442), Point(4.515801, -5.223903, 2.788505),
+ Point(-3.82624, -5.641793, -3.862919), Point(-3.954992, -5.708391, -2.13426),
+ Point(-2.474514, -6.425428, -3.062946), Point(3.954992, 5.708391, -2.13426),
+ Point(2.474514, 6.425428, -3.062946), Point(3.82624, 5.641793, -3.862919),
+ Point(-2.902849, 6.027828, 5.0), Point(-1.349922, -6.833722, -1.678967),
+ Point(-0.984866, -6.82018, -2.989397), Point(1.099974, -6.720175, 3.937961),
+ Point(2.305188, -6.420603, 3.810587), Point(1.488751, -6.522643, 5.0),
+ Point(5.406606, 3.597832, -6.344204), Point(2.087813, 1.129294, -15.99279),
+ Point(2.571543, 0.0, -15.811318), Point(2.918265, 1.379487, -15.084683),
+ Point(-3.431596, 6.084651, 1.089268), Point(-1.133932, -1.981541, -16.070408),
+ Point(-2.139622, -1.089954, -15.968484), Point(-2.048484, -2.398553, -15.176271),
+ Point(-6.027828, 2.902849, 5.0), Point(-6.411347, 2.383643, 3.612678),
+ Point(-5.587544, 3.952833, 3.564672), Point(1.897214, 2.377727, -15.310981),
+ Point(1.046093, 1.815428, -16.220587), Point(1.453357, -6.758674, 2.669302),
+ Point(5.567519, 4.021618, -3.284709), Point(4.761035, 4.978178, -3.023825),
+ Point(5.315178, 4.497747, -1.74984), Point(-6.946014, 0.0, 2.107258),
+ Point(-6.688405, 1.890928, 2.016785), Point(3.189735, 5.528219, -6.981616),
+ Point(-3.54712, -6.031735, -0.461515), Point(-2.359838, -6.581688, 0.814753),
+ Point(3.54712, 6.031735, -0.461515), Point(2.359838, 6.581688, 0.814753),
+ Point(-6.455177, 0.0, -6.575415), Point(-2.533063, -6.224163, -4.761103),
+ Point(-3.685665, -5.51783, -5.413822), Point(-2.586452, -6.024279, -5.957841),
+ Point(6.581422, -2.31542, -1.381745), Point(6.224547, -3.062081, -2.276482),
+ Point(6.019061, -3.551656, -0.961178), Point(6.672435, 0.0, -5.139517),
+ Point(6.455177, 0.0, -6.575415), Point(0.0, -2.571543, -15.811318),
+ Point(0.0, -1.398063, -16.65749), Point(-5.382307, 4.412985, 1.811414),
+ Point(2.361585, -5.756431, -7.789025), Point(2.357461, -5.352553, -9.340551),
+ Point(3.428494, -5.016746, -8.440009), Point(6.729237, 1.062897, 3.906615),
+ Point(-6.729237, -1.062897, 3.906615), Point(4.766464, 4.768086, -4.573313),
+ Point(-2.257959, 4.756075, -11.203421), Point(-1.36684, 4.451445, -12.692949),
+ Point(-1.062504, 5.009573, -11.589977), Point(5.167156, -4.402744, -4.147348),
+ Point(3.777767, -5.527064, -4.965246), Point(4.780363, -4.502621, -5.886234),
+ Point(0.0, -6.845144, 3.556014), Point(4.51624, 1.549567, -12.431443),
+ Point(3.52251, 2.65472, -13.200439), Point(3.861092, 1.253725, -13.849276),
+ Point(2.902849, -6.027828, 5.0), Point(4.211816, -5.554322, 1.555274),
+ Point(4.875184, -5.020319, 0.413436), Point(3.772953, -5.896094, 0.072975),
+ Point(0.0, 5.82959, -9.410796), Point(0.0, 6.176336, -8.000614), Point(6.845144, 0.0, 3.556014),
+ Point(-6.845144, 0.0, 3.556014), Point(-3.084371, 6.148811, -3.146649),
+ Point(-1.265246, 6.767318, -3.074387), Point(-2.36897, 6.328359, -4.438276),
+ Point(6.028063, -3.076393, -4.34318), Point(6.610712, -1.512957, -4.213026),
+ Point(6.390982, -2.5094, -3.310575), Point(-4.041951, 2.756915, -12.157917),
+ Point(-4.539342, 2.930085, -10.808911), Point(-5.061655, 1.766486, -10.931048),
+ Point(2.671025, 6.274763, 3.834513), Point(1.266914, 6.704184, 3.80042),
+ Point(1.488751, 6.522643, 5.0), Point(-1.266914, -6.704184, 3.80042),
+ Point(-2.671025, -6.274763, 3.834513), Point(-1.488751, -6.522643, 5.0),
+ Point(0.0, 2.571543, -15.811318), Point(0.0, 6.455177, -6.575415),
+ Point(6.123598, -3.30198, 1.87878), Point(6.420779, -2.416811, 3.376166),
+ Point(0.0, 6.994789, 0.655811), Point(0.0, 6.992313, -0.796456),
+ Point(1.1619, 6.902212, -0.236159), Point(0.914819, 6.84803, 2.734268),
+ Point(2.360619, 6.516034, 2.390351), Point(-2.360619, -6.516034, 2.390351),
+ Point(-0.914819, -6.84803, 2.734268), Point(5.664101, 3.630364, -4.695638),
+ Point(-4.242077, 5.414186, -3.158422), Point(-6.780302, 0.968207, 3.511015),
+ Point(6.522643, 1.488751, 5.0), Point(-6.522643, -1.488751, 5.0),
+ Point(5.875652, -3.553867, -3.300309), Point(-1.673696, -6.122539, -7.168713),
+ Point(-2.876779, -5.329056, -8.526085), Point(-3.299884, -5.474326, -6.929952),
+ Point(-4.801573, -4.727805, -4.603249), Point(-4.529418, -4.677384, -6.242166),
+ Point(4.335061, 4.523093, -7.582773), Point(-6.070812, 3.481758, -0.365079),
+ Point(-4.997028, -3.92914, -7.118478), Point(0.964632, 0.892293, -16.697784),
+ Point(-5.230753, 4.171387, 5.0), Point(-4.708676, 4.921503, 3.921352),
+ Point(5.776935, -3.948053, 0.48514), Point(6.537817, -2.500334, 0.176435),
+ Point(-5.993219, 3.529543, -1.917895), Point(-3.106151, -1.398539, -14.851248),
+ Point(-3.949241, -0.96255, -13.840086), Point(-3.677936, -2.118345, -13.51853),
+ Point(6.946014, 0.0, 2.107258), Point(6.720688, -1.777257, 1.993361),
+ Point(6.895932, -1.151241, 0.843981), Point(2.068959, 5.840998, -7.907658),
+ Point(0.960975, -1.56304, -16.405613), Point(0.0, 0.0, -17.0),
+ Point(-5.347647, -2.405579, -9.284598), Point(4.368983, 4.826731, -6.246039),
+ Point(6.779151, -0.970394, 3.520552), Point(3.972389, 4.361082, -9.151908),
+ Point(-2.577451, 2.277003, -14.80672), Point(-3.291302, 2.727858, -13.46185),
+ Point(-2.138301, 3.533161, -13.726144), Point(3.032666, 4.829063, -9.859847),
+ Point(3.131431, 5.177879, -8.546264), Point(1.984201, 5.554544, -9.154775),
+ Point(6.092253, -2.56799, -5.585653), Point(5.520836, -3.684754, -5.399389),
+ Point(5.433858, -3.503363, -6.516372), Point(-1.188999, 3.083493, -14.986145),
+ Point(-0.994656, 3.791981, -14.083943), Point(5.392747, -4.369197, -2.210393),
+ Point(4.148639, -5.495195, -3.063884), Point(6.938543, 0.0, -2.247727),
+ Point(-1.034619, -3.05057, -15.093056), Point(-4.761035, -4.978178, -3.023825),
+ Point(-6.86575, -1.231298, 1.427254), Point(5.147379, 4.170317, -5.491441),
+ Point(0.0, 6.832549, -3.696116), Point(-5.689185, -2.677738, -7.470725),
+ Point(-5.462658, -3.640614, -5.901624), Point(0.0, 6.690385, 5.0),
+ Point(-5.186932, 4.684617, -0.940917), Point(-1.236144, -6.676613, -4.131987),
+ Point(1.141034, 6.881578, 1.420093), Point(-6.65833, -2.117126, -1.042925),
+ Point(6.65833, 2.117126, -1.042925), Point(-2.571543, 0.0, -15.811318),
+ Point(4.2245, 3.623008, -10.311399), Point(-1.371061, 2.05757, -15.904195),
+ Point(2.458913, -6.022141, -6.280363), Point(2.393448, -6.239713, -5.05784),
+ Point(0.0, -6.455177, -6.575415), Point(1.261749, -6.227377, -7.133456),
+ Point(-3.944618, 5.738391, -1.735807), Point(2.328664, -4.804379, -10.994556),
+ Point(-1.009156, -0.9582, -16.660683), Point(-5.785114, -3.596067, 3.916523),
+ Point(-6.027828, -2.902849, 5.0), Point(-4.070632, 5.69384, -0.244516),
+ Point(0.0, -4.892602, -12.158012), Point(5.320205, -4.489438, 1.785234),
+ Point(3.060519, -3.976127, -11.853726), Point(3.875497, -2.67401, -12.579484),
+ Point(4.385194, -3.071525, -10.951701), Point(-5.729521, 3.23511, -5.801519),
+ Point(6.690385, 0.0, 5.0), Point(6.832549, 0.0, -3.696116),
+ Point(2.600365, -3.480374, -13.329539), Point(-6.133893, 3.344814, 1.051806),
+ Point(-4.270536, 0.0, -13.469821), Point(-3.511188, 0.0, -14.706712),
+ Point(-3.595697, 1.253851, -14.264384), Point(2.715131, 3.840069, -12.591607),
+ Point(1.103935, 4.106783, -13.503306), Point(1.351252, 4.806118, -11.916182),
+ Point(-2.626274, -2.778315, -14.240554), Point(-1.437571, -3.655417, -14.071326),
+ Point(-2.979576, -3.421968, -12.945053), Point(3.603172, -5.311041, -6.786976),
+ Point(3.511188, 0.0, -14.706712), Point(3.182337, -1.291115, -14.81344),
+ Point(-6.166453, 3.166468, 2.365007), Point(5.82959, 0.0, -9.410796),
+ Point(5.405972, 0.0, -10.799763), Point(5.766663, -1.052471, -9.291597),
+ Point(-5.251275, 2.373881, -9.649958), Point(5.989247, -2.246686, -6.90386),
+ Point(-4.892602, 0.0, -12.158012), Point(-4.456588, -1.011425, -12.877349),
+ Point(-4.047475, 3.883892, -10.169113), Point(-3.592135, 3.849789, -11.201927),
+ Point(-5.189688, 4.693922, 0.449352), Point(-5.315178, -4.497746, -1.749844),
+ Point(-3.16715, 4.596134, -10.259015), Point(-2.367961, 5.271214, -9.594261),
+ Point(-1.190008, 5.475887, -10.188006), Point(2.860538, 4.411988, -11.221915),
+ Point(3.829986, 3.227555, -11.876209), Point(-4.384337, 5.426899, 1.38739),
+ Point(2.347752, 3.206014, -13.995302), Point(1.135316, 3.271845, -14.773639),
+ Point(6.057824, 3.485988, 0.942654), Point(6.597594, 2.323476, 0.657094),
+ Point(-6.057824, -3.485988, 0.942654), Point(-6.597594, -2.323476, 0.657094),
+ Point(-2.733616, 3.945426, -12.374032), Point(-2.751567, -6.399464, -1.675083),
+ Point(2.751567, 6.399464, -1.675083), Point(2.672688, -6.403569, -2.240409),
+ Point(3.617455, -5.961342, -1.490002), Point(0.0, 1.398063, -16.65749),
+ Point(4.600141, -4.296746, -7.436614), Point(5.313399, -3.278468, -7.687278),
+ Point(4.528117, -3.786813, -9.137434), Point(6.161309, -1.160715, -7.560236),
+ Point(-6.690385, 0.0, 5.0), Point(0.0, -5.82959, -9.410796),
+ Point(-6.166127, -2.617645, -4.933534), Point(6.86575, 1.231298, 1.427254),
+ Point(2.209799, 6.635191, -0.732738), Point(-2.209799, -6.635191, -0.732738),
+ Point(1.540189, -4.469592, -12.537303), Point(0.0, -4.270536, -13.469821),
+ Point(-1.171961, 1.010167, -16.579526), Point(-1.398063, 0.0, -16.65749),
+ Point(-2.508088, 1.082042, -15.652266), Point(1.238723, -5.85645, -8.81259),
+ Point(-4.171387, 5.230753, 5.0), Point(-1.441068, -4.966892, -11.456399),
+ Point(-1.454292, -4.417468, -12.705632), Point(-2.454433, -4.313584, -11.988573),
+ Point(4.171387, -5.230753, 5.0), Point(-4.299747, 1.410239, -12.970344),
+ Point(6.427207, -1.246669, -6.016207), Point(-1.220659, -5.961269, -8.40344),
+ Point(-1.64749, -5.430728, -9.952099), Point(6.027828, -2.902849, 5.0),
+ Point(-1.394605, -6.468552, -5.544867), Point(0.0, -6.176336, -8.000614),
+ Point(-4.927717, -2.263377, -10.750298), Point(-5.025668, -1.062983, -11.548547),
+ Point(-5.522375, -1.246777, -9.998204), Point(-5.301568, 3.021948, -8.328692),
+ Point(2.112528, -1.534089, -15.773342), Point(1.398063, 0.0, -16.65749),
+ Point(6.107362, 3.411051, -0.618536), Point(-6.107362, -3.411051, -0.618536),
+ Point(-2.945221, 6.342385, -0.767313), Point(-2.578456, -4.702452, -10.92547),
+ Point(0.0, -6.690385, 5.0), Point(6.229261, 3.058719, -2.226545),
+ Point(1.057366, 5.91407, -8.72454), Point(6.148955, 2.609965, -5.081713),
+ Point(6.750812, -1.518645, -2.570455), Point(-6.229261, -3.058719, -2.226545),
+ Point(1.349922, 6.833722, -1.678967), Point(3.10415, 2.248619, -14.224871),
+ Point(-5.925762, -1.252199, -8.52333), Point(6.522643, -1.488751, 5.0),
+ Point(0.0, 6.938543, -2.247727), Point(0.984866, 6.82018, -2.989396),
+ Point(-5.232724, 3.659023, -6.967144), Point(1.271831, -3.854358, -13.850603),
+ Point(-4.868447, 1.244544, -11.835227), Point(4.13653, -1.154831, -13.424419),
+ Point(4.568106, -5.28243, -1.160473), Point(0.0, -3.511188, -14.706712),
+ Point(0.0, 5.405972, -10.799763), Point(3.416043, -2.23598, -13.808952),
+ Point(2.424344, -2.639003, -14.603414), Point(1.207661, -2.958316, -15.125521),
+ Point(5.598657, -2.35545, -8.450208), Point(1.181101, 6.662658, -4.353845),
+ Point(5.235518, -2.3586, -9.722284), Point(4.270536, 0.0, -13.469821),
+ Point(0.0, 4.270536, -13.469821), Point(5.374799, -4.475919, -0.676876),
+ Point(-2.129133, 6.620033, -1.945822), Point(-1.311519, 6.872717, -0.519053),
+ Point(1.513037, 5.392173, -10.198609), Point(6.027828, 2.902849, 5.0),
+ Point(-2.902849, -6.027828, 5.0), Point(2.902849, 6.027828, 5.0),
+ Point(3.543803, -4.167062, -10.607341), Point(-4.326855, -2.252542, -12.192391),
+ Point(4.653503, -1.604171, -12.087247), Point(5.230753, 4.171387, 5.0),
+ Point(-5.230753, -4.171387, 5.0), Point(-3.697849, -3.67244, -11.349276),
+ Point(0.0, 3.511188, -14.706712), Point(3.308345, -4.701115, -9.700423),
+ Point(5.199357, -1.400357, -10.862445), Point(-6.522643, 1.488751, 5.0),
+ Point(0.0, 4.892602, -12.158012), Point(4.892602, 0.0, -12.158012)
+};
+
+const std::vector<Point> ellipsoidEpicard{
+ Point(4.928, 8.473734, 3.955074), Point(6.340662, 7.317514, 5.0),
+ Point(6.988606, 6.879922, 3.912063), Point(-6.340662, -7.317514, 5.0),
+ Point(-6.988606, -6.879922, 3.912063), Point(-4.928, -8.473734, 3.955074),
+ Point(-6.832802, 6.511147, -6.608413), Point(-8.50195, 4.784208, -4.394632),
+ Point(-7.697445, 5.98417, -4.444797), Point(7.923917, 4.278562, -8.69608),
+ Point(7.207772, 5.426537, -8.625711), Point(6.570754, 5.663218, -9.95051),
+ Point(0.0, -9.822524, -3.751277), Point(1.450248, -9.642715, -4.433885),
+ Point(2.586084, -9.509383, -3.396354), Point(2.321048, 9.211703, -6.247323),
+ Point(3.22564, 8.751511, -7.21285), Point(3.516393, 8.925362, -5.647262),
+ Point(-9.655634, 0.0, -5.203358), Point(-9.586456, 1.581374, -4.732493),
+ Point(-9.39311, 1.358915, -6.300104), Point(1.541872, -7.684792, -12.420403),
+ Point(0.0, -7.89342, -12.279075), Point(1.759295, -8.164936, -10.997946),
+ Point(-8.248962, 4.047046, -7.893298), Point(-9.203081, 1.071802, -7.524505),
+ Point(0.0, 9.945297, 2.089081), Point(0.0, 9.841475, 3.547035),
+ Point(-1.283953, 9.829842, 2.627298), Point(-9.339994, 3.413321, -2.110682),
+ Point(-9.793957, 1.766635, -1.956951), Point(-9.791217, 1.166692, -3.329208),
+ Point(0.0, -9.995065, 0.62828), Point(0.0, -9.945297, 2.089081),
+ Point(1.600841, -9.849996, 1.288222), Point(-5.368326, 7.209325, -8.76509),
+ Point(-5.726311, 7.310576, -7.420199), Point(-4.42404, 8.023131, -8.014298),
+ Point(-1.600844, 9.849995, 1.288256), Point(0.0, -9.991315, -0.833365),
+ Point(0.0, -9.934008, -2.293889), Point(2.019746, -9.779264, -1.070758),
+ Point(-9.991315, 0.0, -0.833365), Point(-9.995065, 0.0, 0.62828),
+ Point(-9.940412, -1.084272, -0.224247), Point(-1.278612, -9.917806, -0.095468),
+ Point(-8.362888, -5.102845, -4.011525), Point(-8.362888, -5.102845, -4.011525),
+ Point(-9.222327, -3.524133, -3.180679), Point(-4.627508, 8.248109, -6.497652),
+ Point(-5.554324, 7.923507, -5.046793), Point(-5.965308, 7.404348, -6.19378),
+ Point(-5.609637, -6.397016, -10.509075), Point(-7.207772, -5.426537, -8.625711),
+ Point(-6.116087, -6.530472, -8.932282), Point(-9.942056, 1.026171, -0.640274),
+ Point(-9.177567, 3.958378, -0.643456), Point(-9.118361, 2.76194, -6.075251),
+ Point(-9.149798, 3.28097, -4.69742), Point(2.308184, -9.729522, 0.186449),
+ Point(-0.970845, 9.754561, 3.952719), Point(0.0, -9.655634, -5.203358),
+ Point(2.172605, -9.354452, -5.576386), Point(3.943607, -9.034725, -3.359585),
+ Point(9.940412, 1.084272, -0.224247), Point(9.991315, 0.0, -0.833365),
+ Point(9.921607, 0.984645, -1.539082), Point(-6.865367, 7.045819, 3.590633),
+ Point(-5.163251, 8.491609, 2.221184), Point(-4.598726, 8.68703, 3.680882),
+ Point(-9.842685, 1.695645, 0.992653), Point(3.468073, -9.374293, -0.616795),
+ Point(3.134015, -9.481336, 1.062484), Point(-9.162213, -1.544635, -7.39404),
+ Point(-9.406783, -1.601729, -5.982272), Point(-8.886464, -3.091493, -6.774492),
+ Point(-1.437837, -9.874303, 1.312648), Point(8.049473, 1.295595, -11.580573),
+ Point(8.049473, 1.295595, -11.580573), Point(6.931408, 3.93547, -12.077691),
+ Point(5.00158, 8.154677, -5.825954), Point(3.516393, 8.925362, -5.647262),
+ Point(-9.374384, 3.461844, 0.739081), Point(-5.047792, 7.049231, -9.965568),
+ Point(-6.489734, 6.430658, -8.131417), Point(-6.73937, 5.331576, -10.228428),
+ Point(9.995065, 0.0, 0.62828), Point(0.0, 9.655634, -5.203358),
+ Point(-3.076371, 9.04, -5.937793), Point(-1.484486, 9.490365, -5.560312),
+ Point(-3.919505, 9.131221, 2.243471), Point(-3.134028, 9.481329, 1.062572),
+ Point(-7.0711, 7.002586, -1.962994), Point(-8.100178, 5.610395, -3.412089),
+ Point(8.012496, -5.660853, 3.875387), Point(6.340662, -7.317514, 5.0),
+ Point(7.317514, -6.340662, 5.0), Point(7.406436, 6.693908, 1.159845),
+ Point(4.608473, 8.8316, 1.74908), Point(6.538408, 7.56573, 0.18913),
+ Point(-6.538408, -7.56573, 0.18913), Point(-4.608473, -8.8316, 1.74908),
+ Point(-7.406436, -6.693908, 1.159845), Point(8.625067, 3.023535, -8.115776),
+ Point(4.581637, -8.878087, 0.867579), Point(3.919496, -9.131229, 2.243396),
+ Point(2.600308, -9.589656, 2.260002), Point(9.147276, 0.0, -8.081422),
+ Point(9.162213, 1.544635, -7.39404), Point(8.855374, 1.485543, -8.803527),
+ Point(8.700675, 4.638636, 3.335458), Point(8.145402, 5.234732, 5.0),
+ Point(-2.727864, 9.290251, 5.0), Point(-3.405379, 9.233434, 3.548003),
+ Point(7.622815, 6.471857, -0.176254), Point(7.035808, 7.07383, -1.353989),
+ Point(-7.622815, -6.471857, -0.176254), Point(-7.035808, -7.07383, -1.353989),
+ Point(-9.608258, -2.178892, -3.425678), Point(-9.384758, -2.603804, -4.537189),
+ Point(7.406436, 6.693908, 1.159845), Point(7.205864, 6.822549, 2.472527),
+ Point(-7.406436, -6.693908, 1.159845), Point(-8.700675, -4.638636, 3.335458),
+ Point(-7.205864, -6.822549, 2.472527), Point(9.35935, 3.0779, 3.422923),
+ Point(9.161986, 3.861816, 2.139516), Point(9.786828, 1.528015, 2.744577),
+ Point(-9.161986, -3.861816, 2.139516), Point(-9.786828, -1.528015, 2.744577),
+ Point(-9.35935, -3.0779, 3.422923), Point(9.222327, 3.524133, -3.180679),
+ Point(9.608258, 2.178892, -3.425678), Point(9.384758, 2.603804, -4.537189),
+ Point(7.389736, 4.362426, -10.2686), Point(8.886464, 3.091493, -6.774492),
+ Point(9.406783, 1.601729, -5.982272), Point(1.217776, 8.982274, -8.446485),
+ Point(-9.934008, 0.0, -2.293889), Point(-9.921607, -0.984645, -1.539082),
+ Point(9.942041, -1.026291, -0.640463), Point(-4.178561, 7.750891, -9.479097),
+ Point(-1.493749, 8.753735, -9.195832), Point(-1.787216, 8.978332, -8.048706),
+ Point(-9.822524, 0.0, -3.751277), Point(-6.402262, -5.224398, -11.263518),
+ Point(-4.400285, -7.167487, -10.819357), Point(-8.38544, 2.822228, -9.320822),
+ Point(-9.147276, 0.0, -8.081422), Point(-8.383643, 0.0, -10.902206),
+ Point(-7.89342, 0.0, -12.279075), Point(-8.41869, 1.475438, -10.382435),
+ Point(4.598732, -8.687048, 3.680679), Point(6.075063, -7.836377, 2.595997),
+ Point(-4.854151, -8.435327, -4.596731), Point(-5.840074, -8.036067, -2.293613),
+ Point(-3.603608, -9.119057, -3.92775), Point(5.840074, 8.036067, -2.293613),
+ Point(3.603608, 9.119057, -3.92775), Point(4.854151, 8.435327, -4.596731),
+ Point(-4.022239, 8.807474, 5.0), Point(-2.247424, -9.697873, -1.897723),
+ Point(-1.263633, -9.814303, -2.886306), Point(0.970846, -9.754562, 3.952713),
+ Point(3.405377, -9.233438, 3.547971), Point(2.727864, -9.290251, 5.0),
+ Point(7.436375, 5.664592, -7.102878), Point(3.822908, 1.76438, -18.14082),
+ Point(4.995853, 0.0, -17.325293), Point(5.295787, 2.384052, -16.281392),
+ Point(-4.581704, 8.878045, 0.867894), Point(-1.3838, -3.769026, -18.317155),
+ Point(-3.822908, -1.76438, -18.14082), Point(-3.63056, -3.896829, -16.927346),
+ Point(-8.807474, 4.022239, 5.0), Point(-9.322278, 3.14467, 3.581158),
+ Point(-8.01251, 5.660827, 3.875421), Point(3.63056, 3.896829, -16.927346),
+ Point(1.3838, 3.769026, -18.317155), Point(1.283954, -9.829843, 2.627281),
+ Point(8.362888, 5.102845, -4.011525), Point(7.030446, 6.98901, -2.627976),
+ Point(7.91231, 5.985194, -2.508227), Point(-9.945297, 0.0, 2.089081),
+ Point(-9.437023, 3.118292, 2.208036), Point(4.482355, 8.219773, -7.026757),
+ Point(-5.588607, -8.278922, -0.952742), Point(-3.847876, -9.228189, 0.370805),
+ Point(5.588607, 8.278922, -0.952742), Point(3.847876, 9.228189, 0.370805),
+ Point(-9.431439, 0.0, -6.647692), Point(-3.516393, -8.925362, -5.647262),
+ Point(-5.00158, -8.154677, -5.825954), Point(-3.516393, -8.925362, -5.647262),
+ Point(9.339895, -3.413559, -2.110911), Point(8.685234, -4.733525, -2.939696),
+ Point(8.42908, -5.360993, -0.9173), Point(9.655634, 0.0, -5.203358),
+ Point(9.431439, 0.0, -6.647692), Point(0.0, -4.995853, -17.325293),
+ Point(0.0, -2.779896, -19.211681), Point(-7.982385, 5.984936, 1.359515),
+ Point(3.595929, -8.207393, -8.878738), Point(3.113525, -7.939146, -10.445271),
+ Point(5.076972, -7.368264, -8.929287), Point(9.680142, 1.4641, 4.074932),
+ Point(-9.680142, -1.4641, 4.074932), Point(6.127597, 7.460083, -5.21525),
+ Point(-3.756919, 6.806889, -12.578047), Point(-2.872075, 6.73563, -13.620937),
+ Point(-1.176309, 7.86629, -12.122339), Point(7.696039, -5.987322, -4.437551),
+ Point(4.761513, -8.369654, -5.395143), Point(6.994977, -6.327682, -6.64251),
+ Point(0.0, -9.841475, 3.547035), Point(6.959392, 2.285125, -13.615442),
+ Point(5.360192, 4.642974, -14.101224), Point(6.290584, 2.421831, -14.773392),
+ Point(4.022239, -8.807474, 5.0), Point(6.001576, -7.965004, 1.469414),
+ Point(6.942243, -7.185644, 0.828909), Point(5.953604, -8.034244, 0.148704),
+ Point(0.0, 8.383643, -10.902206), Point(0.0, 9.147276, -8.081422),
+ Point(9.841475, 0.0, 3.547035), Point(-9.841475, 0.0, 3.547035),
+ Point(-3.944211, 9.034265, -3.361695), Point(-2.586191, 9.509367, -3.39621),
+ Point(-2.841712, 9.309804, -4.583546), Point(8.498357, -4.800328, -4.351912),
+ Point(9.617233, -1.510763, -4.572273), Point(9.120982, -3.714202, -3.471244),
+ Point(-6.142055, 4.485546, -12.985382), Point(-7.074396, 3.644425, -12.111331),
+ Point(-7.619494, 2.69177, -11.780947), Point(3.722538, 9.128716, 3.352165),
+ Point(2.056504, 9.60311, 3.768861), Point(2.727864, 9.290251, 5.0),
+ Point(-2.056504, -9.60311, 3.768861), Point(-3.722538, -9.128716, 3.352165),
+ Point(-2.727864, -9.290251, 5.0), Point(0.0, 4.995853, -17.325293),
+ Point(0.0, 9.431439, -6.647692), Point(8.805112, -4.653675, 1.803671),
+ Point(9.32227, -3.144705, 3.581121), Point(0.0, 9.995065, 0.62828),
+ Point(0.0, 9.991315, -0.833365), Point(1.278612, 9.917806, -0.095468),
+ Point(1.117952, 9.840235, 2.771249), Point(3.722538, 9.128716, 3.352165),
+ Point(-3.722538, -9.128716, 3.352165), Point(-1.117952, -9.840235, 2.771249),
+ Point(8.373136, 4.61219, -5.871384), Point(-6.188205, 7.719229, -2.91179),
+ Point(-9.674476, 1.513073, 4.057154), Point(9.290251, 2.727864, 5.0),
+ Point(-9.290251, -2.727864, 5.0), Point(8.09963, -5.611666, -3.408934),
+ Point(-2.108584, -9.014717, -7.560095), Point(-4.888928, -7.362989, -9.356233),
+ Point(-4.482355, -8.219773, -7.026757), Point(-7.444114, -6.094829, -5.454624),
+ Point(-6.635108, -6.779438, -6.329156), Point(6.371809, 6.680444, -7.686798),
+ Point(-8.428951, 5.361201, -0.917185), Point(-7.436375, -5.664592, -7.102878),
+ Point(2.584154, 2.082214, -18.866535), Point(-7.317514, 6.340662, 5.0),
+ Point(-6.865367, 7.045819, 3.590633), Point(7.95022, -6.065797, 0.02003),
+ Point(9.374344, -3.461958, 0.738975), Point(-8.745408, 4.764895, -1.80401),
+ Point(-5.295787, -2.384052, -16.281392), Point(-6.020811, -1.187948, -15.790962),
+ Point(-5.521678, -3.503721, -15.130764), Point(9.945297, 0.0, 2.089081),
+ Point(9.437005, -3.118358, 2.207972), Point(9.842672, -1.695729, 0.992584),
+ Point(2.468412, 8.687375, -8.58754), Point(1.483419, -3.090854, -18.787878),
+ Point(0.0, 0.0, -20.0), Point(-7.389736, -4.362426, -10.2686),
+ Point(6.635108, 6.779438, -6.329156), Point(9.674474, -1.51309, 4.057145),
+ Point(5.609637, 6.397016, -10.509075), Point(-4.714247, 3.832741, -15.885335),
+ Point(-5.521564, 4.280354, -14.309564), Point(-3.047642, 5.925253, -14.913517),
+ Point(4.400285, 7.167487, -10.819357), Point(4.888928, 7.362989, -9.356233),
+ Point(2.732696, 8.272165, -9.819093), Point(9.054826, -2.932029, -6.136229),
+ Point(7.736149, -5.670009, -5.657916), Point(7.854948, -5.072572, -7.090504),
+ Point(-1.689637, 5.529037, -16.318686), Point(-0.776441, 6.273141, -15.497721),
+ Point(8.055462, -5.814136, -2.285038), Point(6.239651, -7.680884, -2.878045),
+ Point(9.934008, 0.0, -2.293889), Point(-1.217729, -5.448249, -16.593218),
+ Point(-7.030446, -6.98901, -2.627976), Point(-9.89249, -1.290916, 1.374308),
+ Point(7.444114, 6.094829, -5.454624), Point(0.0, 9.822524, -3.751277),
+ Point(-7.923917, -4.278562, -8.69608), Point(-8.373136, -4.61219, -5.871384),
+ Point(0.0, 9.682458, 5.0), Point(-7.608541, 6.459149, -1.248201),
+ Point(-1.330507, -9.683424, -4.224239), Point(1.437837, 9.874303, 1.312648),
+ Point(-9.366259, -3.450924, -1.207167), Point(9.366259, 3.450924, -1.207167),
+ Point(-4.995853, 0.0, -17.325293), Point(6.402262, 5.224398, -11.263518),
+ Point(-2.3527, 3.654664, -18.012022), Point(3.426513, -8.631426, -7.418217),
+ Point(3.470186, -8.913809, -5.831578), Point(0.0, -9.431439, -6.647692),
+ Point(2.089361, -8.89927, -8.108651), Point(-6.033215, 7.945801, -1.36319),
+ Point(2.798414, -7.60218, -11.726165), Point(-2.584154, -2.082214, -18.866535),
+ Point(-8.145402, -5.234732, 5.0), Point(-8.807474, -4.022239, 5.0),
+ Point(-5.95371, 8.034164, 0.149047), Point(0.0, -7.321061, -13.623813),
+ Point(7.982329, -5.985019, 1.359368), Point(4.915513, -6.112277, -12.406096),
+ Point(5.631589, -4.099943, -14.34931), Point(6.328542, -4.984044, -11.850547),
+ Point(-8.515895, 4.254855, -6.123968), Point(9.682458, 0.0, 5.0),
+ Point(9.822524, 0.0, -3.751277), Point(4.339084, -5.137415, -14.802611),
+ Point(-8.805153, 4.653587, 1.803784), Point(-6.656469, 0.0, -14.925337),
+ Point(-5.886745, 0.0, -16.167403), Point(-6.20211, 1.396261, -15.438172),
+ Point(4.061444, 5.551108, -14.517558), Point(1.521275, 6.279054, -15.265543),
+ Point(2.740098, 7.153814, -12.855319), Point(-4.67775, -4.454881, -15.267311),
+ Point(-2.249598, -5.474352, -16.120891), Point(-4.061444, -5.551108, -14.517558),
+ Point(4.615667, -7.95527, -7.850933), Point(5.886745, 0.0, -16.167403),
+ Point(5.295358, -1.544342, -16.682229), Point(-8.775355, 4.55065, 3.023068),
+ Point(8.383643, 0.0, -10.902206), Point(7.89342, 0.0, -12.279075),
+ Point(8.423933, -1.472621, -10.367011), Point(-7.872452, 2.99478, -10.780684),
+ Point(8.525131, -3.881711, -7.001273), Point(-7.321061, 0.0, -13.623813),
+ Point(-6.84593, -1.18465, -14.384693), Point(-5.558497, 6.30843, -10.827155),
+ Point(-4.853359, 6.198252, -12.333139), Point(-6.942289, 7.185593, 0.829179),
+ Point(-7.91231, -5.985194, -2.508227), Point(-4.007736, 7.262206, -11.171109),
+ Point(-2.819087, 8.169654, -10.061711), Point(-2.463612, 7.895329, -11.241779),
+ Point(3.773839, 7.05385, -12.000224), Point(5.967158, 4.933297, -12.6579),
+ Point(-6.001568, 7.964978, 1.470103), Point(3.561016, 4.832041, -15.996318),
+ Point(2.249598, 5.474352, -16.120891), Point(8.744239, 4.781202, 1.647304),
+ Point(9.163301, 3.998621, 0.42399), Point(-8.744239, -4.781202, 1.647304),
+ Point(-9.163301, -3.998621, 0.42399), Point(-4.011247, 5.94065, -13.945406),
+ Point(-3.500931, -9.295193, -2.317653), Point(3.500931, 9.295193, -2.317653),
+ Point(4.050364, -9.089222, -1.980495), Point(4.919632, -8.695004, -0.881183),
+ Point(0.0, 2.779896, -19.211681), Point(6.964654, -5.940871, -8.049756),
+ Point(7.660935, -4.941604, -8.219641), Point(6.664694, -5.336561, -10.412102),
+ Point(8.861878, -1.331606, -8.875572), Point(-9.682458, 0.0, 5.0),
+ Point(0.0, -8.383643, -10.902206), Point(-9.064702, -3.145124, -5.635379),
+ Point(9.89249, 1.290916, 1.374308), Point(2.714435, 9.615937, -0.813884),
+ Point(-2.714435, -9.615937, -0.813884), Point(1.840457, -7.106791, -13.580315),
+ Point(0.0, -6.656469, -14.925337), Point(-2.102126, 1.55927, -19.302823),
+ Point(-2.779896, 0.0, -19.211681), Point(-4.355466, 1.485777, -17.756394),
+ Point(2.00948, -8.548715, -9.566916), Point(-6.340662, 7.317514, 5.0),
+ Point(-2.740098, -7.153814, -12.855319), Point(-2.780235, -6.476139, -14.188716),
+ Point(-3.929611, -6.421026, -13.16489), Point(6.340662, -7.317514, 5.0),
+ Point(-6.900495, 1.652029, -14.093115), Point(9.395315, -1.474676, -6.181709),
+ Point(-1.322352, -8.64356, -9.703661), Point(-2.934265, -7.750527, -11.192752),
+ Point(8.807474, -4.022239, 5.0), Point(-2.321048, -9.211703, -6.247323),
+ Point(0.0, -9.147276, -8.081422), Point(-6.931408, -3.93547, -12.077691),
+ Point(-7.482591, -1.191537, -13.052368), Point(-8.049473, -1.295595, -11.580573),
+ Point(-7.692051, 4.223495, -9.590506), Point(4.345972, -2.438334, -17.339787),
+ Point(2.779896, 0.0, -19.211681), Point(8.777527, 4.773804, -0.814421),
+ Point(-8.777527, -4.773804, -0.814421), Point(-4.919906, 8.694851, -0.881127),
+ Point(-3.773839, -7.05385, -12.000224), Point(0.0, -9.682458, 5.0),
+ Point(8.748052, 4.701532, -2.338539), Point(1.322352, 8.64356, -9.703661),
+ Point(9.064702, 3.145124, -5.635379), Point(9.560134, -2.415388, -3.328507),
+ Point(-8.748052, -4.701532, -2.338539), Point(2.247424, 9.697873, -1.897723),
+ Point(4.592819, 3.658445, -16.189106), Point(-8.855374, -1.485543, -8.803527),
+ Point(9.290251, -2.727864, 5.0), Point(0.0, 9.934008, -2.293889),
+ Point(1.263633, 9.814303, -2.886306), Point(-7.619912, 5.39051, -7.177561),
+ Point(1.489733, -6.13942, -15.503317), Point(-7.350205, 2.354643, -12.716941),
+ Point(6.726934, -1.473238, -14.502127), Point(6.031875, -7.946891, -1.361485),
+ Point(0.0, -5.886745, -16.167403), Point(0.0, 7.89342, -12.279075),
+ Point(4.934104, -4.012772, -15.434024), Point(3.847169, -4.433758, -16.191489),
+ Point(1.52261, -5.373954, -16.589427), Point(8.384342, -2.821454, -9.325706),
+ Point(1.330507, 9.683424, -4.224239), Point(7.906534, -3.014506, -10.65823),
+ Point(6.656469, 0.0, -14.925337), Point(0.0, 6.656469, -14.925337),
+ Point(7.608353, -6.459351, -1.248593), Point(-2.959711, 9.491615, -2.144163),
+ Point(-2.019777, 9.77926, -1.070663), Point(1.6829, 8.234859, -10.836041),
+ Point(8.807474, 4.022239, 5.0), Point(-4.022239, -8.807474, 5.0),
+ Point(4.022239, 8.807474, 5.0), Point(5.4767, -6.270199, -11.079777),
+ Point(-6.412185, -3.514513, -13.642886), Point(7.149459, -2.746444, -12.859593),
+ Point(7.317514, 6.340662, 5.0), Point(-7.317514, -6.340662, 5.0),
+ Point(-5.967158, -4.933297, -12.6579), Point(0.0, 5.886745, -16.167403),
+ Point(4.595934, -7.286997, -10.154223), Point(7.878329, -1.652154, -11.866308),
+ Point(-9.290251, 2.727864, 5.0), Point(0.0, 7.321061, -13.623813),
+ Point(7.321061, 0.0, -13.623813)
+};
+
+//FullEllipsoid
+const std::vector<Point> outer_line_y_full{
+ Point(0, 0, -20),
+ Point(1.4417481, 0, -19.791044),
+ Point(2.7798965, 0, -19.211681),
+ Point(3.9654527, 0, -18.360304),
+ Point(4.9958534, 0, -17.325293),
+ Point(5.886745, 0, -16.167404),
+ Point(6.6564689, 0, -14.925337),
+ Point(7.3210607, 0, -13.623814),
+ Point(7.8934197, 0, -12.279076),
+ Point(8.3836432, 0, -10.902205),
+ Point(8.799571, 0, -9.5010614),
+ Point(9.1472759, 0, -8.0814219),
+ Point(9.4314394, 0, -6.6476922),
+ Point(9.6556339, 0, -5.2033577),
+ Point(9.8225241, 0, -3.751277),
+ Point(9.9340086, 0, -2.2938895),
+ Point(9.9913149, 0, -0.83336478),
+ Point(9.9950647, 0, 0.62827998),
+ Point(9.9452972, 0, 2.0890813),
+ Point(9.8414755, 0, 3.5470352),
+ Point(9.6824579, 0, 5)
+};
+
+//T4_1
+const std::vector<Point> outer_line_y_1{
+ Point(0, 0, -20),
+ Point(-2.7799, 0, -19.2117),
+ Point(-4.99585, 0, -17.3253),
+ Point(-6.65647, 0, -14.9253),
+ Point(-7.89342, 0, -12.2791),
+ Point(-8.79957, 0, -9.50106),
+ Point(-9.43144, 0, -6.64769),
+ Point(-9.82252, 0, -3.75128),
+ Point(-9.99132, 0, -0.833365),
+ Point(-9.9453, 0, 2.08908),
+ Point(-9.68246, 0, 5)
+};
+
+//T4_2
+const std::vector<Point> outer_line_y_2{
+ Point(0, 0, -20),
+ Point(-1.44175, 0, -19.791),
+ Point(-2.7799, 0, -19.2117),
+ Point(-3.96545, 0, -18.3603),
+ Point(-4.99585, 0, -17.3253),
+ Point(-5.88675, 0, -16.1674),
+ Point(-6.65647, 0, -14.9253),
+ Point(-7.32106, 0, -13.6238),
+ Point(-7.89342, 0, -12.2791),
+ Point(-8.38364, 0, -10.9022),
+ Point(-8.79957, 0, -9.50106),
+ Point(-9.14728, 0, -8.08142),
+ Point(-9.43144, 0, -6.64769),
+ Point(-9.65563, 0, -5.20336),
+ Point(-9.82252, 0, -3.75128),
+ Point(-9.93401, 0, -2.29389),
+ Point(-9.99132, 0, -0.833365),
+ Point(-9.99507, 0, 0.62828),
+ Point(-9.9453, 0, 2.08908),
+ Point(-9.841148, 0, 3.54704),
+ Point(-9.68246, 0, 5)
+};
+
+//T4_3
+const std::vector<Point> outer_line_y_3{
+ Point(0, 0, -20),
+ Point(-0.968432, 0, -19.906),
+ Point(-1.90301, 0, -19.6345),
+ Point(-2.7799, 0, -19.2117),
+ Point(-3.58799, 0, -18.6683),
+ Point(-4.32555, 0, -18.0322),
+ Point(-4.99585, 0, -17.3253),
+ Point(-5.60412, 0, -16.5643),
+ Point(-6.15593, 0, -15.7613),
+ Point(-6.65647, 0, -14.9253),
+ Point(-7.11035, 0, -14.0631),
+ Point(-7.52154, 0, -13.1797),
+ Point(-7.89342, 0, -12.2791),
+ Point(-8.22884, 0, -11.3642),
+ Point(-8.53021, 0, -10.4375),
+ Point(-8.79957, 0, -9.50106),
+ Point(-9.03864, 0, -8.5564),
+ Point(-9.24886, 0, -7.60489),
+ Point(-9.43144, 0, -6.64769),
+ Point(-9.58738, 0, -5.68579),
+ Point(-9.71752, 0, -4.72006),
+ Point(-9.82252, 0, -3.75128),
+ Point(-9.90291, 0, -2.78014),
+ Point(-9.95909, 0, -1.8073),
+ Point(-9.99132, 0, -0.833365),
+ Point(-9.99975, 0, 0.141062),
+ Point(-9.98444, 0, 1.1154),
+ Point(-9.9453, 0, 2.08908),
+ Point(-9.88215, 0, 3.06149),
+ Point(-9.79468, 0, 4.03202),
+ Point(-9.68246, 0, 5)
+};
+
+//T4_4
+const std::vector<Point> outer_line_y_4{
+ Point(0, 0, -20),
+ Point(-2.17971, 0, -19.5191),
+ Point(-4.05008, 0, -18.2863),
+ Point(-5.55941, 0, -16.6244),
+ Point(-6.7652, 0, -14.7285),
+ Point(-7.72648, 0, -12.6967),
+ Point(-8.48597, 0, -10.5808),
+ Point(-9.07283, 0, -8.41042),
+ Point(-9.50667, 0, -6.20427),
+ Point(-9.80049, 0, -3.97508),
+ Point(-9.96241, 0, -1.73241),
+ Point(-9.99667, 0, 0.515851),
+ Point(-9.90415, 0, 2.76246),
+ Point(-9.68246, 0, 5)
+};
+
+//T4_5
+const std::vector<Point> outer_line_y_5{
+ Point(0, 0, -20),
+ Point(-2.30478, 0, -19.4615),
+ Point(-4.25345, 0, -18.1006),
+ Point(-5.80168, 0, -16.2899),
+ Point(6.44792, 0, -15.2872),
+ Point(-7.02131, 0, -14.241),
+ Point(-7.97884, 0, -12.0562),
+ Point(-8.7204, 0, -9.78868),
+ Point(-9.38922, 0, -6.88258),
+ Point(-9.74097, 0, -4.52258),
+ Point(-9.94232, 0, -2.14497),
+ Point(-9.99928, 0, 0.240498),
+ Point(-9.91349, 0, 2.6251),
+ Point(-9.68246, 0, 5)
+};
+
+// FullEllipsoid
+const std::vector<Point> mid_line_y_full{
+ Point(0, 0, -18.5),
+ Point(1.4633734, 0, -18.223959),
+ Point(2.7654362, 0, -17.494207),
+ Point(3.8637159, 0, -16.47974),
+ Point(4.7813883, 0, -15.297925),
+ Point(5.5508218, 0, -14.01402),
+ Point(6.1982775, 0, -12.664215),
+ Point(6.7430954, 0, -11.269672),
+ Point(7.1992874, 0, -9.8435955),
+ Point(7.5770636, 0, -8.3947229),
+ Point(7.8839121, 0, -6.9291573),
+ Point(8.1253262, 0, -5.4513831),
+ Point(8.3052893, 0, -3.9648604),
+ Point(8.4265852, 0, -2.4723966),
+ Point(8.4910088, 0, -0.9763937),
+ Point(8.4994879, 0, 0.52097338),
+ Point(8.4521418, 0, 2.0176148),
+ Point(8.3482924, 0, 3.5113959),
+ Point(8.1864214, 0, 5)
+};
+
+//T4_2
+const std::vector<Point> mid_line_y_2{
+ Point(0, 0, -18.5),
+ Point(-1.46337, 0, -18.224),
+ Point(-2.76544, 0, -17.4942),
+ Point(-3.86372, 0, -16.4797),
+ Point(-4.78139, 0, -15.2979),
+ Point(-5.55082, 0, -14.014),
+ Point(-6.19828, 0, -12.6642),
+ Point(-6.7431, 0, -11.2697),
+ Point(-7.19929, 0, -9.8436),
+ Point(-7.57706, 0, -8.39472),
+ Point(-7.88391, 0, -6.92916),
+ Point(-8.12533, 0, -5.45138),
+ Point(-8.30529, 0, -3.96486),
+ Point(-8.42658, 0, -2.4724),
+ Point(-8.49101, 0, -0.976394),
+ Point(-8.49949, 0, 0.520973),
+ Point(-8.45214, 0, 2.01762),
+ Point(-8.34829, 0, 3.5114),
+ Point(-8.18642, 0, 5)
+};
+
+//T4_4
+const std::vector<Point> mid_line_y_4{
+ Point(0, 0, -18.5),
+ Point(-2.1398, 0, -17.9046),
+ Point(-3.86372, 0, -16.4797),
+ Point(-5.18282, 0, -14.666),
+ Point(-6.19828, 0, -12.6642),
+ Point(-6.98156, 0, -10.56),
+ Point(-7.57706, 0, -8.39472),
+ Point(-8.01252, 0, -6.19156),
+ Point(-8.30529, 0, -3.96486),
+ Point(-8.46583, 0, -1.7247),
+ Point(-8.49949, 0, 0.520973),
+ Point(-8.40735, 0, 2.765),
+ Point(-8.18642, 0, 5)
+};
+
+//FullEllipsoid
+const std::vector<Point> inner_line_y_full{
+ Point(0, 0, -17),
+ Point(1.3980633, 0, -16.65749),
+ Point(2.5715427, 0, -15.811318),
+ Point(3.5111878, 0, -14.706712),
+ Point(4.2705364, 0, -13.469821),
+ Point(4.8926024, 0, -12.158011),
+ Point(5.405972, 0, -10.799763),
+ Point(5.8295898, 0, -9.4107952),
+ Point(6.1763358, 0, -8.0006142),
+ Point(6.4551773, 0, -6.5754151),
+ Point(6.6724353, 0, -5.1395173),
+ Point(6.8325491, 0, -3.6961164),
+ Point(6.9385433, 0, -2.2477272),
+ Point(6.9923134, 0, -0.79645586),
+ Point(6.9947896, 0, 0.65581053),
+ Point(6.9460135, 0, 2.1072586),
+ Point(6.8451443, 0, 3.5560143),
+ Point(6.6903844, 0, 5)
+};
+
+//T4_1
+const std::vector<Point> inner_line_y_1{
+ Point(0, 0, -17.0),
+ Point(-2.45361, 0, -15.9215),
+ Point(-4.11491, 0, -13.7526),
+ Point(-5.24561, 0, -11.2565),
+ Point(-6.03107, 0, -8.62953),
+ Point(-6.55904, 0, -5.93835),
+ Point(-6.87378, 0, -3.21376),
+ Point(-6.99728, 0, -0.473752),
+ Point(-6.9374, 0, 2.2684),
+ Point(-6.69039, 0, 5)
+};
+
+//T4_2
+const std::vector<Point> inner_line_y_2{
+ Point(0, 0, -17.0),
+ Point(-1.39806, 0, -16.6575),
+ Point(-2.57154, 0, -15.8113),
+ Point(-3.51119, 0, -14.7067),
+ Point(-4.27054, 0, -13.4698),
+ Point(-4.8926, 0, -12.158),
+ Point(-5.40597, 0, -10.7998),
+ Point(-5.82959, 0, -9.4108),
+ Point(-6.17634, 0, -8.00061),
+ Point(-6.45518, 0, -6.57541),
+ Point(-6.67244, 0, -5.13952),
+ Point(-6.83255, 0, -3.69612),
+ Point(-6.93854, 0, -2.24773),
+ Point(-6.99231, 0, -0.796456),
+ Point(-6.99479, 0, 0.655811),
+ Point(-6.94601, 0, 2.10726),
+ Point(-6.84514, 0, 3.55601),
+ Point(-6.69039, 0, 5)
+};
+
+//T4_3
+const std::vector<Point> inner_line_y_3{
+ Point(0, 0, -17.0),
+ Point(-0.969411, 0, -16.8362),
+ Point(-1.85011, 0, -16.3955),
+ Point(-2.61328, 0, -15.7709),
+ Point(-3.26852, 0, -15.033),
+ Point(-3.83378, 0, -14.2237),
+ Point(-4.3249, 0, -13.3671),
+ Point(-4.7541, 0, -12.4779),
+ Point(-5.13058, 0, -11.565),
+ Point(-5.46125, 0, -10.6345),
+ Point(-5.75136, 0, -9.69053),
+ Point(-6.005, 0, -8.73612),
+ Point(-6.22532, 0, -7.77346),
+ Point(-6.41484, 0, -6.80426),
+ Point(-6.57552, 0, -5.82986),
+ Point(-6.70892, 0, -4.85135),
+ Point(-6.81624, 0, -3.86963),
+ Point(-6.89843, 0, -2.88549),
+ Point(-6.95616, 0, -1.89961),
+ Point(-6.98991, 0, -0.912618),
+ Point(-6.99993, 0, 0.074901),
+ Point(-6.98632, 0, 1.06238),
+ Point(-6.94896, 0, 2.04924),
+ Point(-6.88755, 0, 3.0349),
+ Point(-6.8016, 0, 4.01872),
+ Point(-6.69039, 0, 5)
+};
+
+//T4_4
+const std::vector<Point> inner_line_y_4{
+ Point(0, 0, -17.0),
+ Point(-0.74007, 0, -16.9047),
+ Point(-1.43739, 0, -16.6377),
+ Point(-2.06986, 0, -16.2398),
+ Point(-2.63463, 0, -15.7499),
+ Point(-3.13765, 0, -15.1966),
+ Point(-3.5869, 0, -14.5985),
+ Point(-3.9898, 0, -13.9683),
+ Point(-4.35257, 0, -13.3141),
+ Point(-4.68022, 0, -12.6415),
+ Point(-4.97682, 0, -11.9547),
+ Point(-5.24561, 0, -11.2565),
+ Point(-5.48923, 0, -10.5492),
+ Point(-5.70984, 0, -9.83429),
+ Point(-5.90921, 0, -9.11319),
+ Point(-6.08882, 0, -8.38692),
+ Point(-6.24989, 0, -7.6563),
+ Point(-6.39344, 0, -6.92204),
+ Point(-6.52032, 0, -6.18472),
+ Point(-6.63125, 0, -5.44482),
+ Point(-6.72683, 0, -4.70279),
+ Point(-6.80753, 0, -3.95899),
+ Point(-6.87378, 0, -3.21376),
+ Point(-6.92588, 0, -2.46741),
+ Point(-6.8016, 0, 4.01872),
+ Point(-6.96407, 0, -1.72022),
+ Point(-6.98854, 0, -0.972457),
+ Point(-6.99939, 0, -0.224369),
+ Point(-6.99668, 0, 0.523793),
+ Point(-6.98038, 0, 1.27178),
+ Point(-6.95044, 0, 2.01935),
+ Point(-6.90671, 0, 2.76624),
+ Point(-6.8488, 0, 3.51217),
+ Point(-6.77699, 0, 4.25687),
+ Point(-6.69039, 0, 5)
+};
+
+//T4_5
+const std::vector<Point> inner_line_y_5{
+ Point(0, 0, -17.0),
+ Point(-0.583868, 0, -16.9408),
+ Point(-1.14574, 0, -16.7707),
+ Point(-1.67117, 0, -16.5084),
+ Point(-2.1546, 0, -16.1747),
+ Point(-2.59643, 0, -15.7873),
+ Point(-2.99977, 0, -15.3599),
+ Point(-3.36851, 0, -14.9022),
+ Point(-3.70645, 0, -14.4213),
+ Point(-4.01698, 0, -13.9223),
+ Point(-4.303, 0, -13.4088),
+ Point(-4.56696, 0, -12.8835),
+ Point(-4.81093, 0, -12.3487),
+ Point(-5.03666, 0, -11.806),
+ Point(-5.24561, 0, -11.2565),
+ Point(-5.43903, 0, -10.7014),
+ Point(-5.61799, 0, -10.1415),
+ Point(-5.7834, 0, -9.57741),
+ Point(-5.93605, 0, -9.00973),
+ Point(-6.07662, 0, -8.43895),
+ Point(-6.2057, 0, -7.86545),
+ Point(-6.3238, 0, -7.28959),
+ Point(-6.43136, 0, -6.71167),
+ Point(-6.52877, 0, -6.13195),
+ Point(-6.61636, 0, -5.55066),
+ Point(-6.69442, 0, -4.96802),
+ Point(-6.76321, 0, -4.38421),
+ Point(-6.82294, 0, -3.79941),
+ Point(-6.87378, 0, -3.21376),
+ Point(-6.91589, 0, -2.62743),
+ Point(-6.94939, 0, -2.04053),
+ Point(-6.97438, 0, -1.45322),
+ Point(-6.99092, 0, -0.865601),
+ Point(-6.99906, 0, -0.277809),
+ Point(-6.99884, 0, 0.310039),
+ Point(-6.99023, 0, 0.897825),
+ Point(-6.97323, 0, 1.48543),
+ Point(-6.94778, 0, 2.07272),
+ Point(-6.9138, 0, 2.65959),
+ Point(-6.87122, 0, 3.24589),
+ Point(-6.81989, 0, 3.8315),
+ Point(-6.75968, 0, 4.41625),
+ Point(-6.69039, 0, 5)
+};
\ No newline at end of file
diff --git a/src/elasticity/problem/EllipsoidPoints.hpp b/src/elasticity/problem/EllipsoidPoints.hpp
index 72463020dc379960e434cc556798df42403be05a..6bbd222f7115b0541995e012a01f4a3e5402d769 100644
--- a/src/elasticity/problem/EllipsoidPoints.hpp
+++ b/src/elasticity/problem/EllipsoidPoints.hpp
@@ -3,7 +3,7 @@
#include "Point.hpp"
-std::vector<Point> ellipsoidEndocard = {
+static std::vector<Point> ellipsoidEndocard = {
Point(3.764178, 5.744528, 3.286529), Point(4.171387, 5.230753, 5.0),
Point(5.066103, 4.599616, 3.584207), Point(-4.171387, -5.230753, 5.0),
Point(-5.066103, -4.599616, 3.584207), Point(-3.764178, -5.744528, 3.286529),
@@ -235,7 +235,7 @@ std::vector<Point> ellipsoidEndocard = {
Point(5.199357, -1.400357, -10.862445), Point(-6.522643, 1.488751, 5.0),
Point(0.0, 4.892602, -12.158012), Point(4.892602, 0.0, -12.158012)
};
-std::vector<Point> ellipsoidEpicard = {
+static const std::vector<Point> ellipsoidEpicard{
Point(4.928, 8.473734, 3.955074), Point(6.340662, 7.317514, 5.0),
Point(6.988606, 6.879922, 3.912063), Point(-6.340662, -7.317514, 5.0),
Point(-6.988606, -6.879922, 3.912063), Point(-4.928, -8.473734, 3.955074),
@@ -470,7 +470,7 @@ std::vector<Point> ellipsoidEpicard = {
};
//FullEllipsoid
-const std::vector<Point> outer_line_y_full{
+static const std::vector<Point> outer_line_y_full{
Point(0, 0, -20),
Point(1.4417481, 0, -19.791044),
Point(2.7798965, 0, -19.211681),
@@ -496,7 +496,7 @@ const std::vector<Point> outer_line_y_full{
//T4_1
-const std::vector<Point> outer_line_y_1{
+static const std::vector<Point> outer_line_y_1{
Point(0, 0, -20),
Point(-2.7799, 0, -19.2117),
Point(-4.99585, 0, -17.3253),
@@ -511,7 +511,7 @@ const std::vector<Point> outer_line_y_1{
};
//T4_2
-const std::vector<Point> outer_line_y_2{
+static const std::vector<Point> outer_line_y_2{
Point(0, 0, -20),
Point(-1.44175, 0, -19.791),
Point(-2.7799, 0, -19.2117),
@@ -536,7 +536,7 @@ const std::vector<Point> outer_line_y_2{
};
//T4_3
-const std::vector<Point> outer_line_y_3{
+static const std::vector<Point> outer_line_y_3{
Point(0, 0, -20),
Point(-0.968432, 0, -19.906),
Point(-1.90301, 0, -19.6345),
@@ -572,7 +572,7 @@ const std::vector<Point> outer_line_y_3{
//T4_4
-const std::vector<Point> outer_line_y_4{
+static const std::vector<Point> outer_line_y_4{
Point(0, 0, -20),
Point(-2.17971, 0, -19.5191),
Point(-4.05008, 0, -18.2863),
@@ -590,7 +590,7 @@ const std::vector<Point> outer_line_y_4{
};
//T4_5
-const std::vector<Point> outer_line_y_5{
+static const std::vector<Point> outer_line_y_5{
Point(0, 0, -20),
Point(-2.30478, 0, -19.4615),
Point(-4.25345, 0, -18.1006),
@@ -608,7 +608,7 @@ const std::vector<Point> outer_line_y_5{
};
// FullEllipsoid
-const std::vector<Point> mid_line_y_full{
+static const std::vector<Point> mid_line_y_full{
Point(0, 0, -18.5),
Point(1.4633734, 0, -18.223959),
Point(2.7654362, 0, -17.494207),
@@ -632,7 +632,7 @@ const std::vector<Point> mid_line_y_full{
//T4_2
-const std::vector<Point> mid_line_y_2{
+static const std::vector<Point> mid_line_y_2{
Point(0, 0, -18.5),
Point(-1.46337, 0, -18.224),
Point(-2.76544, 0, -17.4942),
@@ -655,7 +655,7 @@ const std::vector<Point> mid_line_y_2{
};
//T4_4
-const std::vector<Point> mid_line_y_4{
+static const std::vector<Point> mid_line_y_4{
Point(0, 0, -18.5),
Point(-2.1398, 0, -17.9046),
Point(-3.86372, 0, -16.4797),
@@ -672,7 +672,7 @@ const std::vector<Point> mid_line_y_4{
};
//FullEllipsoid
-const std::vector<Point> inner_line_y_full{
+static const std::vector<Point> inner_line_y_full{
Point(0, 0, -17),
Point(1.3980633, 0, -16.65749),
Point(2.5715427, 0, -15.811318),
@@ -694,7 +694,7 @@ const std::vector<Point> inner_line_y_full{
};
//T4_1
-const std::vector<Point> inner_line_y_1{
+static const std::vector<Point> inner_line_y_1{
Point(0, 0, -17.0),
Point(-2.45361, 0, -15.9215),
Point(-4.11491, 0, -13.7526),
@@ -708,7 +708,7 @@ const std::vector<Point> inner_line_y_1{
};
//T4_2
-const std::vector<Point> inner_line_y_2{
+static const std::vector<Point> inner_line_y_2{
Point(0, 0, -17.0),
Point(-1.39806, 0, -16.6575),
Point(-2.57154, 0, -15.8113),
@@ -730,7 +730,7 @@ const std::vector<Point> inner_line_y_2{
};
//T4_3
-const std::vector<Point> inner_line_y_3{
+static const std::vector<Point> inner_line_y_3{
Point(0, 0, -17.0),
Point(-0.969411, 0, -16.8362),
Point(-1.85011, 0, -16.3955),
@@ -760,7 +760,7 @@ const std::vector<Point> inner_line_y_3{
};
//T4_4
-const std::vector<Point> inner_line_y_4{
+static const std::vector<Point> inner_line_y_4{
Point(0, 0, -17.0),
Point(-0.74007, 0, -16.9047),
Point(-1.43739, 0, -16.6377),
@@ -799,7 +799,7 @@ const std::vector<Point> inner_line_y_4{
};
//T4_5
-const std::vector<Point> inner_line_y_5{
+static const std::vector<Point> inner_line_y_5{
Point(0, 0, -17.0),
Point(-0.583868, 0, -16.9408),
Point(-1.14574, 0, -16.7707),
diff --git a/src/elasticity/solvers/IterativePressureSolver.cpp b/src/elasticity/solvers/IterativePressureSolver.cpp
index 0977fb7a7d424e4aeb876ff96843989bdd8c8705..01b48e186f6f0307c6d9850309bb18980c3f8b23 100644
--- a/src/elasticity/solvers/IterativePressureSolver.cpp
+++ b/src/elasticity/solvers/IterativePressureSolver.cpp
@@ -10,7 +10,6 @@ void IterativePressureSolver::Initialize(IElasticity &assemble, Vector &u) {
solver.Initialize(assemble, u);
if(vtk>0){
- //Plotting::Instance().Clear();
assemble.PlotPressureIteration(u, 0);
}
assemble.PrintPressureIteration(u);
diff --git a/src/electrophysiology/MainMonodomain.cpp b/src/electrophysiology/MainMonodomain.cpp
index c77729dbb0b616aee1d5ced03f14eddac8f0a95b..d87affab74f7ef09bbd16fa3b70d728bc4ee8877 100644
--- a/src/electrophysiology/MainMonodomain.cpp
+++ b/src/electrophysiology/MainMonodomain.cpp
@@ -52,6 +52,9 @@ void MainMonodomain::Initialize() {
}
Vector &MainMonodomain::Run() {
+ const Meshes &mesh = (*potential).GetMeshes();
+ LagrangeDiscretization disc0(mesh, 0, 1);
+ const Vector sizeCells(disc0);
auto elphySolver = GetElphySolver(modelName, *cellModel, elphyA->ExcitationData());
elphySolver->Method(*elphyA, *potential);
diff --git a/src/electrophysiology/assemble/IElphyAssemble.cpp b/src/electrophysiology/assemble/IElphyAssemble.cpp
index 476452f02d49286572f5398d6c49499b9913a384..ee300fb2870054e57c04d2e1f71223d42a6ac1b5 100644
--- a/src/electrophysiology/assemble/IElphyAssemble.cpp
+++ b/src/electrophysiology/assemble/IElphyAssemble.cpp
@@ -8,6 +8,9 @@
#include "CardiacInterpolation.hpp"
#include "MCellModel.hpp"
+#include <fstream>
+#include <iostream>
+
IElphyAssemble::IElphyAssemble(ElphyProblem &elphyProblem, int degree, int exactUpTo)
: INonLinearTimeAssemble(), ICardiacAssemble("ElphyVTK"), problem(elphyProblem), degree(degree),
disc(elphyProblem.GetMeshes(), degree, exactUpTo < 0 ? 2 * degree : exactUpTo, DomainPart, 1),
@@ -18,7 +21,6 @@ IElphyAssemble::IElphyAssemble(ElphyProblem &elphyProblem, int degree, int exact
Config::Get("IsExternalCurrentSmooth", iextContinuous);
Config::Get("ThresholdVForActivation", Vthresh);
Config::Get("PrintVSteps", printVSteps);
- Config::Get("PlotVTK", vtksteps);
Config::Get("StartTime", startTime);
}
@@ -101,7 +103,7 @@ void IElphyAssemble::PrintV(const Vector &V) const {
}
cnt = PPM->Sum(cnt);
if (cnt) {
- vout(0) << " v" << P << " = " << PPM->Sum(v) / cnt << endl;
+ vout(0) << " v" <<i<< P << " = " << PPM->Sum(v) / cnt << endl;
}
}
ExchangeBuffer B;
@@ -160,12 +162,45 @@ void IElphyAssemble::PlotExcitation() const {
void IElphyAssemble::PlotActivation() const {
if (problem.UpdateActivationTime == false) return;
-
+ writeActivationTimeFile();
auto &plot = mpp::plot("ActivationTime");
plot.AddData("Activation", *activationTime);
plot.PlotFile("ActivationTime");
}
+void IElphyAssemble::writeActivationTimeFile() const {
+ string lLevel;
+ Config::Get("ElphyLevel", lLevel);
+ string tLevel="0";
+ Config::Get("ElphyTimeLevel", tLevel);
+ string filename("data/AT_l"+lLevel+"j"+ tLevel+".txt");
+ std::fstream file_out;
+ ExchangeBuffer B;
+ for (row r= (*activationTime).rows(); r != (*activationTime).rows_end(); ++r) {
+ if((*activationTime).find_procset(r()) == (*activationTime).procsets_end() || (*activationTime).find_procset(r()).master() == PPM->Proc()) {
+ //if((*activationTime).find_procset(r()).master()==PPM->Proc()) {
+ B.Send(0) << r() << (*activationTime)(r(), 0);
+ }
+ }
+ B.Communicate();
+
+ if (PPM->Master()) {
+ file_out.open(filename, std::ios_base::out);
+ if (!file_out.is_open()) {
+ std::cout << "failed to open " << filename << '\n';
+ } else {
+ for (int i = 0; i < PPM->Size(); i++) {
+ while(B.Receive(i).size()<B.Receive(i).Size()) {
+ Point p;
+ double actiTime;
+ B.Receive(i) >> p >> actiTime;
+ file_out << p << " " << actiTime << endl;
+ }
+ }
+ }
+ file_out.close();
+ }
+}
void IElphyAssemble::PrintActivationTime() const {
auto &AT = *activationTime;
if (problem.UpdateActivationTime) {
@@ -353,13 +388,17 @@ CreateElphyAssemble(ElphyProblem &problem, MCellModel &cellModel, const string &
}*/else if (modelName == "ImplictEuler") {
return std::make_unique<MonodomainFullNewton>(problem, cellModel, discDegree,
hexaQuadExactupTo);
- } else if (modelName == "SemiImplicit") {
+ } else if (modelName == "SemiImplicit" ) {
return std::make_unique<MonodomainSemiImplicit>(problem, cellModel, discDegree,
hexaQuadExactupTo);
} else if (modelName == "SemiImplicitNodes") {
return std::make_unique<MonodomainSemiImplicitNodes>(problem, cellModel, discDegree,
hexaQuadExactupTo);
- } else {
+ }else if (modelName == "SemiImplicitOnCells") {
+ return std::make_unique<MonodomainSemiImplicitCells>(problem, cellModel, discDegree,
+ hexaQuadExactupTo);
+ }
+ else {
return std::make_unique<MonodomainSplitting>(problem, discDegree, hexaQuadExactupTo);
}
}
diff --git a/src/electrophysiology/assemble/IElphyAssemble.hpp b/src/electrophysiology/assemble/IElphyAssemble.hpp
index 4ac02495a1b575b3a3202cb77b6cbdc4ee7162b5..ac2aff06cdd96b17f0676ca75115f0faefe1629a 100644
--- a/src/electrophysiology/assemble/IElphyAssemble.hpp
+++ b/src/electrophysiology/assemble/IElphyAssemble.hpp
@@ -50,12 +50,15 @@ public:
void FinishTimeStep(const Vector &u) override;
- virtual void RHS(const Vector &u, Vector &A) const = 0;
+ virtual void RHS(const Vector &u, Vector &r) const {
+ THROW("Not implemented in the basis class")
+ };
virtual void RHS(const Vectors &u, Vector &r) const {
THROW("Not implemented in the basis class")
};
+
virtual void RHSMAT(const Vector &u, Matrix &A) const {
THROW("Not implemented in the basis class")
};
@@ -116,6 +119,8 @@ public:
void PlotActivation() const;
+ void writeActivationTimeFile() const;
+
void PlotPotential(Vector &V, const std::string &varname, int step) const;
Vector &ActivationTime() const;
diff --git a/src/electrophysiology/assemble/Monodomain.cpp b/src/electrophysiology/assemble/Monodomain.cpp
index 7e6184b3d3670e29876cc424644c31efdc10cf16..c4328b94f3c9fe80b49020704d57d20e4381b447 100644
--- a/src/electrophysiology/assemble/Monodomain.cpp
+++ b/src/electrophysiology/assemble/Monodomain.cpp
@@ -496,6 +496,7 @@ void MonodomainFullNewton::updateVCW(const Vectors &VCW) {
void MonodomainSemiImplicit::RHS(const cell &c, const Vectors &VCW, Vector &r) const {
+ Point P(2.125, 2.75, 2.875);
double J = problem.Jacobian(*c);
double J_old = problem.Jacobian_old(*c);
MultiPartScalarElementWG E(VCW[0], *c, DomainPart(*c));
@@ -506,10 +507,15 @@ void MonodomainSemiImplicit::RHS(const cell &c, const Vectors &VCW, Vector &r) c
for (int j = 0; j < E.size(); ++j) {
E.Values(q, VCW, vcw, E.IndexPart(j));
double Iion = elphyModel.GetIion(vcw);
+ /*if (c()==P ){
+ mout<<"iion "<<Iion<<endl;
+ }*/
+
double iext = 0.0;
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
+
Scalar V_j = E.Value(q, j);
r_c(j) += w * (J_old * vcw[0] -
J * elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
@@ -521,18 +527,18 @@ void MonodomainSemiImplicit::RHS(const cell &c, const Vectors &VCW, Vector &r) c
void MonodomainSemiImplicitNodes::RHS(const cell &c, const Vectors &VCW, Vector &r) const {
MultiPartScalarElementWG E(VCW[0], *c, DomainPart(*c));
PartRowBndValues r_c(r, E, *c, DomainPart(*c));
- std::vector<double> vcw(elphyModel.ElphySize());
+
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- E.Values(q, VCW, vcw, E.IndexPart(j));
+ Scalar v = E.Value(q, VCW[0],E.IndexPart(j));
double Iion = E.Value(q, *IionVEC, E.IndexPart(j));
double iext = 0.0;
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
Scalar V_j = E.Value(q, j);
- r_c(j) += w * (vcw[0] -
+ r_c(j) += w * (v -
elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
(Iion - iext)) * V_j;
}
@@ -551,4 +557,42 @@ void MonodomainSemiImplicitNodes::updateIionVecs(const Vectors &VCW) {
}
}
vout(10).EndBlock();
+}
+void MonodomainSemiImplicitCells::RHS(const cell &c,const Vector &V, Vector &r) const {
+ Point P(2.125, 2.75, 2.875);
+ double J = problem.Jacobian(*c);
+ double J_old = problem.Jacobian_old(*c);
+ MultiPartScalarElementWG E(V, *c, DomainPart(*c));
+ PartRowBndValues r_c(r, E, *c, DomainPart(*c));
+ for (int q = 0; q < E.nQ(); ++q) {
+ double w = E.QWeight(q);
+ for (int j = 0; j < E.size(); ++j) {
+ Scalar v = E.Value(q, V,E.IndexPart(j));
+ double Iion = (*IionCell)(c(), E.IndexPart(j));
+ /*if(c()==P)
+ mout<<"iion "<<Iion<<endl;*/
+ double iext = 0.0;
+ if (!isNormIextZero) {
+ iext = E.Value(q, *Iext, E.IndexPart(j));
+ }
+ Scalar V_j = E.Value(q, j);
+ r_c(j) += w * (J_old * v -
+ J * elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
+ (Iion - iext)) * V_j;
+ }
+ }
+}
+void MonodomainSemiImplicitCells::updateIionVecs(const Vectors &VCW) {
+ vout(10).StartBlock("Update Iion and DVIion");
+ std::vector<double> vcw(elphyModel.ElphySize());
+ for (row r = VCW[0].rows(); r != VCW[0].rows_end(); ++r) {
+ for (int j = 0; j < r.n(); ++j) {
+ for (int i = 0; i < elphyModel.ElphySize(); i++) {
+ vcw[i] = VCW[i](r, j);
+ }
+ (*IionCell)(r, j) = elphyModel.GetIion(vcw);
+ }
+ }
+ vout(10).EndBlock();
+ //PrintVariable(*IionCell,"ionCell nach update");
}
\ No newline at end of file
diff --git a/src/electrophysiology/assemble/Monodomain.hpp b/src/electrophysiology/assemble/Monodomain.hpp
index f9f7892ea8c2f2457fceda4340ad8ad3b11711ab..995679947dcef35d55f4a6959853e478a5424e3f 100644
--- a/src/electrophysiology/assemble/Monodomain.hpp
+++ b/src/electrophysiology/assemble/Monodomain.hpp
@@ -35,6 +35,11 @@ public:
virtual void RHS(const cell &c, const Vectors &u, Vector &r) const;
+ virtual void RHS(const cell &c, const Vector &u, Vector &r) const{
+ THROW("Not implemented in the basis Monodomain")
+ };
+
+
void RHS(const Vectors &u, Vector &r) const override {
r = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c)
@@ -44,7 +49,11 @@ public:
}
void RHS(const Vector &u, Vector &r) const override {
- THROW("Not implemented for this assemble class!")
+ r = 0;
+ for (cell c = u.cells(); c != u.cells_end(); ++c)
+ RHS(c, u, r);
+ //r.ClearDirichletValues();
+ r.Collect();
};
virtual void addReactionSAndRHS(const cell &c, const Vectors &u, Matrix &S, Vector &r) const;
@@ -188,6 +197,25 @@ public:
void updateIionVecs(const Vectors &VCW) override;
+};
+class MonodomainSemiImplicitCells : public Monodomain {
+
+ std::unique_ptr<LagrangeDiscretization> disc0;
+ std::unique_ptr<Vector> IionCell;
+public:
+ MonodomainSemiImplicitCells(ElphyProblem &elphyProblem, MCellModel &cellModel, int degree,
+ int exactUpTo = -1)
+ : Monodomain(elphyProblem, cellModel, degree, exactUpTo) {
+ const Meshes &mesh = GetDisc().GetMeshes();
+ disc0=std::make_unique<LagrangeDiscretization >(mesh, 0, 1);
+ IionCell = std::make_unique<Vector>(0.0,*disc0);;
+
+ };
+
+ void RHS(const cell &c, const Vector &u, Vector &r) const override;
+
+ void updateIionVecs(const Vectors &VCW) override;
+
};
#endif
\ No newline at end of file
diff --git a/src/electrophysiology/problem/CMakeLists.txt b/src/electrophysiology/problem/CMakeLists.txt
index 89819b16feab9cbb308c7ca7489af6df35bb5c04..ad065edc06f4534f93b8986d0cdcab8965fad4bf 100644
--- a/src/electrophysiology/problem/CMakeLists.txt
+++ b/src/electrophysiology/problem/CMakeLists.txt
@@ -1,4 +1,4 @@
-add_library(ELPHY_PROBLEMS STATIC
+add_library(ELPHY_PROBLEMS SHARED
ElphyProblem.cpp
VariableMeshProblem.cpp
ElphyNoDiffusionProblem.cpp
diff --git a/src/electrophysiology/problem/ElphyBenchmarkProblem.cpp b/src/electrophysiology/problem/ElphyBenchmarkProblem.cpp
index f4b2afbb5d07aec26896a4b8a057f176f986bcb2..9328b58f0c31fa6d3055c917949c60f4244da4b8 100644
--- a/src/electrophysiology/problem/ElphyBenchmarkProblem.cpp
+++ b/src/electrophysiology/problem/ElphyBenchmarkProblem.cpp
@@ -14,4 +14,15 @@ void ElphyBenchmarkProblem::InitializeEvaluationPoints() {
evaluationPoints.emplace_back(Point(20.0, 0.0, 3.0));
evaluationPoints.emplace_back(Point(20.0, 7.0, 3.0));
evaluationPoints.emplace_back(Point(10.0, 3.5, 1.5));
+ //Eckpunkte der Zelle P(2.125, 2.75, 2.875)
+ /*evaluationPoints.emplace_back(Point(2, 2.5, 2.5));
+ evaluationPoints.emplace_back(Point(2, 3, 3));
+ evaluationPoints.emplace_back(Point(2, 2.5, 3));
+ evaluationPoints.emplace_back(Point(2.5, 3, 3));*/
+
+ evaluationPoints.emplace_back(Point(0.125, 0.75, 0.875));
+ evaluationPoints.emplace_back(Point(2.125, 2.75, 2.875));
+ evaluationPoints.emplace_back(Point(10.125, 2.75, 1.875));
+ evaluationPoints.emplace_back(Point(19.75, 6.625, 2.875));
+
}
diff --git a/src/electrophysiology/solvers/ElphySolver.cpp b/src/electrophysiology/solvers/ElphySolver.cpp
index 3bb8e60690984b96a62be94bcc622968ce187c8c..85d0dcef1b6c67a3b07400b4b4695c2b3bcd173a 100644
--- a/src/electrophysiology/solvers/ElphySolver.cpp
+++ b/src/electrophysiology/solvers/ElphySolver.cpp
@@ -60,11 +60,13 @@ std::unique_ptr<ElphySolver> GetElphySolver(const string &solverName, MCellModel
return std::make_unique<ImplicitSolver>(cellModel);
} else if (solverName == "SemiImplicit" || solverName == "SemiImplicitNodes") {
return std::make_unique<SemiImplicitSolver>(cellModel);
- }
+ }else if (solverName=="SemiImplicitOnCells")
+ return std::make_unique<SemiImplicitSolverOnCells>(cellModel);
}
THROW(solverName + " No Solver for this Model implemented.")
}
+
std::unique_ptr<ElphySolver>
GetElphySolver(MCellModel &cellModel, const Vector &excitationData) {
std::string solverName{"SemiImplicit"};
diff --git a/src/electrophysiology/solvers/ElphySolver.hpp b/src/electrophysiology/solvers/ElphySolver.hpp
index a8e6edc5de64acbfd751f0642fb9f47977111877..027b30be49e01ae2a3d771f8e86667f95300e03a 100644
--- a/src/electrophysiology/solvers/ElphySolver.hpp
+++ b/src/electrophysiology/solvers/ElphySolver.hpp
@@ -37,15 +37,16 @@ public:
*/
virtual void Step(IElphyAssemble &A, Vectors &values) = 0;
- virtual void Step(IElphyAssemble &A, Vectors &values, Vector &iota_c, Vector& gamma_f_c) { Step(A,values); }
+ virtual void Step(IElphyAssemble &A, Vectors &values, Vector &iota_c, Vector& gamma_f_c) {
+ Exit("This step is only defined for semi-implicit on cells")
+ //Step(A,values);
+ };
virtual void Finalize(Vectors &values) {};
virtual void GatingVectorAti( Vector &values, int i) const {Exit("Gating vector not defined here");};
};
-
std::unique_ptr<ElphySolver>
GetElphySolver(const std::string &solverName, MCellModel &cellModel, const Vector &excitationData);
-
std::unique_ptr<ElphySolver> GetElphySolver(MCellModel &cellModel, const Vector &excitationData);
diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.cpp b/src/electrophysiology/solvers/LinearImplicitSolver.cpp
index 4360316161a3edf3ef592fad6e46951d3c1be6ca..9c183da7ec5b5f2099a7fc116376d2439b67c36c 100644
--- a/src/electrophysiology/solvers/LinearImplicitSolver.cpp
+++ b/src/electrophysiology/solvers/LinearImplicitSolver.cpp
@@ -22,9 +22,24 @@ void LinearImplicitSolver::Initialize(IElphyAssemble &A, Vector &V) {
}
void LinearImplicitSolver::Initialize(IElphyAssemble &A, Vector &V, Vector &gamma_f_c) {
- Initialize(A,V);
- gating_c = std::make_unique<Vectors>(M.Size(), gamma_f_c);
+ //Initialize(A,V);
+ potential=std::make_unique<Vector>(V);
+ (*potential) = M.InitialValue(ELPHY);
+ if ((gating_c) == nullptr)
+ {
+ const Meshes &mesh = (*potential).GetMeshes();
+ disc0=std::make_unique<LagrangeDiscretization >(mesh, 0, 1);
+ cellSizeV=std::make_unique<Vector>(*disc0);
+ }
+ gating_c = std::make_unique<Vectors>(M.Size(), *cellSizeV);
+ //gating_c = std::make_unique<Vectors>(M.Size(), gamma_f_c);
M.Initialize(*gating_c);
+ K = std::make_unique<Matrix>(V);
+ b = std::make_unique<Vector>(V);
+ if (assembleOnce) {
+ MK = std::make_unique<Matrix>(V);
+ }
+ A.SetInitialValue(V);
}
void LinearImplicitSolver::Method(IElphyAssemble &A, Vector &V) {
@@ -35,11 +50,12 @@ void LinearImplicitSolver::Method(IElphyAssemble &A, Vector &V) {
vNew[2] = 1.0;
}
A.PrintIteration(vNew[0]);
+ if(plotCa){
+ A.PlotCa((*gating)[1]);
+ }
while (!A.IsFinished()) {
A.PlotIteration(vNew[0]);
-
- vout(1) << "# Linear implicit step " << A.Step() + 1 << " from " << A.Time() << " to "
- << A.NextTimeStep(false) << " (" << A.StepSize() << ")" << endl;
+ printSolverNameInStep(A);
vout(10).StartBlock("LinearImplicit Step");
Step(A, vNew);
@@ -66,19 +82,6 @@ void LinearImplicitSolver::Step(IElphyAssemble &A, Vectors &values) {
A.updateActivationTime(values[0], A.Time());
}
-void LinearImplicitSolver::Step(IElphyAssemble &A, Vectors &values, Vector &iota_c, Vector &gamma_f_c) {
- A.Initialize(values[0]);
- A.updateExternalCurrent(values[0]);
- A.updateIionVecs(*gating);
-
- StaggeredStepCell(A, values, iota_c, gamma_f_c);
-
- updateMaxMin();
- updateValues(values,gamma_f_c);
-
- A.NextTimeStep();
- A.updateActivationTime(values[0], A.Time());
-}
void LinearImplicitSolver::SolveGating(double t, double dt) {
Point P(0.0, 0.0, 0.0);
@@ -100,25 +103,23 @@ void LinearImplicitSolver::SolveGating(double t, double dt) {
}
}
}
- if (linearImplicitSolverCell) {
- for (row r = (*gating_c)[0].rows(); r != (*gating_c)[0].rows_end(); ++r) {
- for (int j = 0; j < r.n(); ++j) {
- for (int i = 0; i < M.Size(); ++i) {
- vcw[i] = (*gating_c)[i](r, j);
- }
-
- M.ExpIntegratorUpdate(dt * M.TimeScale(), vcw, vcw);
- for (int i = position; i < M.ElphySize(); ++i) {
- (*gating_c)[i](r, j) = vcw[i];
- if (vcw[i] < 0.0 || vcw[i] > 1.0) {
- }
- }
+}
+void LinearImplicitSolver::SolveGatingOnCells(double t,double dt){
+ std::vector<double> vcw(M.Size());
+ int position = M.GatingIndex();
+ for (row r = (*gating_c)[0].rows(); r != (*gating_c)[0].rows_end(); ++r) {
+ for (int j = 0; j < r.n(); ++j) {
+ for (int i = 0; i < M.Size(); ++i) {
+ vcw[i] = (*gating_c)[i](r(), j);
+ }
+ M.ExpIntegratorUpdate(dt * M.TimeScale(), vcw, vcw);
+ for (int i = position; i < M.ElphySize(); ++i) {
+ (*gating_c)[i](r, j) = vcw[i];
}
}
}
}
-
void LinearImplicitSolver::SolvePDE(IElphyAssemble &A) {
if (assembleOnce) {
if (!assembled) {
@@ -149,7 +150,6 @@ void LinearImplicitSolver::SolvePDE(IElphyAssemble &A) {
}
(*gating)[0] = (*S) * (*b);
-
}
void LinearImplicitSolver::updateValues(Vectors &values) {
@@ -159,11 +159,7 @@ void LinearImplicitSolver::updateValues(Vectors &values) {
}
}
-void LinearImplicitSolver::updateValues(Vectors &values, Vector &gamma_f_c) {
- updateValues(values);
- if (M.Type() == TENSION)
- gamma_f_c = (*gating_c)[M.ElphySize()];
-}
+
void LinearImplicitSolver::updateMaxMin() {
for (row r = (*maxgating).rows(); r != (*maxgating).rows_end(); ++r) {
@@ -178,7 +174,10 @@ void LinearImplicitSolver::updateMaxMin() {
}
}
}
-
+void LinearImplicitSolver::printSolverNameInStep(IElphyAssemble &A){
+ vout(1) << "# Linear implicit step " << A.Step() + 1 << " from " << A.Time() << " to "
+ << A.NextTimeStep(false) << " (" << A.StepSize() << ")" << endl;
+}
void LinearImplicitSolver::PrintMaxMinGating() {
std::vector<double> maxGatingOverVertices(M.Size());
std::vector<double> minGatingOverVertices(M.Size());
@@ -213,117 +212,61 @@ void LinearImplicitSolver::PrintMaxMinGating() {
vout(5) << endl;
}
-
-void ImplicitSolver::Method(IElphyAssemble &A, Vector &V) {
- Vectors vNew(M.Type() == TENSION ? 3 : 1, V);
- Initialize(A, vNew[0]);
- if (M.Type() == TENSION) {
- vNew[1] = 0.0;
- vNew[2] = 1.0;
- }
-
- double t = A.FirstTStep();
-
- A.PlotIteration(V);
- A.PrintIteration(V);
-
- while (!A.IsFinished()) {
- vout(1) << "# Implicit step " << A.Step() << " from " << A.Time() << " to "
- << A.NextTimeStep(false) << " (" << A.StepSize() << ")" << endl;
- vout(10).StartBlock("Implicit step");
-
- Step(A, vNew);
-
- vout(10).EndBlock();
-
-
- A.PlotIteration(vNew[0]);
- A.PrintIteration(vNew[0]);
+void LinearImplicitSolver::SolveConcentration(double t, double dt, const Vectors &values) {
+ std::vector<double> vcw(M.Size() + (M.Type() == TENSION));
+ std::vector<vector<double>> G(ionSteps);
+ for (int i = 0; i < ionSteps; ++i) {
+ G[i] = std::vector<double>(M.Size());
}
- V = vNew[0];
- PrintMaxMinGating();
-
-}
-
-void ImplicitSolver::Step(IElphyAssemble &A, Vectors &values) {
- A.NextTimeStep();
- A.Initialize(values[0]);
- A.updateExternalCurrent(values[0]);
+ for (row r = (*gating)[0].rows(); r != (*gating)[0].rows_end(); ++r) {
+ for (int j = 0; j < r.n(); ++j) {
+ for (int i = 0; i < M.Size(); ++i) {
+ vcw[i] = (*gating)[i](r, j);
+ }
+ if (M.Type() == TENSION) {
+ vcw[M.Size()] = values[2](r, j);
+ }
+ SolveConcentrationStep(M, t, M.TimeScale() * dt, vcw, G);
+ for (int i: ConcentrationIndices(M)) {
+ (*gating)[i](r, j) = vcw[i];
+ }
- SolveGating(A.Time(), A.StepSize());
- SolveConcentration(A.Time(), A.StepSize(), values);
- A.updateVCW(*gating);
- vout(10).StartBlock("Solve PDE");
- bool conv = SolvePDE(A, values);
- vout(10).EndBlock();
- if (!conv) {
- THROW("Newton method did not converge")
+ if (vcw[1] < 0.0) {
+ mout << "concentration of Ca is less than zero " << vcw[1] << endl;
+ }
+ }
}
- updateMaxMin();
- updateValues(values);
-
- A.updateActivationTime(values[0], A.Time());
-}
-
-bool ImplicitSolver::SolvePDE(IElphyAssemble &A, Vectors &values) {
- newton(A, values[0]);
- return newton.converged();
}
+void LinearImplicitSolver::SolveConcentrationOnCells(IElphyAssemble &A, const Vectors &values, Vector &iota_c_pot) {
+ double t = A.Time();
+ double dt = A.StepSize();
-void SemiImplicitSolver::Method(IElphyAssemble &A, Vector &V) {
- Vectors vNew(M.Type() == TENSION ? 3 : 1, V);
- Initialize(A, vNew[0]);
- if (M.Type() == TENSION) {
- vNew[1] = 0.0;
- vNew[2] = 1.0;
- }
- A.PrintIteration(vNew[0]);
- if(plotCa){
- A.PlotCa((*gating)[1]);
- }
- while (!A.IsFinished()) {
- A.PlotIteration(vNew[0]);
-
-
- int step = A.Step() + 1;
- vout(1) << "# Semi implicit step " << A.Step() + 1 << " from " << A.Time() << " to "
- << A.NextTimeStep(false) << " (" << A.StepSize() << ")" << endl;
- vout(10).StartBlock("Semi Implicit Step");
-
- Step(A, vNew);
- A.PrintIteration(vNew[0]);
- vout(10).EndBlock();
- }
- V = vNew[0];
- PrintMaxMinGating();
- A.PlotIteration(V);
-
-
-}
-
-void SemiImplicitSolver::SolvePDE(IElphyAssemble &A) {
- if (!assembled || contains(A.GetElphyProblem().Name(),"Deformed")) {
- A.MPlusK((*gating)[0], *K);
- (*S)(*K);
- assembled = true;
+ std::vector<double> vcw(M.Size() + (M.Type() == TENSION));
+ std::vector<vector<double>> G(ionSteps);
+ for (int i = 0; i < ionSteps; ++i) {
+ G[i] = std::vector<double>(M.Size());
}
- vout(10).StartBlock("AssembleRHS ");
- A.RHS((*gating), *b);
- vout(10).EndBlock();
- (*b) -= (*K) * (*gating)[0];
- (*gating)[0] += (*S) * (*b);
-
- if (linearImplicitSolverCell) {
- for (cell c = (*gating)[0].cells(); c != (*gating)[0].cells_end(); ++c) {
- ScalarElement E((*gating)[0],*c);
- (*gating_c)[0] = E.Value(c.LocalCenter(),(*gating)[0]);
+ for (cell c = (*gating_c)[0].cells(); c != (*gating_c)[0].cells_end(); ++c) {
+ if (values[0].find_cell(c()) == values[0].cells_end()) continue;
+ double rate = 1 / A.GetElphyProblem().Jacobian_rate(*c);
+ row r = (*gating_c)[0].find_row(c());
+ for (int j = 0; j < r.n(); ++j) {
+ for (int i = 0; i < M.Size(); ++i) {
+ vcw[i] = (*gating_c)[i](r, j);
+ }
+ if (M.Type() == TENSION) {
+ vcw[M.Size()] = iota_c_pot(r, j);
+ }
+ SolveConcentrationStep(M, t, M.TimeScale() * dt, vcw, G);
+ for (int i: ConcentrationIndices(M)) {
+ (*gating_c)[i](r, j) = vcw[i];
+ if (deformConcentrationStretch) (*gating_c)[i](r, j) *= rate;
+ }
}
}
}
-
-
void LinearImplicitSolver::selectOrder(std::string o) {
if (o == "Vcw") {
StaggeredStep = [this](IElphyAssemble &A, const Vectors &values) {
@@ -361,16 +304,18 @@ void LinearImplicitSolver::selectOrder(std::string o) {
SolveConcentration(A.Time(), A.StepSize(), values);
SolvePDE(A);
if (plotCa) A.PlotCa((*gating)[1]);
+
};
- StaggeredStepCell = [this](IElphyAssemble &A, const Vectors &values, Vector &iota_c, Vector &gamma_f_c) {
- SolveGating(A.Time(), A.StepSize());
- SolveConcentration(A, values, iota_c, gamma_f_c);
- SolvePDE(A);
- if (plotCa) A.PlotCa((*gating)[1]);
- };
}
};
+void LinearImplicitSolver::StaggeredStepOnCells(IElphyAssemble &A, const Vectors &values,Vector &iota_c){
+ SolveGatingOnCells(A.Time(), A.StepSize());
+ SolveConcentrationOnCells(A,values, iota_c);
+ A.updateIionVecs(*gating_c);
+ SolvePDEOnCells(A);
+ if (plotCa) A.PlotCa((*gating_c)[1]);
+}
void LinearImplicitSolver::GatingVectorAti( Vector &values, int i)const{
values=(*gating)[i];
@@ -392,78 +337,192 @@ void LinearImplicitSolver::selectIonScheme(const string &solvingIonScheme) {
}
}
-void LinearImplicitSolver::SolveConcentration(double t, double dt, const Vectors &values) {
- std::vector<double> vcw(M.Size() + (M.Type() == TENSION));
- std::vector<vector<double>> G(ionSteps);
- for (int i = 0; i < ionSteps; ++i) {
- G[i] = std::vector<double>(M.Size());
+void SemiImplicitSolver::printSolverNameInStep(IElphyAssemble &A){
+ vout(1) << "# Semi implicit step " << A.Step() + 1 << " from " << A.Time() << " to "
+ << A.NextTimeStep(false) << " (" << A.StepSize() << ")" << endl;
+}
+
+void SemiImplicitSolver::SolvePDE(IElphyAssemble &A) {
+ if (!assembled || contains(A.GetElphyProblem().Name(),"Deformed")) {
+ A.MPlusK((*gating)[0], *K);
+ (*S)(*K);
+ assembled = true;
}
- for (row r = (*gating)[0].rows(); r != (*gating)[0].rows_end(); ++r) {
- for (int j = 0; j < r.n(); ++j) {
- for (int i = 0; i < M.Size(); ++i) {
- vcw[i] = (*gating)[i](r, j);
- }
- if (M.Type() == TENSION) {
- vcw[M.Size()] = values[2](r, j);
- }
- SolveConcentrationStep(M, t, M.TimeScale() * dt, vcw, G);
- for (int i: ConcentrationIndices(M)) {
- (*gating)[i](r, j) = vcw[i];
- }
+ vout(10).StartBlock("AssembleRHS ");
+ A.RHS((*gating), *b);
+ vout(10).EndBlock();
+ (*b) -= (*K) * (*gating)[0];
+ (*gating)[0] += (*S) * (*b);
- if (vcw[1] < 0.0) {
- mout << "concentration of Ca is less than zero " << vcw[1] << endl;
+}
+
+void SemiImplicitSolverOnCells:: Initialize(IElphyAssemble &A, Vector &V){
+ potential=std::make_unique<Vector>(V);
+ (*potential) = M.InitialValue(ELPHY);
+ if ((gating_c) == nullptr)
+ {
+ const Meshes &mesh = (*potential).GetMeshes();
+ disc0=std::make_unique<LagrangeDiscretization >(mesh, 0, 1);
+ cellSizeV=std::make_unique<Vector>(*disc0);
+ }
+ gating_c = std::make_unique<Vectors>(M.Size(), *cellSizeV);
+ M.Initialize(*gating_c);
+ K = std::make_unique<Matrix>(V);
+ b = std::make_unique<Vector>(V);
+
+ if (assembleOnce) {
+ MK = std::make_unique<Matrix>(V);
+ }
+
+ A.SetInitialValue(V);
+}
+void SemiImplicitSolverOnCells::printSolverNameInStep(IElphyAssemble &A) {
+ vout(1) << "# Semi implicit on cells step " << A.Step() + 1 << " from " << A.Time() << " to "
+ << A.NextTimeStep(false) << " (" << A.StepSize() << ")" << endl;
+}
+void SemiImplicitSolverOnCells::updateValues(Vectors &values) {
+ values[0] = *potential;
+ if (M.Type() == TENSION) {
+ THROW("This method should not be used!!!");
+ }
+}
+
+void SemiImplicitSolverOnCells::updateValues(Vectors &values, Vector &gamma_f_c_pot) {
+ values[0] = *potential;
+ if (M.Type() == TENSION){
+ 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_pot(r(),j) = (*gating_c)[M.ElphySize()](r(),j);
+ }
}
+
}
+}
+void SemiImplicitSolverOnCells::Method(IElphyAssemble &A, Vector &V) {
+ Vectors vNew(M.Type() == TENSION ? 3 : 1, V);
+ Initialize(A, vNew[0]);
+ if (M.Type() == TENSION) {
+ vNew[1] = 0.0;
+ vNew[2] = 1.0;
+ }
+ A.PrintIteration(vNew[0]);
+ if(plotCa){
+ A.PlotCa((*gating_c)[1]);
}
+ while (!A.IsFinished()) {
+ A.PlotIteration(vNew[0]);
+ printSolverNameInStep(A);
+ vout(10).StartBlock("semi implicit on cells Step");
+ Step(A, vNew);
+ A.PrintIteration(vNew[0]);
+ vout(10).EndBlock();
+
+ }
+ V = vNew[0];
+ A.PlotIteration(V);
}
+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]);
+ A.updateIionVecs(*gating_c);
+ SolvePDEOnCells(A);
+ if (plotCa) A.PlotCa((*gating_c)[1]);
+ updateValues(values);
+ A.NextTimeStep();
+ A.updateActivationTime(values[0], A.Time());
+}
+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);
-void LinearImplicitSolver::SolveConcentration(IElphyAssemble &A, const Vectors &values, Vector &iota_c, Vector &gamma_f_c) {
- double t = A.Time();
- double dt = A.StepSize();
- SolveConcentration(t,dt,values);
-
- std::vector<double> vcw(M.Size() + (M.Type() == TENSION));
- std::vector<vector<double>> G(ionSteps);
- for (int i = 0; i < ionSteps; ++i) {
- G[i] = std::vector<double>(M.Size());
+ updateValues(values,gamma_f_c_pot);
+
+ A.NextTimeStep();
+ A.updateActivationTime(values[0], A.Time());
+}
+void SemiImplicitSolverOnCells::SolvePDEOnCells(IElphyAssemble &A){
+ if (!assembled || contains(A.GetElphyProblem().Name(),"Deformed")) {
+ A.MPlusK((*potential), *K);
+ (*S)(*K);
+ assembled = true;
}
-//TODO: Remove (g_min and g_max are for debugging!)
- double g_min = 1;
- double g_max = 1;
-
- for (cell c = (*gating_c)[0].cells(); c != (*gating_c)[0].cells_end(); ++c) {
- if (values[0].find_cell(c()) == values[0].cells_end()) continue;
- double rate = 1 / A.GetElphyProblem().Jacobian_rate(*c);
- g_min = std::min(rate,g_min);
- g_max = std::max(rate,g_max);
+ A.RHS(*potential, *b);
-
- row r = (*gating_c)[0].find_row(c());
- for (int j = 0; j < r.n(); ++j) {
- for (int i = 0; i < M.Size(); ++i) {
- vcw[i] = (*gating_c)[i](r, j);
- }
- if (M.Type() == TENSION) {
- vcw[M.Size()] = iota_c(r, j);
- }
- SolveConcentrationStep(M, t, M.TimeScale() * dt, vcw, G);
- for (int i: ConcentrationIndices(M)) {
- (*gating_c)[i](r, j) = vcw[i];
- // if (deformConcentrationStretch && (i<6)) (*gating_c)[i](r, j) *= rate;
- }
- }
+
+ (*b) -= (*K) * (*potential);
+ (*potential) += (*S) * (*b);
+
+ for (cell c = (*potential).cells(); c != (*potential).cells_end(); ++c) {
+ ScalarElement E((*potential),*c);
+ (*gating_c)[0](c(),0) = E.Value(c.LocalCenter(),(*potential));
+ }
+}
+void SemiImplicitSolverOnCells::SolvePDE(IElphyAssemble &A){
+ SolvePDEOnCells(A);
+}
+
+
+void ImplicitSolver::Method(IElphyAssemble &A, Vector &V) {
+ Vectors vNew(M.Type() == TENSION ? 3 : 1, V);
+ Initialize(A, vNew[0]);
+ if (M.Type() == TENSION) {
+ vNew[1] = 0.0;
+ vNew[2] = 1.0;
+ }
+
+ double t = A.FirstTStep();
+
+ A.PlotIteration(V);
+ A.PrintIteration(V);
+
+ while (!A.IsFinished()) {
+ vout(1) << "# Implicit step " << A.Step() << " from " << A.Time() << " to "
+ << A.NextTimeStep(false) << " (" << A.StepSize() << ")" << endl;
+ vout(10).StartBlock("Implicit step");
+
+ Step(A, vNew);
+
+ vout(10).EndBlock();
+
+
+ A.PlotIteration(vNew[0]);
+ A.PrintIteration(vNew[0]);
}
+ V = vNew[0];
+ PrintMaxMinGating();
+
+}
+
+void ImplicitSolver::Step(IElphyAssemble &A, Vectors &values) {
+ A.NextTimeStep();
+ A.Initialize(values[0]);
+ A.updateExternalCurrent(values[0]);
- g_min = PPM->Min(g_min);
- g_max = PPM->Min(g_max);
+ SolveGating(A.Time(), A.StepSize());
+ SolveConcentration(A.Time(), A.StepSize(), values);
+ A.updateVCW(*gating);
+ vout(10).StartBlock("Solve PDE");
+ bool conv = SolvePDE(A, values);
+ vout(10).EndBlock();
+ if (!conv) {
+ THROW("Newton method did not converge")
+ }
+ updateMaxMin();
+ updateValues(values);
- mout << "J rate in [" << g_min << "," << g_max << "]" << endl;
+ A.updateActivationTime(values[0], A.Time());
+}
-
+bool ImplicitSolver::SolvePDE(IElphyAssemble &A, Vectors &values) {
+ newton(A, values[0]);
+ return newton.converged();
}
+
void QuadraticIonStep(MCellModel &cellModel, double t, double dt, std::vector<double> &vcw,
std::vector<std::vector<double>> &g) {
cellModel.ConcentrationUpdate(vcw, g[0], cellModel.IExt(t));
diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.hpp b/src/electrophysiology/solvers/LinearImplicitSolver.hpp
index 866c3cc997a53cf1965a88b9d434672a1b9a4efa..ef81917b06c9bc2ba889be4ff0dfd266b1536757 100644
--- a/src/electrophysiology/solvers/LinearImplicitSolver.hpp
+++ b/src/electrophysiology/solvers/LinearImplicitSolver.hpp
@@ -31,6 +31,7 @@ protected:
std::unique_ptr<Matrix> K;
std::unique_ptr<Vector> b;
std::unique_ptr<Matrix> MK;
+ std::unique_ptr<Vector> potential;
std::unique_ptr<Vectors> gating;
std::unique_ptr<Vectors> gating_c;
std::unique_ptr<Vectors> maxgating;
@@ -47,6 +48,8 @@ protected:
int ionSteps{1};
bool plotCa{false};
bool linearImplicitSolverCell = false;
+ std::unique_ptr<Vector> cellSizeV;
+ std::unique_ptr<LagrangeDiscretization> disc0;
void selectIonScheme(const string &solvingIonScheme);
@@ -92,30 +95,31 @@ public:
void SolveConcentration(double t, double dt, const Vectors &values);
- void SolveConcentration(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);
virtual void SolvePDE(IElphyAssemble &A);
+ virtual void SolvePDEOnCells(IElphyAssemble &A){
+ Exit("only defined for semi-implicit solver");
+ };
void Step(IElphyAssemble &A, Vectors &values) override;
- void Step(IElphyAssemble &A, Vectors &values, Vector &iota_c, Vector &gamma_f_c) override;
-
std::function<void(IElphyAssemble &A, const Vectors &values)> StaggeredStep;
-
- std::function<void(IElphyAssemble &A, const Vectors &values, Vector &iota_c, Vector &gamma_f_c)> StaggeredStepCell;
+ void StaggeredStepOnCells(IElphyAssemble &A, const Vectors &values,Vector &iota_c);
void selectOrder(std::string o);
void updateMaxMin();
void PrintMaxMinGating();
+ virtual void printSolverNameInStep(IElphyAssemble &A);
+ virtual void updateValues(Vectors &values);
- void updateValues(Vectors &values);
-
- void updateValues(Vectors &values, Vector &gamma_f_c);
+ virtual void updateValues(Vectors &values, Vector &gamma_f_c){Exit("only defined for semi-implicit solver");};
void GatingVectorAti( Vector &values, int i) const override;
};
@@ -145,9 +149,30 @@ class SemiImplicitSolver : public LinearImplicitSolver {
public:
SemiImplicitSolver(MCellModel &cellModel) : LinearImplicitSolver(cellModel) {};
- void Method(IElphyAssemble &A, Vector &V) override;
+ void SolvePDE(IElphyAssemble &A) override;
+ void printSolverNameInStep(IElphyAssemble &A) override;
+};
+class SemiImplicitSolverOnCells : public LinearImplicitSolver {
+ //std::unique_ptr<Vector> cellSizeV;
+ //std::unique_ptr<LagrangeDiscretization> disc0;
+public:
+ //SemiImplicitSolverOnCells(MCellModel &cellModel,const Vector &cSV) : LinearImplicitSolver(cellModel) {
+ // cellSizeV=std::make_unique<Vector>(cSV);
+ //};
+ SemiImplicitSolverOnCells(MCellModel &cellModel) : LinearImplicitSolver(cellModel){};
+
+
+ void Initialize(IElphyAssemble &A, Vector &potential) override;
+ void printSolverNameInStep(IElphyAssemble &A) override;
+ void Method(IElphyAssemble &A, Vector &V) override;
+ void Step(IElphyAssemble &A, Vectors &values, Vector &iota_c, Vector &gamma_f_c) override;
+ void Step(IElphyAssemble &A, Vectors &values) override;
void SolvePDE(IElphyAssemble &A) override;
+ void SolvePDEOnCells(IElphyAssemble &A) override;
+ void updateValues(Vectors &values) override;
+ void updateValues(Vectors &values, Vector &gamma_f_c) override;
+
};
#endif //ELPHYSOLVERS_H
\ No newline at end of file
diff --git a/tools/evaluation/simpleTests.py b/tools/evaluation/simpleTests.py
new file mode 100644
index 0000000000000000000000000000000000000000..7fa00ca344b0ebe11ac84ba96dcc52711efd19fb
--- /dev/null
+++ b/tools/evaluation/simpleTests.py
@@ -0,0 +1,16 @@
+from utility import reading as read
+import matplotlib.pyplot as plt
+
+# def plotPotential(self,dt,T,V)
+def plotPotential(filename ='../../build/log/logfile.log'):
+ print('plot from file', filename)
+ (T,dt)=read.getTimeDataFromLog(filename)
+ read.extractDataFromLog('../../build/','logfile.log',False)
+ (t,V)=read.getDataFromTXT('../../build/log/','logfile',1)
+ plt.plot(t,V)
+ (t,V)=read.getDataFromTXT('../../build/log/','logfile',9)
+ plt.plot(t,V)
+ plt.show()
+
+if __name__=="__main__":
+ plotPotential()
diff --git a/tools/evaluation/thesisEvaluations.py b/tools/evaluation/thesisEvaluations.py
index c10b3a5820098bb4a6faf6bbacfce1396b756bf6..329ff192a94f83a44cdf54ef00d7ace9331f31b1 100644
--- a/tools/evaluation/thesisEvaluations.py
+++ b/tools/evaluation/thesisEvaluations.py
@@ -17,8 +17,8 @@ def Histplot(Path,language):
def MeshInfoTab(path,savepath,tablename,language):
- lL=[0,1,2]
- tL=[0,1,2]
+ lL=[0,1,2,3]
+ tL=[0,1,2,3]
sdt=0.0004
endTime=0.1
SI=path
@@ -67,29 +67,33 @@ def plotPotential(path):
def writeConvergenceBiventricleActivationTime(language='eng'):
path='../../../data/biventricleTT/SI0004/'
SI=path
- lL=[0,1,2]
- tL=[0,1,2]
+ OnCellsSI ='../../../data/biventricleTT/OnCells/SI0004/'
+ lL=[0,1,2,3]
+ tL=[0,1,2,3]
sdt=0.0004
endTime=0.3
listEvaluationPoints=[]
- l0Path=path+'vtu/AT_l0j0.vtu'
+ saveInDiss=True
- #read.writeActivationDataFiles(path,lL,tL,True)
+ #read.writeActivationDataFiles(OnCellsSI,lL,tL,True)
-
- plot=Plot(lL,tL,sdt,endTime,1,listEvaluationPoints)
- ValueArray=plot.MeanTactSpace(SI,l0Path,[],language)
- print('Plot finished')
- savePath='../../data/biventricleTT/SI0004/tex/'#'images/'
- tabWriter = TableWriter(lL,tL,sdt,endTime,1,[SI],listEvaluationPoints)
- L=0
- J=0
- tabWriter.writeMeanTactTable(SI,l0Path,savePath,L,J,SI,ValueArray,language)
+ experimentslist=[SI,OnCellsSI]
+ nameList=['SIVert','SICells']
+ for i in range(len(experimentslist)):
+ l0Path=experimentslist[i]+'vtu/AT_l0j0.vtu'
+ plot=Plot(lL,tL,sdt,endTime,1,listEvaluationPoints)
+ ValueArray=plot.MeanTactSpace(experimentslist[i],l0Path,nameList[i],[],saveInDiss,language)
+ print('Plot finished')
+ savePath='images/'
+ tabWriter = TableWriter(lL,tL,sdt,endTime,1,[experimentslist[i]],listEvaluationPoints)
+ L=0
+ J=0
+ tabWriter.writeMeanTactTable(experimentslist[i],l0Path,savePath,L,J,experimentslist[i],ValueArray,nameList[i],saveInDiss,language)
if __name__=="__main__":
language='d'
- #MeshInfoTab(pathTT+'SI0004/','../../../thesis-lindner/thesis/tables/','MeshInfoSpaceBiventricle',language)
+ MeshInfoTab(pathTT+'SI0004/','../../../thesis-lindner/thesis/tables/','MeshInfoSpaceBiventricle',language)
#Histplot(pathTT,language)
- writeConvergenceBiventricleActivationTime(language)
+ #writeConvergenceBiventricleActivationTime(language)
diff --git a/tools/evaluation/utility/latexTables.py b/tools/evaluation/utility/latexTables.py
index 51ddb97334f850009f1904f0d8f5cd85d215128b..bd18d35b13fee6ce287283854103cc34560c8c2e 100644
--- a/tools/evaluation/utility/latexTables.py
+++ b/tools/evaluation/utility/latexTables.py
@@ -1308,7 +1308,7 @@ class Table:
for l in self.lL:
path = eP+ 'vtu/AT_l'+str(l)+'j'+str(j)+'.vtu'
tact=-10
- if len(varr)!=0:
+ if l<len(varr) and j<len(varr[l]):
tact=varr[l][j]
else:
tact=comp.computeMeanTact(path,l0Path)
@@ -1869,22 +1869,25 @@ class TableWriter:
output_file.write(self.tab.EndTabular())
output_file.write(self.tab.writeEndTable())
- def writeMeanTactTable(self,fn,l0Path,sP,L,J,refPath, vArr,language='eng'):
+ def writeMeanTactTable(self,fn,l0Path,sP,L,J,refPath, vArr,saveName='',saveInDiss=False,language='eng'):
if language=='eng':
cap='The activation time $\\tact(\V)$ and the difference with respect to the reference solution on different levels in space and time for the semi-implicit (\\SISVI) scheme (the numbers are displayed in ms)'
elif language=='d':
cap='Die Aktivierungszeit $\\tact(\V)$ und die Differenz zur Referenzlösung $\\rmse(\V^{j,l},\V^{'+str(J)+','+str(L)+'})$ für verschiedene Zeit- und Ortsdiskretisierungen mit dem semi-impliciten Verfahren(\\SISVI)'
else:
print('language not known')
- label='tab:ActiAndError'
- filename=fn+'tex/'+'ActivationVentricles'
+ label='tab:ActiAndError'+saveName
+ if saveInDiss:
+ filename='../../../thesis-lindner/thesis/figure/'+'ActivationVentricles'+saveName
+ else:
+ filename=fn+'tex/'+'ActivationVentricles'+saveName
with open(filename+'.tex', 'w') as output_file:
output_file.write(self.tab.writeBeginTable('h'))
output_file.write(self.tab.writeCaption(cap,label))
#output_file.write('\\hspace*{-2mm}\n')
output_file.write('\\vspace*{7mm}\n')
output_file.write('\\begin{tabular}{cc}\n')
- output_file.write('\\resizebox{9.3cm}{!}{\\input{'+sP+'ActivationSpace}}\n')
+ output_file.write('\\resizebox{9.3cm}{!}{\\input{'+sP+'ActivationSpace'+saveName+'}}\n')
output_file.write('\\\\[-6.2cm]\n')
output_file.write('&\n')
output_file.write(self.tab.writeMeanTactTabular(fn,l0Path,vArr))
diff --git a/tools/evaluation/utility/plotting.py b/tools/evaluation/utility/plotting.py
index 32ca57eda34fa576dab224005dee7e53569f5256..43cbfdefff6d5e029c750e40f347409a4bf25fba 100644
--- a/tools/evaluation/utility/plotting.py
+++ b/tools/evaluation/utility/plotting.py
@@ -45,7 +45,7 @@ class Plot:
self.plotData()
if plotLegend:
plt.legend(bbox_to_anchor=(1.1,1), loc="upper left")
- self.saveTikz(sP,plotname)
+ self.saveTikz(sP+'/tex/', plotname)
plt.close()
def plotSchemesInP(self, p,l,j,L,J,fL,nL,refF,refN,sP,plotLegend=False,resize=1):
for index in range(len(fL)):
@@ -65,7 +65,7 @@ class Plot:
self.plotData()
if plotLegend:
plt.legend(bbox_to_anchor=(1.1,1), loc="upper left")
- self.saveTikz(sP,plotname)
+ self.saveTikz(sP+'/tex/', plotname)
plt.close()
def plotAllPoints(self,fL):
for file in fL:
@@ -81,7 +81,7 @@ class Plot:
plotname ='BRl' +str(l)+ 'j'+str(j)+'m'+str(m)+'AllPoints'
self.plotData(True)
plt.legend(bbox_to_anchor=(0.95,0.05), loc="lower right")
- self.saveTikz(path,plotname)
+ self.saveTikz(path+'/tex/', plotname)
plt.close()
def plotSpaceForPointp(self,path,p,j,m,n):
@@ -96,7 +96,7 @@ class Plot:
plotname ='BRSpaceConvergence'+ 'j'+str(j)+'m'+str(m)+'P'+str(p)
self.plotData()
plt.legend(loc='lower right')
- self.saveTikz(path,plotname)
+ self.saveTikz(path+'/tex/', plotname)
plt.close()
def plotTimeForPointp(self,path,p,l,m,resize):
for j in self.tL:
@@ -110,7 +110,7 @@ class Plot:
plotname ='BRTimeConvergence'+ 'l'+str(l)+'m'+str(m)+'P'+str(p)
self.plotData()
plt.legend()
- self.saveTikz(path,plotname)
+ self.saveTikz(path+'/tex/', plotname)
plt.close()
def plotSpaceForPointpWithRef(self,path,p,j,m,L,J,resize,interval=None):
levelList=[]
@@ -144,7 +144,7 @@ class Plot:
#plotname ='BRSpaceConvergenceWithRef'+ 'j'+str(j)+'m'+str(m)+'P'+str(p)
self.plotData(True)
plt.legend()
- self.saveTikz(path,plotname)
+ self.saveTikz(path+'/tex/', plotname)
plt.close()
def plotTimeForPointpWithRef(self,path,p,l,m,L,J,resize,interval=None):
for j in self.tL:
@@ -171,7 +171,7 @@ class Plot:
plotname ='BRTimeConvergenceWithRef1015'+ 'l'+str(l)+'m'+str(m)+'P'+str(p)
self.plotData(True)
plt.legend()
- self.saveTikz(path,plotname)
+ self.saveTikz(path+'/tex/', plotname)
plt.close()
def valuePlot(self,p,filepath,nL):
for j in self.tL:
@@ -329,7 +329,7 @@ class Plot:
self.saveTikz(eP,plotname)
plt.close()
return plotname
- def MeanTactSpace(self,eP,l0Path,vArr=[],language='eng'):
+ def MeanTactSpace(self,eP,l0Path,saveName='',vArr=[],saveInDissFolder=False,language='eng'):
vArray=[]
for l in self.lL:
tactList=[]
@@ -350,8 +350,11 @@ class Plot:
vArray.append(tactList)
self.DataTact(language)
# plt.show()
- plotname = 'ActivationSpace'
- self.saveTikz(eP, plotname)
+ plotname = 'ActivationSpace'+saveName
+ if saveInDissFolder:
+ self.saveTikz('../../../thesis-lindner/thesis/images/', plotname)
+ else:
+ self.saveTikz(eP+'/tex/', plotname)
plt.close()
return vArray
def ActivationTimePlot(self):
@@ -405,7 +408,7 @@ class Plot:
plt.grid(alpha=0.5)
plotname = 'NiedBMInSpace'
- self.saveTikz(path, plotname)
+ self.saveTikz(path+'/tex/', plotname)
#plt.show()
plt.close()
@@ -416,7 +419,7 @@ class Plot:
plt.plot(lengthList, actiList, label=labelname, linewidth=1.5)
self.DiagonalActivationTimeLabel()
plotname='NiedBMInTime'
- self.saveTikz(path, plotname)
+ self.saveTikz(path+'/tex/', plotname)
#plt.show()
plt.close()
@@ -524,7 +527,7 @@ class Plot:
print('time unit not implemented')
def saveTikz(self, p,pname):
import tikzplotlib
- tikzplotlib.save(p+'tex/' + pname + ".tex")
+ tikzplotlib.save(p+ pname + ".tex")
def plotData(self,inMS=False):
xmin, xmax, ymin, ymax = 0.0, self.T+0.05, -100.0,70.0
@@ -605,7 +608,7 @@ class Plot:
plt.plot(time,order,'r--',linewidth=3.5,label ='order 2')
self.plotDataConPlot()
plotname ='ConvergenceTime' + 'l'+ str(l)+'m'+str(self.m)
- self.saveTikz(path,plotname)
+ self.saveTikz(path+'/tex/', plotname)
plt.close()
def convergencePlotSpace(self, path,j,numberofEPwithNorms):
print("Plot Space convergence")
@@ -622,7 +625,7 @@ class Plot:
plt.plot(space,order,'k--',linewidth=3.5, label='order 1')
self.plotDataConPlot()
plotname ='ConvergenceSpace' + 'j'+ str(j)+'m'+str(self.m)
- self.saveTikz(path,plotname)
+ self.saveTikz(path+'/tex/', plotname)
plt.close()
def plotSpaceExtrapolates(self,fL,nL,savePath,interval=None):
for point in self.evalP:
@@ -747,7 +750,7 @@ class Plot:
self.plotData()
plt.legend(bbox_to_anchor=(1.1,1), loc="upper left")
#plt.legend()
- self.saveTikz(sP,pN)
+ self.saveTikz(sP+'/tex/', pN)
plt.close()
def plotTwoPotentials(self,dt,T,V,W):
diff --git a/tools/evaluation/utility/reading.py b/tools/evaluation/utility/reading.py
index 187d7c0e069ca9f8f8812270df5a636e49c947f2..e30bd072d0b9abb8e5650662fd9f831f3f46969f 100644
--- a/tools/evaluation/utility/reading.py
+++ b/tools/evaluation/utility/reading.py
@@ -114,7 +114,20 @@ def getTimeData(path,l_entry,t_entry,m):
dt = m.group(1)
f.close()
return (float(EndTime),float(dt))
-
+def getTimeDataFromLog(logfile):
+ f = open(logfile)
+ EndTime=0.0
+ dt = 0.0
+ for l in f:
+ m = re.search("EndTime:\s*\.*\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)", l)
+ if m:
+ EndTime =m.group(1)
+ m = re.search("DeltaTime:\s*\.*\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)", l)
+ if m:
+ dt = m.group(1)
+ f.close()
+ return (float(EndTime),float(dt))
+
def getPoints(path,l_entry,t_entry,m,number):
f = open(path+ 'log/'+'log_' +createFilename(l_entry,t_entry,m))
eP = [None] * (number-3)
@@ -144,7 +157,7 @@ def writeDataInFile(path,l_List,t_List,m,option=True):
def extractDataFromLog(path,fn,option):
- if fn=='log':
+ if fn=='logfile.log':
f=open(path+ 'log/'+fn)
else:
if option:
@@ -175,7 +188,10 @@ def extractDataFromLog(path,fn,option):
m=re.search("ElphyModel: [\.]+ SemiImplicit",l)
if m:
model='SemiImplicit'
-
+
+ m=re.search("ElphyModel: [\.]+ SemiImplicitOnCells",l)
+ if m:
+ model='SemiImplicitOnCells'
if model=='LinearImplicit':
m = re.search("# Linear implicit step (\d+)\s*from\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)\s*to\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)", l)
if m:
@@ -191,6 +207,12 @@ def extractDataFromLog(path,fn,option):
if m:
steps[m.group(3)] = OrderedDict()
current_step = m.group(3)
+ elif model == 'SemiImplicitOnCells':
+ m=re.search("# Semi implicit on cells step (\d+)\s*from\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)\s*to\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)", l)
+ if m:
+ steps[m.group(3)] = OrderedDict()
+ current_step = m.group(3)
+
else:
m = re.search("Solving Step (\d+)\s*at time\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)\.",l)
#m = re.search("# Euler step (\d+)\s*from\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)\s*to\s*([+\-eE\d]+\.{0,1}[+\-eE\d]*)", l)
@@ -217,9 +239,9 @@ def extractDataFromLog(path,fn,option):
steps[current_step][int(m.group(1))] = m.group(5)
#print(steps[current_step])
- f.close()
- if fn=='log':
- out_f = open(path + 'log/' +fn +'.txt', "w")
+ f.close()
+ if fn=='logfile.log':
+ out_f = open(path + 'log/logfile.txt', "w")
else:
if option:
out_f = open(path + 'data/' +fn +'.txt', "w")
@@ -245,9 +267,7 @@ def extractDataFromLog(path,fn,option):
listVValues.append(v[11])
#listVValues[5],listVValues[6]=listVValues[6],listVValues[5]
-
-
- out_f.write("{} {}\n".format('{:.8f}'.format(float(k)), " ".join(listVValues)))
+ out_f.write("{} {}\n".format('{:.8f}'.format(float(k)), " ".join(listVValues)))
out_f.close()
diff --git a/tools/solving/biventricleThesis.py b/tools/solving/biventricleThesis.py
index 90ebef0a5f3f242d3e1a011be27f20c2a73ed339..c01b6dc38ab8b78a8955da02f0a41338a4706559 100644
--- a/tools/solving/biventricleThesis.py
+++ b/tools/solving/biventricleThesis.py
@@ -4,32 +4,33 @@ import os
import shutil
import mpp.python.mppy as mppy
-lList = [0,1,2,3]
-lList = [0,1,2]
-tList = [0,1,2]
+lList = [3]#,1,2,3]
+tList = [0,1,2]#,1,2,3]
#start_dt=0.0001
start_dt=0.0004
path='../../../data/biventricleTT/'
pathToMakeFolders=path #'../'+path
-def moveActivationtime(l,j, expPath):
- shutil.move('../../build/data/vtu/ActivationTime.vtu', pathToMakeFolders+expPath+'AT_l'+str(l)+'j'+str(j)+'.vtu')
+def moveActivationtime(l,j, expPath,moveVtu=False):
+ shutil.move('../../build/data/AT_l'+str(l)+'j'+str(j)+'.txt', pathToMakeFolders+expPath+'AT_l'+str(l)+'j'+str(j)+'.txt')
+ if moveVtu:
+ shutil.move('../../build/data/vtu/ActivationTime.vtu', pathToMakeFolders+expPath+'AT_l'+str(l)+'j'+str(j)+'.vtu')
-def startExperiment(procs,startconf,pE,added=None):
+def startExperiment(procs,startconf,pE,added=None,moveVTU=False):
checkIfFolderExists(pE)
mpp = mppy.Mpp(project_name='CardMech',executable='Elphy-M++',kernels=4,mute=False)
for l in lList:
for j in tList:
file = path+pE+'log_l'+str(l)+'j'+str(j)+'m1'
- kwargs={"ElphyLevel": l,"DeltaTime":start_dt*2**(-j),"logfile":file,"EndTime":0.3,"ElphyVTK":0}
+ kwargs={"ElphyLevel": l,"DeltaTime":start_dt*2**(-j),"logfile":file,"EndTime":0.3,"ElphyVTK":0,"ElphyTimeLevel":j}
if added!=None:
kwargs.update(added)
mpp.run(procs, config=startconf, kwargs=kwargs) #replace 64 by hosts if necessary
- moveActivationtime(l,j,pE)
+ moveActivationtime(l,j,pE,moveVTU)
def checkIfFolderExists(path):
folderPath =pathToMakeFolders+path
@@ -39,10 +40,14 @@ def checkIfFolderExists(path):
if __name__=="__main__":
procs=64
startconf="electrophysiology/BiVentricleTest/biventricelTest"
-
+ moveVTU=True
pathExperiment="SI0004/"
- startExperiment(procs,startconf,pathExperiment,{"ElphyModelName": "TenTusscher"})
+ startExperiment(procs,startconf,pathExperiment,{"ElphyModelName": "TenTusscher","ParallelPlotting":False,"Compression":"ascii"},moveVTU)
+
+ pathExperiment="OnCells/SI0004/"
+
+ #startExperiment(procs,startconf,pathExperiment,{"ElphyModelName": "TenTusscher","ElphyModel":"SemiImplicitOnCells","ParallelPlotting":False,"Compression":"ascii"},moveVTU)