diff --git a/conf/elasticity/benchmarks/LandBenchmark.conf b/conf/elasticity/benchmarks/LandBenchmark.conf
index 92164b58cae47c59a894bbe466e47f8465b73a73..6efcfcc6fc4dd03062d4c4e8daf430f5948df646 100644
--- a/conf/elasticity/benchmarks/LandBenchmark.conf
+++ b/conf/elasticity/benchmarks/LandBenchmark.conf
@@ -2,28 +2,44 @@
ClearData = true
# ========================================================
# = Location of .geo files ===============================
-GeoPath = ../geo/
+#GeoPath = ../geo/
# ========================================================
h0=1e-4
h_min=1e-6
+
#logfile = AS_T4_1_deg1_gamma07
-activeDeformation = Quarteroni2
-g = 0.75 #gamma for updateStretch
-scale_gamma = -1.7
+# CAUTION: right now there is no possibility to use homotopy for different meshes and/or different degrees
+# enables the possibility to start homotopy for different discretizations or to test penalty choices
+multiDiscs = 0
+# CAUTION: if useDifferentMeshes = 0, it will use the entry in Mesh!!!
+useDifferentMeshes = 0 # possibility to change meshes isn't implemented right now
+MechDiscretizationList = Conforming,DG # chosing different discs for homotopy
+MechPolynomialDegreeList = 1 #if there is only one entry, it usess the same degree for all discs
+MeshList = OrientedEllipsoid_2 #if there is only one entry, it usess the same mesh for all discs
+# If |DGPenaltyList| != |MechDiscretizationList| DGPenalty = 0 für fehlende Einträge
+DGPenaltyList= 0,1000,2000
+
+
+#activeDeformation = Quarteroni2
+#g = 0.75 #gamma for updateStretch
+#scale_gamma = -1.7
#StretchFactorNormal =1
-activeStress = 0#60
+activeStress = 60 #kPa
Model = ActiveStrainElasticity
MechRuntype = Default
MechDynamics = Static
-MechDiscretization = Conforming
+MechDiscretization = Conforming#DG
+
+DGSign=-1.0
+DGPenalty= 1000# Penalty for EG=200
CalculateReference = 0
-ReferenceLevel = 3
-InterpolateStartVector = true
+ReferenceLevel = 5
+InterpolateStartVector = false #true
#Mesh=FullEllipsoid
-Mesh = OrientedEllipsoid_4
+Mesh = OrientedEllipsoid_2
MechProblem = LandEllipsoid
MechMeshPart = leftventricle
@@ -56,14 +72,14 @@ TractionC = 0
#------------------------
Incompressible = false
QuasiCompressiblePenalty = Ciarlet
-VolumetricPenalty = 1000
+VolumetricPenalty = 1000 # 1 MPa
######################
# Pressure Parameters#
######################
#ConstantPressure = -10
-LVPressure = 15
+LVPressure = 15 #kPa
RVPressure = 0
RAPressure = 0
LAPressure = 0
@@ -78,17 +94,17 @@ PlotCellPotential = 0
# === Mechanics Solver ============================================
MechSolver = gmres
-MechPreconditioner = SuperLU;
-MechEpsilon = 1e-8;
-MechReduction = 1e-12;
+MechPreconditioner = GaussSeidel;
+MechEpsilon = 1e-8#1e-8;
+MechReduction = 1e-5;
MechSteps = 10000;
# === Newton Method === #
-NewtonEpsilon = 1e-6
-NewtonReduction = 1e-12;
+NewtonEpsilon = 1e-4
+NewtonReduction = 1e-5;
NewtonSteps = 200
NewtonLineSearchSteps = 10;
-NewtonLinearizationReduction = 1e-12
+NewtonLinearizationReduction = 1e-3
#NewtonJacobiUpdate = 8
#NewtonDamping = 0.75
@@ -137,21 +153,21 @@ DynamicSolverVerbose = -1
CoupledSolverVerbose = 2
CellModelVerbose = -10
-MechVTK = 1
+MechVTK = 0
# = Main Verbose =
ElphyVerbose = -1
MechVerbose = 1
# = Plot Verbose =
-PlotVTK = 1
+PlotVTK = 0
PlottingSteps = 1
ElphySolverVTK = -1
-PressureSolverVTK = 1
+PressureSolverVTK = 0
DynamicSolverVTK = -1
CoupledSolverVTK = -1
PrestressPlotVerbose = 0
-PressurePlotVerbose = 1
+PressurePlotVerbose = 0
DynamicPlotVerbose = -1
diff --git a/conf/elasticity/benchmarks/beam.conf b/conf/elasticity/benchmarks/beam.conf
index a606258f2e41815a832b0d77b554ba5444b16c55..56fdc2c578a929be0e321a21f64486476a43d185 100644
--- a/conf/elasticity/benchmarks/beam.conf
+++ b/conf/elasticity/benchmarks/beam.conf
@@ -7,10 +7,9 @@ h_min=1e-6
Model = ActiveStrainElasticity
MechRuntype = Default
MechDynamics = Static
-#MechDiscretization = DG
-MechDiscretization = Conforming
+MechDiscretization = DG #EG#Conforming
-CalculateReference = 1
+CalculateReference = 0
ReferenceLevel = 1
InterpolateStartVector = false
@@ -23,13 +22,13 @@ ProblemDimension = 3
ProblemGeometry = Tet
# Amount of pressure steps per iteration
-PressureSteps = 1
+PressureSteps = 10
PressureIterations = 10
PressureDepth = 0
DebugPressure = 1
-MechLevel = 0
+MechLevel = 1
MechPolynomialDegree = 1
WithPrestress=false
@@ -181,14 +180,16 @@ ElphyVerbose = -1
MechVerbose = -1
# = Plot Verbose =
-MechPlot=-1
+MechPlot=1
ElphySolverVTK = -1
-PressureSolverVTK = -1
+PressureSolverVTK = 1
DynamicSolverVTK = -1
CoupledSolverVTK = -1
PrestressPlotVerbose = -1
-PressurePlotVerbose = -1
+PressurePlotVerbose = 1
DynamicPlotVerbose = -1
-PlotVTK = -1
+PlotVTK = 1
+
+ParallelPlotting = 0
diff --git a/conf/elasticity/m++.conf b/conf/elasticity/m++.conf
index 97de08d7c08eb6a75a601a12443d6bf4b5d55e86..c100113caf546b637bda1b47e8eea51b50feddbe 100644
--- a/conf/elasticity/m++.conf
+++ b/conf/elasticity/m++.conf
@@ -11,8 +11,8 @@ ClearData = true
#loadconf = elasticity/benchmarks/dg_beam.conf
#loadconf = elasticity/benchmarks/eg_beam.conf
#loadconf = elasticity/benchmarks/ellipsoid.conf
-#loadconf = elasticity/benchmarks/LandBenchmark.conf
-loadconf = elasticity/benchmarks/LandBenchmarkDG.conf
+loadconf = elasticity/benchmarks/LandBenchmark.conf
+#loadconf = elasticity/benchmarks/LandBenchmarkDG.conf
#loadconf = elasticity/benchmarks/LandBenchmarkBiVentricle.conf
#loadconf = elasticity/elasticity.conf
# ======================================
diff --git a/pipeline/.cardmech-benchmark-on-horeka.yml b/pipeline/.cardmech-benchmark-on-horeka.yml
index 94f45f5c82b68a2ed738d25a2c98fc85cb05051e..8182cb668a9f7aea0f450cb4d456896fc3c0b71f 100644
--- a/pipeline/.cardmech-benchmark-on-horeka.yml
+++ b/pipeline/.cardmech-benchmark-on-horeka.yml
@@ -26,3 +26,61 @@ ellipsoid-test:
CMAKE_ARGS: '-DMPP_BUILD_TYPE=MppRelease'
script:
- salloc -p cpuonly -t 00:01:00 -n 4 mpirun Elphy-M++ electrophysiology/TetraTest/start
+
+.land-benchmark_on_horeka:
+ extends: .benchmark_on_horeka
+ only:
+ variables:
+ - $BENCHMARK == $CI_JOB_NAME || $BENCHMARK == 'all' || $BENCHMARK == 'on-horeka' || $BENCHMARK == 'land-on-horeka'
+ variables:
+ CMAKE_ARGS: '-DMPP_BUILD_TYPE=MppRelease'
+
+
+
+# extends: .land-benchmark_on_horeka
+ # script:
+# - salloc -p cpuonly -t 00:05:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=1 CalculateReference=0
+# - salloc -p cpuonly -t 00:15:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=1 MechPolynomialDegree=1 CalculateReference=0
+# - salloc -p cpuonly -t 01:00:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=2 MechPolynomialDegree=1 CalculateReference=0
+# - salloc -p cpuonly -t 15:00:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=3 MechPolynomialDegree=1 CalculateReference=0
+ # - salloc -p cpuonly -t 72:00:00 -n 512 -N 16 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=4 MechPolynomialDegree=1 CalculateReference=0
+# - salloc -p cpuonly -t 100:00:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=5 MechPolynomialDegree=1 CalculateReference=0
+
+# alle auf mehr Kernen testen
+
+#landbenchmark-deg1_ref:
+# extends: .land-benchmark_on_horeka
+# script:
+# - salloc -p cpuonly -t 72:00:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=1 CalculateReference=1 ReferenceLevel=4
+# - salloc -p cpuonly -t 72:00:00 -n 1024 -N 16 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=1 CalculateReference=1 ReferenceLevel=4
+# - salloc -p cpuonly -t 72:00:00 -n 512 -N 8 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=1 CalculateReference=1 ReferenceLevel=4
+# - salloc -p cpuonly -t 72:00:00 -n 512 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=1 CalculateReference=1 ReferenceLevel=4
+# - salloc -p cpuonly -t 72:00:00 -n 1024 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=1 CalculateReference=1 ReferenceLevel=4
+# - salloc -p cpuonly -t 72:00:00 -n 1024 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=1 CalculateReference=1 ReferenceLevel=5
+
+#TODO: Zeiten anpassen
+landbenchmark-deg2:
+ extends: .land-benchmark_on_horeka
+ script:
+# - salloc -p cpuonly -t 00:05:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=2 CalculateReference=0
+# - salloc -p cpuonly -t 00:15:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=1 MechPolynomialDegree=2 CalculateReference=0
+# - salloc -p cpuonly -t 01:00:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=2 MechPolynomialDegree=2 CalculateReference=0
+# - salloc -p cpuonly -t 72:00:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=3 MechPolynomialDegree=2 CalculateReference=0
+ - salloc -p cpuonly -t 72:00:00 -n 1024 -N 40 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=4 MechPolynomialDegree=2 CalculateReference=0 PressureSteps=10
+# - salloc -p cpuonly -t 72:00:00 -n 256 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=5 MechPolynomialDegree=2 CalculateReference=0
+
+#TODO: Zeiten anpassen
+#landbenchmark-deg2_ref:
+# extends: .land-benchmark_on_horeka
+# script:
+# - salloc -p cpuonly -t 72:00:00 -n 1024 -N 20 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=2 CalculateReference=1 ReferenceLevel=4
+# - salloc -p cpuonly -t 72:00:00 -n 512 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=2 CalculateReference=1 ReferenceLevel=4
+# - salloc -p cpuonly -t 72:00:00 -n 1024 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=2 CalculateReference=1 ReferenceLevel=4
+# - salloc -p cpuonly -t 72:00:00 -n 1024 mpirun M++ elasticity/benchmarks/LandBenchmark MechLevel=0 MechPolynomialDegree=2 CalculateReference=1 ReferenceLevel=5
+
+
+EG_landbenchmark-deg1_ref:
+ extends: .land-benchmark_on_horeka
+ script:
+ # - salloc -p cpuonly -t 72:00:00 -n 1024 -N 30 mpirun M++ elasticity/benchmarks/LandBenchmark Mesh=OrientedEllipsoid_3 MechLevel=0 MechPolynomialDegree=1 CalculateReference=1 ReferenceLevel=5
+ - salloc -p cpuonly -t 72:00:00 -n 1024 -N 30 mpirun M++ elasticity/benchmarks/LandBenchmark Mesh=OrientedEllipsoid_3 MechLevel=3 MechPolynomialDegree=1 CalculateReference=0
diff --git a/src/CardMechIProblem.hpp b/src/CardMechIProblem.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1a53c04f484b9478c0291b98f86703201e63297f
--- /dev/null
+++ b/src/CardMechIProblem.hpp
@@ -0,0 +1,64 @@
+#ifndef CARDMECH_CARDMECHIPROBLEM_H
+#define CARDMECH_CARDMECHIPROBLEM_H
+
+#include "Config.hpp"
+#include "Meshes.hpp"
+#include "MeshesCreator.hpp"
+#include <memory>
+
+class CardMechIProblem {
+protected:
+ int verbose = 1;
+
+ std::string meshesName = "";
+
+ std::shared_ptr<Meshes> meshes = nullptr;
+
+ std::shared_ptr<CoarseGeometry> cGeo = nullptr;
+
+public:
+ explicit CardMechIProblem(std::string meshName = "") : meshesName(meshName) {
+ if (meshesName.empty())
+ Config::Get("Mesh", meshesName);
+ if (!meshesName.empty()) {
+ CreateMeshes(MeshesCreator(meshesName));
+ }
+ Config::Get("ProblemVerbose", verbose);
+ }
+
+ explicit CardMechIProblem(std::shared_ptr<Meshes> meshes)
+ : meshes(std::move(meshes)) {
+ Config::Get("ProblemVerbose", verbose);
+ }
+
+ explicit CardMechIProblem(CardMechIProblem &otherProblem)
+ : meshes(otherProblem.meshes), verbose(otherProblem.verbose) {}
+
+ virtual ~CardMechIProblem() = default;
+
+ virtual bool HasExactSolution() const { return false; }
+
+ void CreateMeshes(MeshesCreator meshesCreator) {
+ std::string injectedMeshName = meshesCreator.GetMeshName();
+ cGeo = CreateCoarseGeometryShared(
+ injectedMeshName.empty() ? meshesName : injectedMeshName);
+ meshes = meshesCreator.WithCoarseGeometry(cGeo).CreateShared();
+ }
+
+ const CoarseGeometry &Domain() const { return *cGeo; };
+
+ const Meshes &GetMeshes() const {
+ if (meshes == nullptr)
+ Exit("Meshes not initialized") return *meshes;
+ }
+
+ bool IsMeshInitialized() const { return meshes != nullptr; }
+
+ const Mesh &GetMesh(int level) const { return (*meshes)[level]; }
+
+ const Mesh &GetMesh(LevelPair level) const { return (*meshes)[level]; }
+
+ virtual std::string Name() const = 0;
+};
+
+#endif // CARDMECH_CARDMECHIPROBLEM_H
diff --git a/src/CardiacData.cpp b/src/CardiacData.cpp
index 6fe36143278bace00cf176bf0f210fb870d4b432..7bfc2a845a9fe6bb7a979f5b620d34e1b3e6e842 100644
--- a/src/CardiacData.cpp
+++ b/src/CardiacData.cpp
@@ -99,7 +99,7 @@ std::vector<int> getIncludedChambers(const string &chamberName) {
#include "LagrangeDiscretization.hpp"
void PlotFibreOrientation(const Vector &U, VtuPlot& plot) {
- LagrangeDiscretization fibreDisc(U.GetDisc().GetMeshes(), 0, SpaceDimension);
+ auto fibreDisc = std::make_shared<LagrangeDiscretization>(U.GetDisc().GetMeshes(), 0, SpaceDimension);
Vector f(0.0, fibreDisc);
Vector s(0.0, fibreDisc);
Vector n(0.0, fibreDisc);
diff --git a/src/CardiacInterpolation.cpp b/src/CardiacInterpolation.cpp
index 1dc9d17ed4cd675fa5c9467981d208a43728958f..c1f89bed82a00173c0f470ae6537b55912d529a4 100644
--- a/src/CardiacInterpolation.cpp
+++ b/src/CardiacInterpolation.cpp
@@ -1,6 +1,5 @@
#include <DGDiscretization.hpp>
#include "LagrangeDiscretization.hpp"
-#include "EGDiscretization.hpp"
#include "MixedEGDiscretization.hpp"
#include "MultiPartDiscretization.hpp"
#include "CardiacInterpolation.hpp"
@@ -411,7 +410,7 @@ void smoothMeshData(const Mesh &mesh, Vector &vec) {
void interpolateLagrangeCellData(const Meshes &meshes, Vector &vec) {
- const auto &coarse = meshes.coarseMesh();
+ const auto &coarse = meshes.coarse();
for (cell c = coarse.cells(); c != coarse.cells_end(); ++c) {
// Fetch coarse cell corners and data
@@ -446,13 +445,13 @@ void interpolateLagrangeData(Vector &vec, int degree) {
std::string interpolation{};
Config::Get("ExternalCurrentSmoothingInSpace", interpolation);
if (interpolation != "discrete" && meshes[0].Name() != "BiVentricle_smoothed" && meshes[0].Name() != "KovachevaBiventricle"&& meshes[0].Name() != "BiventricleDirichlet" && meshes[0].Name()!="KovacheveBiventricleSimpleExcitation_smooth") {
- smoothMeshData(meshes.coarseMesh(), vec);
+ smoothMeshData(meshes.coarse(), vec);
}
InterpolateCellData(meshes);
if (degree > 0) {
- interpolateLagrangeVertexData(meshes.coarseMesh(), vec, pow(2, meshes.Level()) * degree);
+ interpolateLagrangeVertexData(meshes.coarse(), vec, pow(2, meshes.Level()) * degree);
} else {
interpolateLagrangeCellData(meshes, vec);
}
@@ -464,22 +463,24 @@ void InterpolateCellData(const Meshes &meshes) {
}
}
-Vector InterpolateMeshData(const IDiscretization& disc){
+Vector InterpolateMeshData(std::shared_ptr<const IDiscretization> disc){
Vector interpolated(0.0, disc);
- const auto &mesh = disc.GetMeshes().coarseMesh();
+ const auto &mesh = disc->GetMeshes().coarse();
- if (typeid(disc) == typeid(LagrangeDiscretization)) {
- int d = ((const LagrangeDiscretization &) disc).Degree();
+ if (typeid(*disc) == typeid(LagrangeDiscretization)) {
+ const auto &lagrangeDisc = std::dynamic_pointer_cast<const LagrangeDiscretization>(disc);
+ int d = lagrangeDisc->Degree();
interpolated = -100.0;
interpolateLagrangeData(interpolated, d);
- } else if (typeid(disc) == typeid(MultiPartDiscretization)) {
- int d = ((const MultiPartDiscretization &) disc).Degree();
+ } else if (typeid(*disc) == typeid(MultiPartDiscretization)) {
+ const auto &multiPartDisc = std::dynamic_pointer_cast<const MultiPartDiscretization>(disc);
+ int d = multiPartDisc->Degree();
interpolated = -100.0;
interpolateLagrangeData(interpolated, d);
- } else if(typeid(disc) == typeid(DGDiscretization)) {
+ } else if(typeid(*disc) == typeid(DGDiscretization)) {
interpolated = -100.0;
interpolateLagrangeData(interpolated, 0);
- } else if(typeid(disc) == typeid(MixedEGDiscretization)) {
+ } else if(typeid(*disc) == typeid(MixedEGDiscretization)) {
interpolated = -100.0;
interpolateLagrangeData(interpolated, 0);
} else {
@@ -487,4 +488,4 @@ Vector InterpolateMeshData(const IDiscretization& disc){
}
return interpolated;
-}
\ No newline at end of file
+}
diff --git a/src/CardiacInterpolation.hpp b/src/CardiacInterpolation.hpp
index adf837f7e1c70f66a4e8f5d6998d1c9e0ea7ef0d..3f19fb1d36a3c14f718485be0dec389e1370f85f 100644
--- a/src/CardiacInterpolation.hpp
+++ b/src/CardiacInterpolation.hpp
@@ -11,7 +11,7 @@ void InterpolateCellData(const Meshes &meshes);
* @param disc Discretization which is used by the assemble requesting the data
* @return A vector containing interpolated vertex data
*/
-Vector InterpolateMeshData(const IDiscretization &disc);
+Vector InterpolateMeshData(std::shared_ptr<const IDiscretization> disc);
#endif //CARDIACINTERPOLATION_HPP
diff --git a/src/MainCardMech.cpp b/src/MainCardMech.cpp
index abb0d20f039a7c61cf29d1e4ab81ce7ef2789a2c..c37015aaa1f4fbf65a58f621619fe3d3d4325428 100644
--- a/src/MainCardMech.cpp
+++ b/src/MainCardMech.cpp
@@ -2,6 +2,7 @@
#include "electrophysiology/MainMonodomain.hpp"
#include "elasticity/MultiElasticity.hpp"
+#include "elasticity/MultiDiscElasticity.hpp"
#include "coupled/MainCoupled.hpp"
@@ -10,10 +11,7 @@ int main(int argc, char **argv) {
Config::SetConfPath(std::string(ProjectSourceDir) + "/conf/");
Config::SetConfigFileName("m++.conf");
Mpp::initialize(&argc, argv);
-
-
Config::PrintInfo();
-
bool clear = false;
Config::Get("ClearData", clear);
if (clear && PPM->master()) {
@@ -26,7 +24,9 @@ int main(int argc, char **argv) {
Config::Get("Model", model);
bool reference{false};
+ bool multiDiscs{false};
Config::Get("CalculateReference", reference);
+ Config::Get("multiDiscs", multiDiscs);
MainCardMech *cmMain;
if (model == "Monodomain") {
@@ -36,6 +36,8 @@ int main(int argc, char **argv) {
mout << "Starting Elasticity Suite of M++." << endl;
if (reference)
cmMain = new MultiElasticity();
+ else if(multiDiscs)
+ cmMain = new MultiDiscElasticity();
else
cmMain = new MainElasticity();
} else if (model.find("Coupled") != std::string::npos) {
diff --git a/src/coupled/MainCoupled.cpp b/src/coupled/MainCoupled.cpp
index a72ea2195d4514e7cbda6d08516f1002e112ccc0..8215a43c44eea7e85dc13e811ce45ab70d7cd243 100644
--- a/src/coupled/MainCoupled.cpp
+++ b/src/coupled/MainCoupled.cpp
@@ -65,9 +65,9 @@ Vector &MainCoupled::Run() {
}
void MainCoupled::generateStartVectors() {
- potential = std::make_unique<Vector>(elphyA->GetDisc());
- stretch = std::make_unique<Vector>(elphyA->GetDisc());
- displacement = std::make_unique<Vector>(0.0, mechA->GetDisc());
+ potential = std::make_unique<Vector>(elphyA->GetSharedDisc());
+ stretch = std::make_unique<Vector>(elphyA->GetSharedDisc());
+ displacement = std::make_unique<Vector>(0.0, mechA->GetSharedDisc());
// Gets initial value of cellmodel at src/cellmodesl/*.cpp
*potential = tensionModel->InitialValue(ELPHY);
diff --git a/src/coupled/problem/CoupledProblem.cpp b/src/coupled/problem/CoupledProblem.cpp
index e75f9677371465062e4052466ac0b6c6ac69fc17..9657a51485abc3eb68c52d34161f37fba724e435 100644
--- a/src/coupled/problem/CoupledProblem.cpp
+++ b/src/coupled/problem/CoupledProblem.cpp
@@ -66,16 +66,14 @@ const std::string &CoupledProblem::MechMeshName() const {
void CoupledProblem::CreateElphyMesh(MeshesCreator &creator) {
{
elphyGeo = CreateCoarseGeometryShared(elphyMeshName);
- elphyMeshes = creator.WithCoarseGeometry(elphyGeo).
- WithCommSplit(ElphyProblem::commSplit).CreateShared();
+ elphyMeshes = creator.WithCoarseGeometry(elphyGeo).CreateShared();
}
}
void CoupledProblem::CreateMechMesh(MeshesCreator &creator) {
{
mechGeo = CreateCoarseGeometryShared(mechMeshName);
- mechMeshes = creator.WithCoarseGeometry(mechGeo).
- WithCommSplit(ElasticityProblem::commSplit).CreateShared();
+ mechMeshes = creator.WithCoarseGeometry(mechGeo).CreateShared();
InterpolateCellData(*mechMeshes);
}
}
diff --git a/src/coupled/problem/DeformedProblem.hpp b/src/coupled/problem/DeformedProblem.hpp
index edcaa009a10b7c26230e80fcb6c4691597507901..f08c98e0090c21e7f1a920ec121ed9df462b2bd7 100644
--- a/src/coupled/problem/DeformedProblem.hpp
+++ b/src/coupled/problem/DeformedProblem.hpp
@@ -8,7 +8,7 @@
class DeformedCoupledProblem : public CoupledProblem {
std::unique_ptr<CoupledProblem> originalProblem{nullptr};
std::unique_ptr<MultiPartDeformationTransfer> transfer;
- std::unique_ptr<IDiscretization> defDisc;
+ std::shared_ptr<IDiscretization> defDisc;
Vector* displacement; //TODO: has to be more efficient
Vector* displacement_old;
std::unique_ptr<Vector> velocity;
@@ -156,11 +156,11 @@ public:
if (!transfer) {
if (contains(v.GetDisc().DiscName(), "MultiPart")) {
const auto &elphyDisc = (const MultiPartDiscretization &) v.GetDisc();
- defDisc = std::make_unique<MultiPartDiscretization>(elphyDisc, u.GetMesh().dim());
+ defDisc = std::make_shared<MultiPartDiscretization>(elphyDisc, u.GetMesh().dim());
}
- displacement = new Vector(0.0, *defDisc);
- displacement_old = new Vector(0.0, *defDisc);
- velocity = std::make_unique<Vector>(0.0, *defDisc);
+ displacement = new Vector(0.0, defDisc);
+ displacement_old = new Vector(0.0, defDisc);
+ velocity = std::make_unique<Vector>(0.0, defDisc);
transfer = GetCardiacTransfer(u, *displacement);
transfer->Prolongate(u, *displacement);
}
diff --git a/src/coupled/solvers/SegregatedSolver.cpp b/src/coupled/solvers/SegregatedSolver.cpp
index 06053ab6dbb8bfd5e0dd091139b026a337d189ba..ea9fd8f3bd286f302714406220c0c02b8873d193 100644
--- a/src/coupled/solvers/SegregatedSolver.cpp
+++ b/src/coupled/solvers/SegregatedSolver.cpp
@@ -186,8 +186,6 @@ bool SegregatedSolver::Method(IElphyAssemble &elphyAssemble, IElasticity &mechAs
return true;
}
-
-
bool SegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &mechAssemble,
Vectors &elphyValues,
Vectors &mechValues) {
@@ -319,7 +317,7 @@ bool adaptiveSegregatedSolver::Step(IElphyAssemble &elphyAssemble, IElasticity &
// Projects stretch from elphy to mech
cardiacTransfer->Project(u_stretch(mechValues), v_stretch(elphyValues));
mechAssemble.UpdateStretch(u_stretch(mechValues));
-
+
vout(2) << "Solving MechStep from " << mechAssemble.Time() << " to "
<< mechAssemble.NextTimeStep(false) << endl;
diff --git a/src/coupled/transfers/MultiPartTransfer.cpp b/src/coupled/transfers/MultiPartTransfer.cpp
index 3b3b680fba390b28b95d07c20b381bf34e9b470b..3f5345816fb60f04436969707299dfc1645be79b 100644
--- a/src/coupled/transfers/MultiPartTransfer.cpp
+++ b/src/coupled/transfers/MultiPartTransfer.cpp
@@ -81,9 +81,9 @@ MultiPartLagrangeTransfer::MultiPartLagrangeTransfer(const Vector &coarse, const
mechDim = mechDisc.Size();
elphyDim = elphyDisc.Size();
- discs = std::pair<std::unique_ptr<LagrangeDiscretization>, std::unique_ptr<MultiPartDiscretization>>{
- std::make_unique<LagrangeDiscretization>(mechDisc, mechDisc.Degree(), mechDisc.Size()),
- std::make_unique<MultiPartDiscretization>(elphyDisc, mechDisc.Size()),
+ discs = std::pair<std::shared_ptr<LagrangeDiscretization>, std::shared_ptr<MultiPartDiscretization>>{
+ std::make_shared<LagrangeDiscretization>(mechDisc, mechDisc.Degree(), mechDisc.Size()),
+ std::make_shared<MultiPartDiscretization>(elphyDisc, mechDisc.Size()),
};
constructRowList(coarse, fine);
@@ -252,7 +252,7 @@ RMatrix MultiPartLagrangeTransfer::AsMatrix() const {
}
std::unique_ptr<Vector> MultiPartLagrangeTransfer::CreateFineDeformation() const {
- return std::make_unique<Vector>(0.0, *discs.second);
+ return std::make_unique<Vector>(0.0, discs.second);
}
std::unique_ptr<MultiPartDeformationTransfer>
diff --git a/src/coupled/transfers/MultiPartTransfer.hpp b/src/coupled/transfers/MultiPartTransfer.hpp
index 973176d2e028dbdc1f32ad5d66140da16f1a73b6..51b253a3d5fbcdc0e6837b70a652a0dc352dcafb 100644
--- a/src/coupled/transfers/MultiPartTransfer.hpp
+++ b/src/coupled/transfers/MultiPartTransfer.hpp
@@ -47,7 +47,7 @@ class MultiPartLagrangeTransfer : public MultiPartDeformationTransfer {
int mechDim{1};
int elphyDim{1};
- std::pair<std::unique_ptr<LagrangeDiscretization>, std::unique_ptr<MultiPartDiscretization>> discs{};
+ std::pair<std::shared_ptr<LagrangeDiscretization>, std::shared_ptr<MultiPartDiscretization>> discs{};
RowList rowList;
void constructRowList(const Vector &coarse, const Vector &fine);
diff --git a/src/elasticity/CMakeLists.txt b/src/elasticity/CMakeLists.txt
index 9b14dc9f38f745320cc4ab26e7b2e90445eb979f..c6c299770570e1530d879036c2a61a2f806cea8a 100644
--- a/src/elasticity/CMakeLists.txt
+++ b/src/elasticity/CMakeLists.txt
@@ -6,7 +6,9 @@ add_subdirectory(solvers)
add_library(ELASTICITY
MainElasticity.cpp
- MultiElasticity.cpp)
+ MultiElasticity.cpp
+ MultiDiscElasticity.cpp
+)
set(CM_LOG data_processing/VectorLogger.cpp
diff --git a/src/elasticity/MainElasticity.cpp b/src/elasticity/MainElasticity.cpp
index 249207644ba2a31e4812af9adbf019fb53deb45a..920d99f6836e12eec2d612130749be713bd4d549 100644
--- a/src/elasticity/MainElasticity.cpp
+++ b/src/elasticity/MainElasticity.cpp
@@ -25,13 +25,28 @@ void MainElasticity::Initialize() {
generateMesh(elasticityModel);
initAssemble();
+ generateStartVectors();
+}
+void MainElasticity::Initialize(std::unique_ptr<ElasticityProblem> problem) {
+ elasticityProblem = std::move(problem);
+ generateMesh(elasticityProblem->Model());
+ initAssemble();
generateStartVectors();
}
Vector &MainElasticity::Run(const Vector &start) {
- // load Interpolation from used discretization model from src/elasticity/assemble/IElasticity.hpp
- *displacement = mechA->Interpolation(start, *displacement);
+ if (start.GetDisc().DiscName() == (*displacement).GetDisc().DiscName() && start.Level() < (*displacement).Level()){
+ // load Interpolation from used discretization model from src/elasticity/assemble/IElasticity.hpp
+ *displacement = mechA->Interpolation(start, *displacement);
+ }
+ else {
+ if(start.GetDisc().DiscName()==(*displacement).GetDisc().DiscName())
+ *displacement = start;
+ else {
+ *displacement = mechA->UseLagrangeAsStartVector(start, *displacement);
+ }
+ }
return Run();
}
@@ -55,10 +70,12 @@ Vector &MainElasticity::Run() {
PrintInfoEntry("Poisson's Ratio",
material::poissonInvariant(elasticityProblem->GetMaterials().first),
1),
+ PrintInfoEntry("MechPolynomialDegree:", discDegree, 1),
PrintInfoEntry("Degrees of Freedom:", displacement->pSize(), 1)/*,
PrintInfoEntry("Matrix Entries", solution.pMatrixSize(), 2)*/);
mout << "===========================================================" << endl << endl;
-
+// displacement->GetMatrixGraph().PrintMatrixMemoryInfo();
+// displacement->GetMatrixGraph().PrintVectorMemoryInfo();
displacement->GetMesh().PrintInfo();
mout << endl << "===========================================================" << endl << endl;
@@ -100,11 +117,6 @@ void MainElasticity::generateMesh(const string &model) const {
WithMeshName(elasticityProblem->GetMeshName()).
WithIncludeVector(chambers).
WithLevel(mlevel);
- if (model.find("DG") != std::string::npos) {
- mCreator = mCreator.WithOverlap("dG1");
- } else if (model.find("EG") != std::string::npos) {
- mCreator.WithOverlap("EG");
- } else mCreator.WithoutOverlap();
elasticityProblem->CreateMeshes(mCreator);
InterpolateCellData(elasticityProblem->GetMeshes());
@@ -115,7 +127,7 @@ void MainElasticity::generateStartVectors() {
vectorDisc = mechA->CreateDisc(3);
scalarDisc = mechA->CreateDisc(1);
- displacement = std::make_unique<Vector>(*vectorDisc);
+ displacement = std::make_unique<Vector>(vectorDisc);
*displacement = 0.0;
//displacement->Clear();
}
@@ -138,11 +150,9 @@ Vector &MainElasticity::GetDisplacement() const {
}
void MainElasticity::initAssemble() {
- int degree = 1;
bool isStatic = false;
- Config::Get("MechPolynomialDegree", degree);
// Create assemble for used discretization in src/elasticity/assemble/IElasticity.cpp
- mechA = CreateFiniteElasticityAssemble(*elasticityProblem, degree, isStatic);
+ mechA = CreateFiniteElasticityAssemble(*elasticityProblem, discDegree, isStatic);
auto times = ReadTimeFromConfig();
mechA->ResetTime(times[0], times[1], times[2]);
}
@@ -221,6 +231,11 @@ void MainElasticity::ResetLevel(int newLevel) {
mlevel = newLevel;
}
+
+void MainElasticity::SetDegree(int degree){
+ discDegree = degree;
+}
+
void MainElasticity::EvaluateReferece(const Vector &solution, const Vector &reference) {
mout << "============== Comparing grid level " << solution.GetMesh().Level() << " with "
<< reference.GetMesh().Level() << " =========" << endl << endl;
diff --git a/src/elasticity/MainElasticity.hpp b/src/elasticity/MainElasticity.hpp
index 889eb56e81261ccdfcb39fa1516de8a9d04b576e..32d8f189bc37e8e87291290cdafd6425208a2f38 100644
--- a/src/elasticity/MainElasticity.hpp
+++ b/src/elasticity/MainElasticity.hpp
@@ -9,7 +9,7 @@
/// Is used to move ElasticityProblem, Discretization and displacement in MoveVector() function.
struct VectorStruct {
std::unique_ptr<ElasticityProblem> problem{nullptr};
- std::unique_ptr<IDiscretization> disc{nullptr};
+ std::shared_ptr<IDiscretization> disc{nullptr};
std::unique_ptr<Vector> vec{nullptr};
Vector &operator()() {
@@ -26,6 +26,7 @@ struct VectorStruct {
*/
class MainElasticity : public MainCardMech {
friend class MultiElasticity;
+ friend class MultiDiscElasticity;
void generateMesh(const string &model) const;
@@ -34,8 +35,8 @@ class MainElasticity : public MainCardMech {
std::unique_ptr<ElasticityProblem> elasticityProblem = nullptr;
- std::unique_ptr<IDiscretization> vectorDisc = nullptr;
- std::unique_ptr<IDiscretization> scalarDisc = nullptr;
+ std::shared_ptr<IDiscretization> vectorDisc = nullptr;
+ std::shared_ptr<IDiscretization> scalarDisc = nullptr;
std::unique_ptr<IElasticity> mechA = nullptr;
std::unique_ptr<IElasticity> prestressedEla = nullptr;
@@ -65,10 +66,14 @@ public:
void Initialize() override;
+ void Initialize(std::unique_ptr<ElasticityProblem> problem);
+
void ResetTime(double start, double end, int steps) override;
void ResetLevel(int newLevel);
+ void SetDegree(int degree);
+
Vector &Run() override;
std::vector<double> Evaluate(const Vector &solution) override;
diff --git a/src/elasticity/MultiDiscElasticity.cpp b/src/elasticity/MultiDiscElasticity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..24603cbfca46754cdf2857fb2a7cd39711651e89
--- /dev/null
+++ b/src/elasticity/MultiDiscElasticity.cpp
@@ -0,0 +1,93 @@
+#include "MultiDiscElasticity.hpp"
+
+
+
+MultiDiscElasticity::MultiDiscElasticity() : main() {
+ Config::Get("multiDiscs", multiDiscs);
+ Config::Get("useDifferentMeshes", useDifferentMeshes);
+ Config::Get("useDifferentDegs", useDifferentDegs);
+ Config::Get("MechDiscretizationList", discs);
+ solutions.resize(discs.size());
+ meshes.resize(discs.size());
+ degrees.resize(discs.size());
+ DGPenalties.resize(discs.size());
+
+ Config::Get("MechPolynomialDegreeList", degreeList);
+ Config::Get("MeshList", meshesList);
+ Config::Get("DGPenaltyList", DGPenaltyList);
+}
+
+
+void MultiDiscElasticity::Initialize() {
+
+ if (degreeList.size() == 1){
+ for(int i = 0; i <= discs.size()-1; ++i) {
+ degrees[i] = degreeList[0];
+ }
+ } else {
+ if (degreeList.size() != degrees.size()) THROW("List of MechPolynomialDegrees is not matching!")
+ degrees = degreeList;
+ }
+
+
+ if (meshesList.size() == 1){
+ for(int i = 0; i <= discs.size()-1; ++i) {
+ meshes[i] = meshesList[0];
+ }
+ } else {
+ if (meshesList.size() != meshes.size()) THROW("List of MechPolynomialDegrees is not matching!")
+ meshes = meshesList;
+ }
+
+ DGPenalties = DGPenaltyList;
+
+ for (int j = 0; j < degrees.size(); ++j){
+ if (degrees[j] == degrees[j+1]){
+ useDifferentDegs = false;
+ }
+ }
+ auto elasticityProblemLocal = GetElasticityProblem();
+ elasticityProblemLocal->SetModel(discs[0]); // sets disc
+ elasticityProblemLocal->SetDGPenalty(DGPenalties[0]);
+ if (useDifferentMeshes) {
+ elasticityProblemLocal->SetMeshName(meshes[0]);
+ }
+ main.SetDegree(degrees[0]);
+ main.Initialize(std::move(elasticityProblemLocal));
+}
+
+Vector &MultiDiscElasticity::Run() {
+ for (int i = 0; i <= discs.size()-1; ++i) {
+ mout << "\n=================== Running discretization " << discs[i] << " ================ " << endl;
+ auto elasticityProblemLocal = GetElasticityProblem();
+ elasticityProblemLocal->SetModel(discs[i]);
+
+ elasticityProblemLocal->SetDGPenalty(DGPenalties[i]);
+
+ if (useDifferentMeshes) {
+ elasticityProblemLocal->SetMeshName(meshes[i]);
+ }
+ main.SetDegree(degrees[i]);
+ main.Initialize(std::move(elasticityProblemLocal));
+
+ bool usePrevDisc = (multiDiscs && i > 0 && !useDifferentMeshes && !useDifferentDegs);
+
+ if ((usePrevDisc && useDifferentMeshes) || (usePrevDisc && useDifferentDegs)){
+ THROW("Not implemented right now!")
+ }
+
+ main.pressureSteps = 1;
+ const auto &solution = usePrevDisc ? main.Run(solutions[i - 1]()) : main.Run();
+
+ main.Evaluate(solution);
+ solutions[i] = main.MoveVector();
+ }
+
+ return solutions[discs.size() - 1]();
+}
+
+vector<double> MultiDiscElasticity::Evaluate(const Vector &solution) {
+ mout << endl << "=============================================================" << endl << endl;
+ return {};
+}
+
diff --git a/src/elasticity/MultiDiscElasticity.hpp b/src/elasticity/MultiDiscElasticity.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b12cc4ecfcb2955a125f92e667e942c1e032f70c
--- /dev/null
+++ b/src/elasticity/MultiDiscElasticity.hpp
@@ -0,0 +1,46 @@
+//
+// Created by laura on 30.01.24.
+//
+
+#ifndef CARDMECH_MULTIDISCELASTICITY_HPP
+#define CARDMECH_MULTIDISCELASTICITY_HPP
+
+#endif //CARDMECH_MULTIDISCELASTICITY_HPP
+
+#include "MainCardMech.hpp"
+#include "MainElasticity.hpp"
+
+/**
+ * Running different discretizations one after the other.
+ * Homotopy (for different discs or e.g. to test different penalty parameters for DG/EG faster).
+ * Possibilty to choose different meshes and degrees is not implemented yet!
+ * "meshesList" and "degreeList": if there is only one entry in the config, it uses the same entries
+ * for all discretizations in "discs"
+ * homotopy possibilities: CG -> DG, CG -> EG, DG -> DG, EG -> EG
+ */
+class MultiDiscElasticity : public MainCardMech {
+ MainElasticity main;
+
+ bool multiDiscs{false};
+ bool useDifferentMeshes{false};
+ bool useDifferentDegs{false};
+ std::vector<VectorStruct> solutions{};
+ std::vector<string> discs{};
+ std::vector<string> meshes;
+ std::vector<string> meshesList;
+ std::vector<int> degrees;
+ std::vector<int> degreeList;
+ std::vector<int> DGPenalties;
+ std::vector<int> DGPenaltyList;
+
+public:
+ MultiDiscElasticity();
+
+ void Initialize() override;
+
+ Vector &Run() override;
+
+ vector<double> Evaluate(const Vector &solution) override;
+
+};
+
diff --git a/src/elasticity/assemble/DGElasticity.cpp b/src/elasticity/assemble/DGElasticity.cpp
index 5261b199ddeccdcaf764c9ba0bef7278ebebdd22..ab9428a22cb218db3ce73129999623e2629a7cd7 100644
--- a/src/elasticity/assemble/DGElasticity.cpp
+++ b/src/elasticity/assemble/DGElasticity.cpp
@@ -2,11 +2,15 @@
#include "DGDoF.hpp"
#include "DGVectorFieldElement.hpp"
#include "VectorAccess.hpp"
+#include "VectorFieldElement.hpp"
DGElasticity::DGElasticity(ElasticityProblem &eP, int degree, bool isStatic)
- : IElasticity(eP), disc(eP.GetMeshes(), degree, dim()) {
+ : IElasticity(eP), disc(std::make_shared<const DGDiscretization>(eP.GetMeshes(), degree, dim())) {
Config::Get("DGSign", sign);
- Config::Get("DGPenalty", penalty);
+
+ penalty = eP.GetDGPenalty();
+ if (penalty == 0) Config::Get("DGPenalty", penalty);
+
this->degree = degree;
}
@@ -193,9 +197,9 @@ void DGElasticity::Residual(const cell &c, const Vector &U, Vector &r) const {
}
} else if (E.Bnd(f) == -1) {
- cell cf = r.GetMesh().find_neighbour_cell(c, f);
+ cell cf = r.find_neighbour_cell(c, f);
if (cf() < c()) continue;
- int f1 = r.GetMesh().find_neighbour_face_id(c.Face(f), cf);
+ int f1 = r.find_neighbour_face_id(c.Face(f), cf);
DGVectorFieldElement E_nghbr(U, *cf);
DGVectorFieldFaceElement FE_nghbr(r, *cf, f1);
DGSizedRowBndValues r_cf(r, *cf, E_nghbr.shape_size());
@@ -481,14 +485,14 @@ void DGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
}
}
if (E.Bnd(face) == -1) {
- cell cf = U.GetMesh().find_neighbour_cell(c, face);
+ cell cf = U.find_neighbour_cell(c, face);
if (cf() < c()) continue;
DGVectorFieldElement E_nghbr(U, *cf);
DGSizedRowEntries A_cf(A, *c, E.shape_size(), *cf, E_nghbr.shape_size());
DGSizedRowEntries A_fc(A, *cf, E_nghbr.shape_size(), *c, E.shape_size());
DGSizedRowEntries A_ff(A, *cf, E_nghbr.shape_size(), *cf, E_nghbr.shape_size());
- int f1 = U.GetMesh().find_neighbour_face_id(c.Face(face), cf);
+ int f1 = U.find_neighbour_face_id(c.Face(face), cf);
DGVectorFieldFaceElement FE_nghbr(U, *cf, f1);
for (int q = 0; q < FE.nQ(); ++q) {
@@ -654,6 +658,24 @@ void DGElasticity::Interpolate(const Vector &coarse, Vector &u) const {
u.Accumulate();
}
+void DGElasticity::UseLagrangeValuesAsStartVector(const Vector &u_conforming, Vector &u) const {
+ for (cell c = u.cells(); c != u.cells_end(); ++c) {
+ DGVectorFieldElement E(u, *c);
+ VectorFieldElement EC(u_conforming, *c);
+ for (int i = 0; i < E.shape_size(); ++i) {
+ Point x = E.NodalPoint(i);
+ row r = u_conforming.find_row(x);
+ if (r == u_conforming.rows_end()) continue;
+ for (int k = 0; k < E.Dim(); ++k)
+ u(c(), k * E.shape_size() + i) = u_conforming(r,k);
+ }
+ auto x = c.GlobalToLocal(c());
+ VectorField val = E.VectorValue(x, u);
+ VectorField valC = EC.VectorValue(x, u_conforming);
+ }
+ u.MakeAdditive();
+ u.Accumulate();
+}
double DGElasticity::StrainEnergyError(const Vector &u) const {
double err = 0.0;
@@ -1070,7 +1092,7 @@ double DGElasticity::DeformedVolume(const Vector &U) const {
void DGElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {
if (vtkplot < 1) return;
- Vector bnd(0.0, U.GetDisc());
+ Vector bnd(0.0, U.GetSharedDisc());
auto &plot = mpp::plot("Boundary");
for (cell c = U.cells(); c != U.cells_end(); c++) {
@@ -1081,17 +1103,36 @@ void DGElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {
for (int face = 0; face < c.Faces(); ++face) {
DGVectorFieldFaceElement FE(U, *c, face);
int bc = u_c.bc(face);
-
- for (int i = 0; i < FE.shape_size(); ++i) {
- bnd(c(), i) = max(bnd(c(), i), (double) bc);
- }
-
+ bnd(c(), 0) = max(bnd(c(), 0), (double) bc);
+ bnd(c(), 1) = max(bnd(c(), 1), bc == 199 ? 1.0 : 0.0); // Fixation Boundary bc == 199
}
}
- plot.AddData("Boundary", bnd);
+ plot.AddData("Boundary", bnd, 0);
+ plot.AddData("Fixation", bnd, 1);
plot.PlotFile("Boundary");
}
+Vector makePlottableDGVector(const Vector &U) {
+ Vector plotU(0.0, U);
+ for (cell c = U.cells(); c != U.cells_end(); c++) {
+ row r = plotU.find_row(c());
+ DGVectorFieldElement elem(U, *c);
+ VectorField disp = elem.VectorValue(c.LocalCenter(), U);
+ plotU(r, 0) = disp[0];
+ plotU(r, 1) = disp[1];
+ plotU(r, 2) = disp[2];
+ }
+ return plotU;
+}
+
+
+void DGElasticity::PlotDisplacement(const Vector &U, int step,
+ const string &varname) const {
+ // Plotting the deformation in src/CardiacAssemble.hpp
+ PlotDeformed(U, {{"Displacement", makePlottableDGVector(U)}}, step, varname);
+
+}
+
void DGElasticity::SetDisplacement(Vector &U) const {
for (cell c = U.cells(); c != U.cells_end(); ++c) {
auto index = U.find_row(c());
@@ -1206,8 +1247,6 @@ void DGElasticity::ProblemEvaluation(const Vector &u) const {
}
-
-
void DGElasticity::MassMatrix(Matrix &massMatrix) const {
THROW("MassMatrix is not implemented for DG")
}
@@ -1218,4 +1257,24 @@ void DGElasticity::SystemMatrix(Matrix &systemMatrix) const {
void DGElasticity::Initialize(Vectors &vecs) const {
THROW("Initialize is not implemented for DG")
-}
\ No newline at end of file
+}
+
+
+std::pair<double, double> DGElasticity::detF(const Vector &u) const {
+ double Jmax = 1.0;
+ double Jmin = 1.0;
+
+ for (cell c = u.cells(); c != u.cells_end(); ++c) {
+ DGVectorFieldElement E(u, *c);
+ for (int q = 0; q < E.nQ(); ++q) {
+ Tensor DU = E.VectorGradient(q, u);
+ double J = det(One + DU);
+ Jmin = min(Jmin, J);
+ Jmax = max(Jmax, J);
+ }
+ }
+ Jmin = PPM->Min(Jmin);
+ Jmax = PPM->Max(Jmax);
+
+ return {Jmin, Jmax};
+}
diff --git a/src/elasticity/assemble/DGElasticity.hpp b/src/elasticity/assemble/DGElasticity.hpp
index 751565d87754f5fdbb7b72c1a78b364b38dec973..010a62e350c6efea650988639a5979707aef3445 100644
--- a/src/elasticity/assemble/DGElasticity.hpp
+++ b/src/elasticity/assemble/DGElasticity.hpp
@@ -10,18 +10,20 @@ class DGElasticity : public IElasticity {
double penalty{1};
int size{3};
- const DGDiscretization disc;
+ std::shared_ptr<const DGDiscretization> disc;
protected:
void prestress(const Vector &u, Vector &pStress) override;
public:
DGElasticity(ElasticityProblem &eP, int degree, bool isStatic = false);
- std::unique_ptr<IDiscretization> CreateDisc(int sizes) const override {
- return std::make_unique<DGDiscretization>(disc.GetMeshes(), degree, sizes);
+ std::shared_ptr<IDiscretization> CreateDisc(int sizes) const override {
+ return std::make_shared<DGDiscretization>(disc->GetMeshes(), degree, sizes);
}
- const IDiscretization &GetDisc() const override { return disc; }
+ const IDiscretization &GetDisc() const override { return *disc; }
+
+ std::shared_ptr<const IDiscretization> GetSharedDisc() const override {return disc; }
void Initialize(Vector &u) const override {};
@@ -52,7 +54,9 @@ public:
void Interpolate(const Vector &coarse, Vector &u) const override;
- //void Project(Vector &coarse, const Vector &fine) const;
+ void UseLagrangeValuesAsStartVector(const Vector &u_conforming, Vector &u) const override;
+
+ //void Project(Vector &coarse, const Vector &fine) const;
double L2Error(const Vector &u) const override;
@@ -86,6 +90,8 @@ public:
void PlotBoundary(const Vector &U, const Vector &scal) const override;
+ void PlotDisplacement(const Vector &u, int step, const string &varname) const override;
+
//std::array<double, 4> ChamberVolume(const Vector &U) const override;
void FixationBoundary(const cell &c, int face, int bc, const Vector &U, Vector &r) const;
@@ -118,8 +124,10 @@ public:
void Scale (double s) const {
Material &cellMat = eProblem.GetMaterial();
- cellMat.ScalePenalty(s*s);
+// cellMat.ScalePenalty(s*s);
}
+
+ std::pair<double, double> detF(const Vector &u) const override;
};
#endif //DGELASTICITY_HPP
diff --git a/src/elasticity/assemble/EGElasticity.cpp b/src/elasticity/assemble/EGElasticity.cpp
index d50a0c6c789251dab64eab59889c3b539c496a2f..7b835427f0d6bd76535268443522fd2de1f12ac1 100644
--- a/src/elasticity/assemble/EGElasticity.cpp
+++ b/src/elasticity/assemble/EGElasticity.cpp
@@ -1,9 +1,12 @@
#include "EGElasticity.hpp"
EGElasticity::EGElasticity(ElasticityProblem &eP, int degree, bool isStatic)
- : IElasticity(eP), disc(eP.GetMeshes(), degree, dim()){
+ : IElasticity(eP), disc(std::make_shared<const MixedEGDiscretization>(eP.GetMeshes(), degree, dim())){
Config::Get("DGSign", sign);
- Config::Get("DGPenalty", penalty);
+
+ penalty = eP.GetDGPenalty();
+ if (penalty == 0) Config::Get("DGPenalty", penalty);
+
this->degree = degree;
}
@@ -26,13 +29,23 @@ void EGElasticity::Energy(const cell &c, const Vector &U, double &energy) const
double w = E.QWeight(q);
Tensor F = DeformationGradient(E, q, U);
-
+ Tensor FT = transpose(F);
+ Tensor E = 0.5 * (FT * F - One);
// VectorField v = timeScheme->Velocity<MixedEGVectorFieldElement>(U, E, q); // noch überlegen
// energy += w * rho * (v * v);
- energy += w * cellMat.TotalEnergy(F, Q);
+ // energy += w * cellMat.TotalEnergy(F, Q);
+
+
+ //energy = min(energy,det(F));
+
+ energy = max(energy,cellMat.q(E, Q[0]));
+
}
+ return;
+
+
BFParts bnd(U.GetMesh(), *c);
if (!bnd.onBnd()) return;
@@ -197,9 +210,9 @@ void EGElasticity::Residual(const cell &c, const Vector &U, Vector &r) const {
}
} else if (E.Bnd(f) == -1) {
- cell cf = r.GetMesh().find_neighbour_cell(c, f);
+ cell cf = r.find_neighbour_cell(c, f);
if (cf() < c()) continue;
- int f1 = r.GetMesh().find_neighbour_face_id(c.Face(f), cf);
+ int f1 = r.find_neighbour_face_id(c.Face(f), cf);
MixedEGVectorFieldElement E_nghbr(U, *cf);
MixedEGVectorFieldFaceElement FE_nghbr(r, *cf, f1);
MixedRowValues r_cf(r, *cf);
@@ -444,9 +457,9 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
w * (localPressure) * (Cross(F, Dtheta) * N) * theta;
}
} else if (E.Bnd(face) == -1) {
- cell cf = U.GetMesh().find_neighbour_cell(c, face);
+ cell cf = U.find_neighbour_cell(c, face);
if (cf() < c()) continue;
- int f1 = U.GetMesh().find_neighbour_face_id(c.Face(face), cf);
+ int f1 = U.find_neighbour_face_id(c.Face(face), cf);
MixedEGVectorFieldElement E_nghbr(U, *cf);
MixedEGVectorFieldFaceElement FE_nghbr(U, *cf, f1);
MixedRowEntries A_cf(A, *c, *cf);
@@ -927,30 +940,54 @@ void EGElasticity::GetInvariant(const Vector &displacement, Vector &iota4) const
iota4.Collect();
}
-/*void EGElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {
- if (vtkplot < 1) return;
- Vector bnd(0.0, U.GetDisc());
-
- auto &plot = mpp::plot("Boundary");
- for (cell c = U.cells(); c != U.cells_end(); c++) {
- MixedEGVectorFieldElement E(U, *c);
- RowBndValues u_c(bnd, *c);
- if (!u_c.onBnd()) continue;
-
- for (int face = 0; face < c.Faces(); ++face) {
- MixedEGVectorFieldFaceElement FE(U, *c, face);
- int shape_size = FE.dimension()-1;
- int bc = u_c.bc(face);
+void EGElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {
+ if (vtkplot < 1) return;
+ Vector bnd(scal);
+ bnd.Clear();
+ Vector count(bnd);
+ count.Clear();
+ auto &plot = mpp::plot("Boundary");
+ for (cell c = U.cells(); c != U.cells_end(); c++) {
+ RowBndValues u_c(bnd, *c);
+ RowBndValues u_count(count, *c);
+
+ if (!u_c.onBnd()) continue;
+
+ for (int i = 0; i < c.Faces(); ++i) {
+ for (int j = 0; j < U.NumberOfNodalPointsOnFace(*c, i); ++j) {
+ int id = U.IdOfNodalPointOnFace(*c, i, j);
+ if (u_c.bc(i) > u_c(id, 0)) {
+ u_c(id, 0) = (double) u_c.bc(i);
+ u_count(id, 0) = 1;
+ }
+ }
+ }
+ }
+ bnd.Accumulate();
+ count.Accumulate();
+ for (size_t i = 0; i < count.size(); i++) {
+ if (count[i] > 0) {
+ bnd[i] /= count[i];
+ }
+ }
+ plot.AddData("Boundary", bnd, 0);
- for (int i = 0; i < shape_size; ++i) {
- bnd(c(), i) = max(bnd(c(), i), (double) bc);
- }
+ bnd = 0;
+ for (cell c = U.cells(); c != U.cells_end(); c++) {
+ RowBndValues u_c(bnd, *c);
+ if (!u_c.onBnd()) continue;
- }
+ for (int i = 0; i < c.Faces(); ++i) {
+ for (int j = 0; j < U.NumberOfNodalPointsOnFace(*c, i); ++j) {
+ int id = U.IdOfNodalPointOnFace(*c, i, j);
+ u_c(id, 0) = max(u_c(id, 0), u_c.bc(i) == 199 ? 1.0 : 0.0);
+ }
}
- plot.AddData("Boundary", bnd);
- plot.PlotFile("Boundary");
-}*/
+ }
+ plot.AddData("Fixation", bnd);
+
+ plot.PlotFile("Boundary");
+}
void EGElasticity::SetDisplacement(Vector &U) const {
for (row r = U.rows(); r != U.rows_end(); ++r) {
@@ -1036,6 +1073,22 @@ void EGElasticity::Interpolate(const Vector &coarse, Vector &u) const {
}
+void EGElasticity::UseLagrangeValuesAsStartVector(const Vector &u_conforming, Vector &u) const{
+ for (row r = u.rows(); r != u.rows_end(); ++r){
+ row r2 = u_conforming.find_row(r());
+ if (r2 == u_conforming.rows_end()) continue;
+ int n = r.NumberOfDofs();
+ if (n == 1) continue;
+ for (int j = 0; j < n ; ++j){
+ u(r,j) = u_conforming(r2, j);
+ }
+ }
+ u.MakeAdditive();
+ u.Accumulate();
+}
+
+
+
void EGElasticity::MassMatrix(Matrix &massMatrix) const {
THROW("MassMatrix is not implemented for DG")
}
@@ -1170,4 +1223,24 @@ void EGElasticity::Initialize(const cell &c, Vector &U) const {
break;
}
}
-}
\ No newline at end of file
+}
+
+std::pair<double, double> EGElasticity::detF(const Vector &u) const {
+ double Jmax = 1.0;
+ double Jmin = 1.0;
+
+ for (cell c = u.cells(); c != u.cells_end(); ++c) {
+ MixedEGVectorFieldElement E(u, *c);
+ for (int q = 0; q < E.nQ(); ++q) {
+ Tensor DU = E.VectorGradient(q, u);
+ double J = det(One + DU);
+ Jmin = min(Jmin, J);
+ Jmax = max(Jmax, J);
+ }
+ }
+ Jmin = PPM->Min(Jmin);
+ Jmax = PPM->Max(Jmax);
+
+ return {Jmin, Jmax};
+}
+
diff --git a/src/elasticity/assemble/EGElasticity.hpp b/src/elasticity/assemble/EGElasticity.hpp
index e3c5fb6b56ef255b67c53554e8d839b09042f026..24f70cd21d9c7be64f6518fc1baf29feae70a826 100644
--- a/src/elasticity/assemble/EGElasticity.hpp
+++ b/src/elasticity/assemble/EGElasticity.hpp
@@ -14,7 +14,7 @@ class EGElasticity : public IElasticity {
double penalty{1};
int size{3};
- const MixedEGDiscretization disc;
+ std::shared_ptr<const MixedEGDiscretization> disc;
void fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const;
@@ -26,11 +26,13 @@ protected:
public:
EGElasticity(ElasticityProblem &eP, int degree, bool isStatic = false);
- std::unique_ptr<IDiscretization> CreateDisc(int sizes) const override {
- return std::make_unique<MixedEGDiscretization>(disc.GetMeshes(), degree, sizes);
+ std::shared_ptr<IDiscretization> CreateDisc(int sizes) const override {
+ return std::make_shared<MixedEGDiscretization>(disc->GetMeshes(), degree, sizes);
}
- const IDiscretization &GetDisc() const override { return disc; }
+ const IDiscretization &GetDisc() const override { return *disc; }
+
+ std::shared_ptr<const IDiscretization> GetSharedDisc() const override {return disc; }
//void Initialize(Vector &u) const override {};
@@ -92,7 +94,7 @@ public:
double NormPrestress();
- /*void PlotBoundary(const Vector &U, const Vector &scal) const override;*/
+ void PlotBoundary(const Vector &U, const Vector &scal) const override;
/*std::array<double, 4> CalculateVolume(const Vector &U) const override;*/
@@ -106,8 +108,12 @@ public:
void Interpolate(const Vector &coarse, Vector &u) const override;
+ void UseLagrangeValuesAsStartVector(const Vector &u_conforming, Vector &u) const override;
+
std::array<double, 4> ChamberVolume(const Vector &u) const override;
+ std::pair<double, double> detF(const Vector &u) const override;
+
};
diff --git a/src/elasticity/assemble/IElasticity.cpp b/src/elasticity/assemble/IElasticity.cpp
index 5d7242a134e5ebfd59c5ba5d2d7e402c25521abe..6e4ed0dcd9fc478d09ebe11051de6e28535d0347 100644
--- a/src/elasticity/assemble/IElasticity.cpp
+++ b/src/elasticity/assemble/IElasticity.cpp
@@ -24,24 +24,13 @@ void IElasticity::SetInitialValue(Vector &u) {
// Initialize ElasticityProblem in src/elasticity/problem/ElasticityProblem.hpp
eProblem.Initialize(StepSize(), Time(), ChamberVolume(u));
const Meshes &M = u.GetMeshes();
- disc03 = new LagrangeDiscretization(M, 0, 3);
+ auto disc03 = std::make_shared<LagrangeDiscretization>(M, 0, 3);
gamma_c= !Prestress ? std::make_unique<Vector>(*disc03,prestressLevel) : std::make_unique<Vector>(*disc03);
if (!Prestress) {
Prestress = std::make_unique<Vector>(0.0, u);
}
gamma = std::make_unique<Vector>(0.0, u);
-
-
- /*int cellcounterGamma=0;
- int cellcounterU=0;
- for(cell c = (*gamma_c).cells(); c != (*gamma_c).cells_end(); ++c)
- cellcounterGamma+=1;
- for(cell c = (u).cells(); c != (u).cells_end(); ++c)
- cellcounterU+=1;
- std::cout<<"gamma cells "<<cellcounterGamma<<" u cells "<<cellcounterU<<endl;
-*/
-
}
void IElasticity::SetInitialValue(Vectors &values) {
@@ -182,6 +171,9 @@ void IElasticity::PrintPressureIteration(const Vector &u) const {
}
void IElasticity::PlotIteration(const Vectors &mechValues, int step) const {
+ if (mechValues[0].GetDisc().DiscName().find("DG") != std::string::npos) {
+ THROW("TODO: Override in DGElasticity as in PlotDisplacement!")
+ }
step = step < 0 ? Step() : step;
if (step == 0) PlotBoundary(mechValues[0], mechValues[0]);
diff --git a/src/elasticity/assemble/IElasticity.hpp b/src/elasticity/assemble/IElasticity.hpp
index dcc38cdd194408f2e6245798a7e6daa21f74809f..789f74d4d82ff315a5e2db674964253ed5f50d62 100644
--- a/src/elasticity/assemble/IElasticity.hpp
+++ b/src/elasticity/assemble/IElasticity.hpp
@@ -50,8 +50,10 @@ public:
/// Gets the corresponding discretization.
virtual const IDiscretization &GetDisc() const = 0;
+ virtual std::shared_ptr<const IDiscretization> GetSharedDisc() const = 0;
+
/// Creates a new discretization with a different size.
- virtual std::unique_ptr<IDiscretization> CreateDisc(int size) const = 0;
+ virtual std::shared_ptr<IDiscretization> CreateDisc(int size) const = 0;
const char *Name() const override { return "Cardiac Elasticity"; }
@@ -143,7 +145,7 @@ public:
virtual void GetInvariant(const Vector &u, Vector &iota4, Vector &iota) const { return GetInvariant(u,iota4); }
- void PlotDisplacement(const Vector &u, int step, const string &varname) const;
+ virtual void PlotDisplacement(const Vector &u, int step, const string &varname) const;
/**
* Creates a new vector with the discretization of target (coarse vector) and projects fine
@@ -175,6 +177,17 @@ public:
Warning("No Interpolation implemented")
}
+ Vector UseLagrangeAsStartVector(const Vector &u_conforming, const Vector &target) const {
+ Vector v(target);
+ UseLagrangeValuesAsStartVector(u_conforming, v);
+ return v;
+ }
+
+
+ virtual void UseLagrangeValuesAsStartVector(const Vector &u_conforming, Vector &u) const {
+ Warning("Not implemented for this discretization")
+ }
+
virtual double L2Norm(const Vector &u) const { return 0.0; };
diff --git a/src/elasticity/assemble/LagrangeElasticity.cpp b/src/elasticity/assemble/LagrangeElasticity.cpp
index f8d72b3be2248b72c0cad929e6cfae7367f477d4..9c262f1e462d5fa821972d231c981bfe55a0e252 100644
--- a/src/elasticity/assemble/LagrangeElasticity.cpp
+++ b/src/elasticity/assemble/LagrangeElasticity.cpp
@@ -4,7 +4,7 @@
#include "VectorFieldElement.hpp"
LagrangeElasticity::LagrangeElasticity(ElasticityProblem &eP, int degree, bool isStatic)
- : IElasticity(eP), disc(eP.GetMeshes(), degree, dim()) {
+ : IElasticity(eP), disc(std::make_shared<LagrangeDiscretization>(eP.GetMeshes(), degree, dim())) {
this->degree = degree;
Config::Get("withPressure", withPressure);
@@ -1222,7 +1222,7 @@ void LagrangeElasticity::SetInitialValue(Vector &u) {
void LagrangeElasticity::Initialize(const cell &c, Vector &U) const {
RowBndValues u_c(U, *c);
- if (!U.GetMesh().OnBoundary(*c)) return;
+ if (!U.OnBoundary(*c)) return;
for (int face = 0; face < c.Faces(); ++face) {
switch (u_c.bc(face)) {
case 1:
@@ -1290,7 +1290,7 @@ void LagrangeElasticity::TractionStiffness(Matrix &matrix) const {
for (int f = 0; f < c.Faces(); ++f) {
VectorFieldFaceElement faceElem(matrix, *c, f);
- if (matrix.GetMesh().OnBoundary(*c, f)) {
+ if (matrix.OnBoundary(*c, f)) {
int bc = matrix.GetMesh().BoundaryFacePart(c.Face(f));
if ((bc >= ROBIN_BC_EPI) && (bc <= ROBIN_BC_BASE)) {
@@ -1329,7 +1329,7 @@ void LagrangeElasticity::EvaluateTractionStiffness(Vector &u) const {
for (int f = 0; f < c.Faces(); ++f) {
VectorFieldFaceElement faceElem(u, *c, f);
- if (u.GetMesh().OnBoundary(*c, f)) {
+ if (u.OnBoundary(*c, f)) {
for (int face = 0; face < c.Faces(); ++face) {
int bc = u.GetMesh().BoundaryFacePart(c.Face(f));
@@ -1356,7 +1356,7 @@ void LagrangeElasticity::TractionViscosity(Matrix &matrix) const {
for (int f = 0; f < c.Faces(); ++f) {
VectorFieldFaceElement faceElem(matrix, *c, f);
- if (matrix.GetMesh().OnBoundary(*c, f)) {
+ if (matrix.OnBoundary(*c, f)) {
int bc = matrix.GetMesh().BoundaryFacePart(c.Face(f));
if ((bc >= ROBIN_BC_EPI) && (bc <= ROBIN_BC_BASE)) {
diff --git a/src/elasticity/assemble/LagrangeElasticity.hpp b/src/elasticity/assemble/LagrangeElasticity.hpp
index 84645ad83247d454b7d8e20f8bebd065b8733dd6..71e2589c17250a78ab6be577a402a74c100ff527 100644
--- a/src/elasticity/assemble/LagrangeElasticity.hpp
+++ b/src/elasticity/assemble/LagrangeElasticity.hpp
@@ -6,7 +6,7 @@
class LagrangeElasticity : public IElasticity {
int size{3};
- const LagrangeDiscretization disc;
+ std::shared_ptr<const LagrangeDiscretization> disc;
bool withPressure = true;
@@ -23,9 +23,13 @@ protected:
public:
LagrangeElasticity(ElasticityProblem &eP, int degree, bool isStatic = false);
- const IDiscretization &GetDisc() const override { return disc; }
+ const IDiscretization &GetDisc() const override { return *disc; }
- std::unique_ptr<IDiscretization> CreateDisc(int sizes) const override {return std::make_unique<LagrangeDiscretization>(disc, degree,sizes);}
+ std::shared_ptr<IDiscretization> CreateDisc(int sizes) const override {
+ return std::make_shared<LagrangeDiscretization>(*disc, degree, sizes);
+ }
+
+ std::shared_ptr<const IDiscretization> GetSharedDisc() const override {return disc; }
void Initialize(Vectors &vecs) const override;
@@ -138,8 +142,9 @@ public:
std::pair<double, double> detF(const Vector &u) const override;
void Scale (double s) const {
+ //TODO warum konvergiert der LandBenchmark mit mehreren PressureSteps nicht, wenn einkommentiert?
Material &cellMat = eProblem.GetMaterial();
- cellMat.ScalePenalty(s*s);
+// cellMat.ScalePenalty(s*s);
}
};
diff --git a/src/elasticity/materials/GuccioneMaterial.cpp b/src/elasticity/materials/GuccioneMaterial.cpp
index 1eff5ec2ebfdaa42e32e87b5ca93d017dc4e948d..bb7c6127277e5a2cf4ad19a5e80e235d3898fb94 100644
--- a/src/elasticity/materials/GuccioneMaterial.cpp
+++ b/src/elasticity/materials/GuccioneMaterial.cpp
@@ -140,22 +140,20 @@ Tensor GuccioneBonetMaterial::ddq(const Tensor &E, const VectorField &f, const T
return 2*QE;
}
-
double GuccioneBonetMaterial::energy(const Tensor &E, const VectorField &f) const {
- return 0.5 * C() * (mat::cutexp(q(E, f)) - 1);
+ return 0.5 * C() * (mat::cutexp(q(E, f),exp_max) - 1);
}
Tensor GuccioneBonetMaterial::firstPK(const Tensor &E, const VectorField &f) const {
- Tensor DPsi = 0.5 * C() * mat::cutexp(q(E, f)) * dq(E, f) + activeStress * Product(f,f);
+ Tensor DPsi = 0.5 * C() * mat::cutexp(q(E, f),exp_max) * dq(E, f) + activeStress * Product(f,f);
return DPsi;
}
Tensor
GuccioneBonetMaterial::constitutiveMatrix(const Tensor &E, const VectorField &f,
const Tensor &H) const {
- double expQ = mat::cutexp(q(E, f));
-
+ double expQ = mat::cutexp(q(E, f),exp_max);
Tensor DDPsi = Frobenius(dq(E, f), H) * dq(E, f);
DDPsi += ddq(E, f, H);
return 0.5 * C() * expQ * DDPsi;
diff --git a/src/elasticity/materials/GuccioneMaterial.hpp b/src/elasticity/materials/GuccioneMaterial.hpp
index 174b82bc8aba1e01f56411054344a6e87f12a61a..626b2bbb82191cb37c86c2219aad9b9f25fa9c33 100644
--- a/src/elasticity/materials/GuccioneMaterial.hpp
+++ b/src/elasticity/materials/GuccioneMaterial.hpp
@@ -90,7 +90,8 @@ class GuccioneBonetMaterial : public HyperelasticMaterial {
double b{1.0};
double c{1.0};
double activeStress{0.0};
-
+ double exp_max = 10;
+
void calculateParams() {
a = 0.25 * (bf() - 2.0 * bfs() + bs());
b = 0.5 * (bfs() - bs());
@@ -130,18 +131,21 @@ public:
Config::Get("GuccioneMat_bf", parameters[1]);
Config::Get("GuccioneMat_bs", parameters[2]);
Config::Get("GuccioneMat_bfs", parameters[3]);
+ Config::Get("exp_max", exp_max);
calculateParams();
Config::Get("activeStress", activeStress);
}
explicit GuccioneBonetMaterial(const std::vector<double> ¶ms) : HyperelasticMaterial("Bonet",
params) {
+ Config::Get("exp_max", exp_max);
calculateParams();
Config::Get("activeStress", activeStress);
};
explicit GuccioneBonetMaterial(std::vector<double> &¶ms) : HyperelasticMaterial("Bonet",
std::move(params)) {
+ Config::Get("exp_max", exp_max);
calculateParams();
Config::Get("activeStress", activeStress);
};
diff --git a/src/elasticity/materials/HolzapfelMaterial.cpp b/src/elasticity/materials/HolzapfelMaterial.cpp
index 2cdfcad1df8c7f5c5bf18c725c526cc3e6b28d8a..79374d162dd5c90b474492d4a7ceef8fcbcd87df 100644
--- a/src/elasticity/materials/HolzapfelMaterial.cpp
+++ b/src/elasticity/materials/HolzapfelMaterial.cpp
@@ -2,44 +2,44 @@
double HolzapfelMaterial::energyIso(const Tensor &C) const {
double I_1 = trace(C);
- return 0.5 * a() / b() * (mat::cutexp(b() * (I_1 - 3)) - 1);
+ return 0.5 * a() / b() * (mat::cutexp(b() * (I_1 - 3),exp_max) - 1);
}
Tensor HolzapfelMaterial::firstPKIso(const Tensor &C) const {
double I_1 = trace(C);
- Tensor Psi = 0.5 * a() * (mat::cutexp(b() * (I_1 - 3))) * One;
+ Tensor Psi = 0.5 * a() * (mat::cutexp(b() * (I_1 - 3),exp_max)) * One;
return Psi;
}
Tensor
HolzapfelMaterial::constitutiveMatrixIso(const Tensor &C, const Tensor &H) const {
double I_1 = trace(C);
- Tensor Psi = 0.5 * a() * b()* (mat::cutexp(b() * (I_1 - 3))) * trace(H) * One;
+ Tensor Psi = 0.5 * a() * b()* (mat::cutexp(b() * (I_1 - 3),exp_max)) * trace(H) * One;
return Psi;
}
double HolzapfelMaterial::energyAniso(const Tensor &C, const Tensor &Q) const {
double I_4f = macaulay(Q[0] * (C * Q[0]) - 1.0);
- double energy = af() / (2.0 * bf()) * (mat::cutexp(bf() * I_4f * I_4f) - 1);
+ double energy = af() / (2.0 * bf()) * (mat::cutexp(bf() * I_4f * I_4f,exp_max) - 1);
double I_4s = macaulay(Q[1] * (C * Q[1]) - 1.0);
- energy += as() / (2.0 * bs()) * (mat::cutexp(bs() * I_4s * I_4s) - 1);
+ energy += as() / (2.0 * bs()) * (mat::cutexp(bs() * I_4s * I_4s,exp_max) - 1);
double I_8fs = Q[0] * (C * Q[1]);
- energy += afs() / (2.0 * bfs()) * (mat::cutexp(bfs() * I_8fs * I_8fs) - 1);
+ energy += afs() / (2.0 * bfs()) * (mat::cutexp(bfs() * I_8fs * I_8fs,exp_max) - 1);
return energy;
}
Tensor HolzapfelMaterial::firstPKAniso(const Tensor &C, const Tensor &Q) const {
double I_4f = macaulay(Q[0] * (C * Q[0]) - 1.0);
- Tensor Psi = af() * mat::cutexp(bf() * I_4f * I_4f) * I_4f * Product(Q[0], Q[0]);
+ Tensor Psi = af() * mat::cutexp(bf() * I_4f * I_4f,exp_max) * I_4f * Product(Q[0], Q[0]);
double I_4s = macaulay(Q[1] * (C * Q[1]) - 1.0);
- Psi += as() * mat::cutexp(bs() * I_4s * I_4s) * I_4s * Product(Q[1], Q[1]);
+ Psi += as() * mat::cutexp(bs() * I_4s * I_4s,exp_max) * I_4s * Product(Q[1], Q[1]);
double I_8fs = Q[0] * (C * Q[1]);
- Psi += afs() * mat::cutexp(bfs() * I_8fs * I_8fs) * I_8fs * sym(Product(Q[0], Q[1]));
+ Psi += afs() * mat::cutexp(bfs() * I_8fs * I_8fs,exp_max) * I_8fs * sym(Product(Q[0], Q[1]));
return Psi;
}
@@ -48,15 +48,15 @@ Tensor
HolzapfelMaterial::constitutiveMatrixAniso(const Tensor &C, const Tensor &Q, const Tensor &H) const {
double I_4f = macaulay(Q[0] * (C * Q[0]) - 1.0);
double Hf = Q[0] * (H * Q[0]);
- Tensor Psi = af() * mat::cutexp(bf() * I_4f * I_4f) * Hf * (2 * bf() * I_4f * I_4f + 1) * Product(Q[0], Q[0]);
+ Tensor Psi = af() * mat::cutexp(bf() * I_4f * I_4f,exp_max) * Hf * (2 * bf() * I_4f * I_4f + 1) * Product(Q[0], Q[0]);
double I_4s = macaulay(Q[1] * (C * Q[1]) - 1.0);
double Hs = Q[1] * (H * Q[1]);
- Psi += as() * mat::cutexp(bs() * I_4s * I_4s) * Hs * (2 * bs() * I_4s * I_4s + 1) * Product(Q[1], Q[1]);
+ Psi += as() * mat::cutexp(bs() * I_4s * I_4s,exp_max) * Hs * (2 * bs() * I_4s * I_4s + 1) * Product(Q[1], Q[1]);
double I_8fs = Q[0] * (C * Q[1]);
double Hfs = 0.5 * (Q[0] * (H * Q[1]) + Q[1] * (H * Q[0]));
- Psi += afs() * mat::cutexp(bfs() * I_8fs * I_8fs) * Hfs * (2 * bfs() * I_8fs * I_8fs + 1) * sym(Product(Q[0], Q[1]));
+ Psi += afs() * mat::cutexp(bfs() * I_8fs * I_8fs,exp_max) * Hfs * (2 * bfs() * I_8fs * I_8fs + 1) * sym(Product(Q[0], Q[1]));
return Psi;
}
diff --git a/src/elasticity/materials/HolzapfelMaterial.hpp b/src/elasticity/materials/HolzapfelMaterial.hpp
index 77117a8349863375e51c942d95a573124882b7bf..f46d801815ff9234d27eeb160e2d77c934e3325b 100644
--- a/src/elasticity/materials/HolzapfelMaterial.hpp
+++ b/src/elasticity/materials/HolzapfelMaterial.hpp
@@ -10,6 +10,8 @@
* See i.e. doi:10.1016/j.ijnonlinmec.2006.02.001
*/
class HolzapfelMaterial : public HyperelasticMaterial {
+ double exp_max = 10;
+
double a() const { return GetParameter(0); }
double af() const { return GetParameter(1); }
@@ -57,12 +59,17 @@ public:
Config::Get("HolzapfelMat_bf", parameters[5]);
Config::Get("HolzapfelMat_bs", parameters[6]);
Config::Get("HolzapfelMat_bfs", parameters[7]);
+ Config::Get("exp_max", exp_max);
}
- explicit HolzapfelMaterial(const std::vector<double> ¶ms) : HyperelasticMaterial("Holzapfel", params) {};
+ explicit HolzapfelMaterial(const std::vector<double> ¶ms) : HyperelasticMaterial("Holzapfel", params) {
+ Config::Get("exp_max", exp_max);
+ };
explicit HolzapfelMaterial(std::vector<double> &¶ms) : HyperelasticMaterial("Holzapfel",
- std::move(params)) {};
+ std::move(params)) {
+ Config::Get("exp_max", exp_max);
+ };
Tensor FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const Tensor &Q) const override;
diff --git a/src/elasticity/materials/Material.hpp b/src/elasticity/materials/Material.hpp
index ba618cbf11407d70bdfd672e63a4a7f9b9a16b21..cdb1f6a517a7b73da6841af0fab689c505981423 100644
--- a/src/elasticity/materials/Material.hpp
+++ b/src/elasticity/materials/Material.hpp
@@ -24,6 +24,9 @@ public:
HyperelasticMaterial(std::string &&name, std::vector<double> &¶ms) noexcept: name(std::move(name)), parameters(
std::move(params)) {}
+ virtual double q(const Tensor &E, const VectorField &f) const { return 0; }
+
+
/**
* Returns the strain energy of the material given the deformation Gradient F
* @param F = (Du + One) for the derivative Du of a given displacement u.
@@ -165,7 +168,10 @@ public:
IsoSecondDerivative(const Tensor &F, const T &H, const T &G) const {
return hyperelMat->IsotropicSecondDerivative(F, H, G);
}
-
+
+ double q(const Tensor &E, const VectorField &f) const {
+ return hyperelMat->q(E,f);
+ }
double AnisoEnergy(const Tensor &F, const Tensor &Q) const {
return hyperelMat->AnisotropicEnergy(F, Q);
diff --git a/src/elasticity/materials/VolumetricPenalty.hpp b/src/elasticity/materials/VolumetricPenalty.hpp
index 57bf5699757fb1a1f6cf1fb64bd699be4e9f6e67..fa70c96b14cbe7fb032a076698f01496ebc5426c 100644
--- a/src/elasticity/materials/VolumetricPenalty.hpp
+++ b/src/elasticity/materials/VolumetricPenalty.hpp
@@ -45,14 +45,22 @@ public:
CiarletPenalty(double p) : VolumetricPenalty(p) {}
- CiarletPenalty(double p, double m) : VolumetricPenalty(p), mu(p > 0.0 ? m / p : 1.0) {
- }
+ CiarletPenalty(double p, double m) : VolumetricPenalty(p), mu(p > 0.0 ? m / p : 1.0) {}
- double Functional(double J) const override { return 0.5 * (J - 1.0) * (J - 1.0) - mu * log(J); };
+ double Functional(double J) const override {
+ if (J <= 0) {
+ return 0;
+ }
+ return 0.5 * (J - 1.0) * (J - 1.0) - mu * log(J);
+ }
- double Derivative(double J) const override {return (J - 1.0) - mu / J; };
+ double Derivative(double J) const override {
+ return (J - 1.0) - mu / J;
+ }
- double SecondDerivative(double J) const override { return 1.0 + mu / (J * J); };
+ double SecondDerivative(double J) const override {
+ return 1.0 + mu / (J * J);
+ }
};
class LogarithmicCompressiblePenalty : public VolumetricPenalty {
@@ -61,7 +69,8 @@ public:
LogarithmicCompressiblePenalty(double p) : VolumetricPenalty(p) {}
- double Functional(double J) const override { return 0.5 * (J - 1.0) * log(J); };
+ double Functional(double J) const override { return 0.5 * (J - 1.0) * log(J);
+ }
double Derivative(double J) const override { return 0.5 * (log(J) + (J - 1.0) / J); };
diff --git a/src/elasticity/materials/materials.hpp b/src/elasticity/materials/materials.hpp
index 3db20c3e59ec52809c3b6d173b057e5c67eec095..00611e6bd923cb0cd50861ed81f22d46ac04678c 100644
--- a/src/elasticity/materials/materials.hpp
+++ b/src/elasticity/materials/materials.hpp
@@ -1,12 +1,12 @@
#ifndef MATERIALS_HPP
#define MATERIALS_HPP
-
#include "Tensor.hpp"
namespace mat {
- constexpr double cutexp(double argument, double max = 35.0) {
- return exp(min(argument, max));
+ constexpr double cutexp(double argument, double max = 35) {
+ if (argument<max) return exp(argument);
+ return exp(max) * (1 + argument-max);
}
static Tensor cofactor(const Tensor &F) {
diff --git a/src/elasticity/problem/ElasticityProblem.cpp b/src/elasticity/problem/ElasticityProblem.cpp
index e812f9aa13d33b93ced7d7e6c738fe82b6cd09d4..9be4e99f3be0888b2dcfb51f1eb90fa351875533 100644
--- a/src/elasticity/problem/ElasticityProblem.cpp
+++ b/src/elasticity/problem/ElasticityProblem.cpp
@@ -8,9 +8,8 @@
std::vector<double> findDisplacements(const Vector &u, const vector<Point> &points, int index) {
std::vector<double> displacements(points.size());
- std::string disc = "";
- Config::Get("MechDiscretization",disc);
- if (contains(disc,"Conforming") || contains(disc,"EG")) {
+
+ if (!contains(u.GetDisc().DiscName(),"DG")) {
for (row r = u.rows(); r != u.rows_end(); ++r) {
auto [i, isNear] = isIn(r(), points);
if (isNear) {
@@ -26,7 +25,7 @@ std::vector<double> findDisplacements(const Vector &u, const vector<Point> &poin
displacements[i] = PPM->Sum(displacements[i]);
}
}
- else if (contains(disc,"DG")) {
+ else {
std::vector<VectorField> disp(points.size());
disp = findDisplacements(u,points);
for (int l = 0; l < displacements.size(); ++l) {
@@ -38,9 +37,7 @@ std::vector<double> findDisplacements(const Vector &u, const vector<Point> &poin
std::vector<VectorField> findDisplacements(const Vector &u, const vector<Point> &points) {
std::vector<VectorField> displacements(points.size());
- std::string disc = "";
- Config::Get("MechDiscretization",disc);
- if (contains(disc,"Conforming") || contains(disc,"EG")) {
+ if (!contains(u.GetDisc().DiscName(),"DG")) {
for (row r = u.rows(); r != u.rows_end(); ++r) {
auto [i, isNear] = isIn(r(), points);
if (isNear) {
@@ -58,7 +55,7 @@ std::vector<VectorField> findDisplacements(const Vector &u, const vector<Point>
displacements[i] = PPM->Sum(displacements[i]);
}
}
- else if (contains(disc,"DG")) {
+ else {
std::vector<int> cellCount(points.size());
for (cell c = u.cells(); c != u.cells_end(); ++c) {
@@ -84,7 +81,7 @@ std::vector<VectorField> findDisplacements(const Vector &u, const vector<Point>
}
-ElasticityProblem::ElasticityProblem() : IProblem(), IElasticityProblem(), IPressureProblem() {
+ElasticityProblem::ElasticityProblem() : CardMechIProblem(), IElasticityProblem(), IPressureProblem() {
Config::Get("MechProblemVerbose", verbose);
std::string aMatName{"Linear"};
@@ -98,7 +95,7 @@ ElasticityProblem::ElasticityProblem() : IProblem(), IElasticityProblem(), IPres
ElasticityProblem::ElasticityProblem(const std::string &matName,
const std::vector<double> &matParameters)
- : IProblem(), IElasticityProblem(), IPressureProblem() {
+ : CardMechIProblem(), IElasticityProblem(), IPressureProblem() {
Config::Get("MechProblemVerbose", verbose);
activeMaterial = Material(matName, matParameters);
passiveMaterial = Material(matName, matParameters);
diff --git a/src/elasticity/problem/ElasticityProblem.hpp b/src/elasticity/problem/ElasticityProblem.hpp
index 21b0351f702fe094d887148099090206dbf06f1f..14e32b3ad53d3eb047a69982b5b8d9095e79c055 100644
--- a/src/elasticity/problem/ElasticityProblem.hpp
+++ b/src/elasticity/problem/ElasticityProblem.hpp
@@ -1,22 +1,22 @@
#ifndef ELASTICITYPROBLEM_HPP
#define ELASTICITYPROBLEM_HPP
-#include "IProblem.hpp"
+#include "CardMechIProblem.hpp"
#include "Vector.hpp"
#include "CardiacProblem.hpp"
#include "materials/Material.hpp"
namespace mpp {
- template<class T, class Alloc, class Pred>
- constexpr typename std::vector<T, Alloc>::size_type
- erase_if(std::vector<T, Alloc> &c, Pred pred) {
- auto it = std::remove_if(c.begin(), c.end(), pred);
- auto r = std::distance(it, c.end());
- c.erase(it, c.end());
- return r;
- }
+template <class T, class Alloc, class Pred>
+constexpr typename std::vector<T, Alloc>::size_type
+erase_if(std::vector<T, Alloc> &c, Pred pred) {
+ auto it = std::remove_if(c.begin(), c.end(), pred);
+ auto r = std::distance(it, c.end());
+ c.erase(it, c.end());
+ return r;
}
+} // namespace mpp
std::array<double, 4> pressureFromConfig();
double pressureFromConfig(short chamber);
@@ -24,14 +24,15 @@ double pressureFromConfig(short chamber);
std::vector<double>
findDisplacements(const Vector &u, const std::vector<Point> &points, int index);
-std::vector<VectorField>
-findDisplacements(const Vector &u, const std::vector<Point> &points);
+std::vector<VectorField> findDisplacements(const Vector &u,
+ const std::vector<Point> &points);
class IElasticityProblem {
bool useFlorySplit{false};
bool isStatic{false};
std::string elasticityModel{"Conforming"};
+
public:
IElasticityProblem() {
Config::Get("MechDiscretization", elasticityModel);
@@ -40,121 +41,176 @@ public:
Config::Get("VolumetricSplit", useFlorySplit);
}
- [[nodiscard]] virtual const std::string &Model() const { return elasticityModel; }
+ [[nodiscard]] virtual const std::string &Model() const {
+ return elasticityModel;
+ }
+
+ virtual const void SetModel(std::string modelName) {
+ elasticityModel = modelName;
+ }
[[nodiscard]] virtual bool IsStatic() const { return isStatic; }
- [[nodiscard]] virtual VectorField
- Load(double t, const Point &x, const Cell &c) const { return zero; };
+ [[nodiscard]] virtual VectorField Load(double t, const Point &x,
+ const Cell &c) const {
+ return zero;
+ };
- [[nodiscard]] virtual double ActiveStress(double t, const Point &x) const { return 0.0; };
+ [[nodiscard]] virtual double ActiveStress(double t, const Point &x) const {
+ return 0.0;
+ };
- [[nodiscard]] virtual double ActiveStretch(double t, const Point &x) const { return 0.0; };
+ [[nodiscard]] virtual double ActiveStretch(double t, const Point &x) const {
+ return 0.0;
+ };
- [[nodiscard]] virtual double
- TractionStiffness(double t, const Point &x, int sd) const { return 0.0; };
+ [[nodiscard]] virtual double TractionStiffness(double t, const Point &x,
+ int sd) const {
+ return 0.0;
+ };
- [[nodiscard]] virtual double
- TractionViscosity(double t, const Point &x, int sd) const { return 0.0; };
+ [[nodiscard]] virtual double TractionViscosity(double t, const Point &x,
+ int sd) const {
+ return 0.0;
+ };
[[nodiscard]] Tensor IsochoricPart(const Tensor &F) const {
return useFlorySplit ? pow(det(F), -1.0 / 3.0) * F : F;
}
- [[nodiscard]] virtual VectorField Deformation(double time, const Point &x) const { return zero; };
+ [[nodiscard]] virtual VectorField Deformation(double time,
+ const Point &x) const {
+ return zero;
+ };
- [[nodiscard]] virtual VectorField Velocity(double time, const Point &x) const { return zero; };
+ [[nodiscard]] virtual VectorField Velocity(double time,
+ const Point &x) const {
+ return zero;
+ };
- [[nodiscard]] virtual VectorField
- Acceleration(double time, const Point &x) const { return zero; };
+ [[nodiscard]] virtual VectorField Acceleration(double time,
+ const Point &x) const {
+ return zero;
+ };
- [[nodiscard]] virtual Tensor
- DeformationGradient(double time, const Point &x) const { return Zero; };
+ [[nodiscard]] virtual Tensor DeformationGradient(double time,
+ const Point &x) const {
+ return Zero;
+ };
};
class IPressureProblem {
mutable double pressureScale{1.0};
+
protected:
- virtual double pressure(double t, const Point &x, int bndCondition) const { return 0.0; }
+ virtual double pressure(double t, const Point &x, int bndCondition) const {
+ return 0.0;
+ }
+
public:
IPressureProblem() = default;
- [[nodiscard]] virtual double Scale() const{
- return pressureScale;
- }
+ [[nodiscard]] virtual double Scale() const { return pressureScale; }
virtual void Scale(double s) const {
- pressureScale = ((s>1.0) + (0.0<s<=1.0) * s);
+ pressureScale = ((s > 1.0) + (0.0 < s <= 1.0) * s);
}
virtual double Pressure(double t, const Point &x, int bndCondition) const {
return Scale() * pressure(t, x, bndCondition);
}
- virtual double InitialPressure(const Point &x, int bndCondition) const { return 0.0; }
+ virtual double InitialPressure(const Point &x, int bndCondition) const {
+ return 0.0;
+ }
- virtual void Initialize(double dt, double t, const std::array<double, 4> &mechVolumes) {}
+ virtual void Initialize(double dt, double t,
+ const std::array<double, 4> &mechVolumes) {}
- virtual void Update(double dt, double t, const std::array<double, 4> &mechVolumes) {}
+ virtual void Update(double dt, double t,
+ const std::array<double, 4> &mechVolumes) {}
- virtual void Finalize(double dt, double t, const std::array<double, 4> &mechVolumes) {}
+ virtual void Finalize(double dt, double t,
+ const std::array<double, 4> &mechVolumes) {}
- virtual void Iterate(double dt, double t, const std::array<double, 4> &mechVolumes) {}
+ virtual void Iterate(double dt, double t,
+ const std::array<double, 4> &mechVolumes) {}
- virtual bool Converged(const std::array<double, 4> &mechVolumes) { return true; }
+ virtual bool Converged(const std::array<double, 4> &mechVolumes) {
+ return true;
+ }
};
-
-class ElasticityProblem : public IProblem, public IElasticityProblem, public IPressureProblem {
+class ElasticityProblem : public CardMechIProblem,
+ public IElasticityProblem,
+ public IPressureProblem {
protected:
std::string elasticityModel{"Conforming"};
Material activeMaterial;
Material passiveMaterial;
+ double DGpenalty = 0;
+
bool meshFromConfig() const { return meshesName != ""; }
+
public:
ElasticityProblem();
- ElasticityProblem(const std::string& matName, const std::vector<double>& matParameters);
+ ElasticityProblem(const std::string &matName,
+ const std::vector<double> &matParameters);
[[nodiscard]] const string &GetMeshName() const { return meshesName; }
+ virtual const void SetMeshName(std::string meshName) {
+ meshesName = meshName;
+ }
[[nodiscard]] const string &GetMaterialName(bool active) const {
return (active ? activeMaterial.Name() : passiveMaterial.Name());
};
[[nodiscard]] const Material &GetMaterial(const cell &c) const {
- //TODO: enable different materials
+ // TODO: enable different materials
bool isActive = true;
return (isActive ? activeMaterial : passiveMaterial);
}
- [[nodiscard]] Material &GetMaterial() {
- return activeMaterial;
- }
+ [[nodiscard]] Material &GetMaterial() { return activeMaterial; }
- [[nodiscard]] std::pair<const Material &, const Material &> GetMaterials() const {
- return std::pair<const Material &, const Material &>{activeMaterial, passiveMaterial};
+ [[nodiscard]] std::pair<const Material &, const Material &>
+ GetMaterials() const {
+ return std::pair<const Material &, const Material &>{activeMaterial,
+ passiveMaterial};
}
- [[nodiscard]] virtual std::string Evaluate(const Vector &solution) const { return ""; };
+ [[nodiscard]] virtual std::string Evaluate(const Vector &solution) const {
+ return "";
+ };
- [[nodiscard]] virtual std::vector<std::string> EvaluationMetrics() const { return {}; }
+ [[nodiscard]] virtual std::vector<std::string> EvaluationMetrics() const {
+ return {};
+ }
[[nodiscard]] virtual std::vector<double>
- EvaluationResults(const Vector &solution) const { return std::vector<double>{}; }
+ EvaluationResults(const Vector &solution) const {
+ return std::vector<double>{};
+ }
-};
+ [[nodiscard]] virtual const double &GetDGPenalty() const { return DGpenalty; }
+ virtual const void SetDGPenalty(double p) { DGpenalty = p; }
+};
class DefaultElasticityProblem : public ElasticityProblem {
std::array<double, 4> chamberPressure{0.0, 0.0, 0.0, 0.0};
+
protected:
- [[nodiscard]] double pressure(double t, const Point &x, int bndCondition) const override {
+ [[nodiscard]] double pressure(double t, const Point &x,
+ int bndCondition) const override {
return chamberPressure[bndCondition % 230];
}
+
public:
DefaultElasticityProblem()
: ElasticityProblem(), chamberPressure(pressureFromConfig()) {}
@@ -165,25 +221,26 @@ public:
explicit DefaultElasticityProblem(std::array<double, 4> pressure)
: ElasticityProblem(), chamberPressure(pressure) {}
-
[[nodiscard]] std::string Name() const override {
return "Default Elasticity Problem";
}
};
-
class InitialPrestressProblem : public ElasticityProblem {
const ElasticityProblem &refProblem;
+
protected:
- [[nodiscard]] double pressure(double t, const Point &x, int bndCondition) const override {
+ [[nodiscard]] double pressure(double t, const Point &x,
+ int bndCondition) const override {
return refProblem.Pressure(0.0, x, bndCondition);
}
-public:
- InitialPrestressProblem(const ElasticityProblem &reference) : refProblem(reference) {
- }
+public:
+ InitialPrestressProblem(const ElasticityProblem &reference)
+ : refProblem(reference) {}
- [[nodiscard]] VectorField Load(double t, const Point &x, const Cell &c) const override {
+ [[nodiscard]] VectorField Load(double t, const Point &x,
+ const Cell &c) const override {
return refProblem.Load(0.0, x, c);
}
@@ -191,27 +248,33 @@ public:
return refProblem.ActiveStretch(0.0, x);
}
- [[nodiscard]] double TractionStiffness(double t, const Point &x, int sd) const override {
+ [[nodiscard]] double TractionStiffness(double t, const Point &x,
+ int sd) const override {
return refProblem.TractionStiffness(t, x, sd);
}
- [[nodiscard]] double TractionViscosity(double t, const Point &x, int sd) const override {
+ [[nodiscard]] double TractionViscosity(double t, const Point &x,
+ int sd) const override {
return refProblem.TractionViscosity(t, x, sd);
}
- [[nodiscard]] VectorField Deformation(double time, const Point &x) const override {
+ [[nodiscard]] VectorField Deformation(double time,
+ const Point &x) const override {
return refProblem.Deformation(0.0, x);
}
- [[nodiscard]] VectorField Velocity(double time, const Point &x) const override {
+ [[nodiscard]] VectorField Velocity(double time,
+ const Point &x) const override {
return refProblem.Velocity(0.0, x);
}
- [[nodiscard]] VectorField Acceleration(double time, const Point &x) const override {
+ [[nodiscard]] VectorField Acceleration(double time,
+ const Point &x) const override {
return refProblem.Acceleration(0.0, x);
}
- [[nodiscard]] Tensor DeformationGradient(double time, const Point &x) const override {
+ [[nodiscard]] Tensor DeformationGradient(double time,
+ const Point &x) const override {
return refProblem.DeformationGradient(0.0, x);
}
@@ -219,10 +282,7 @@ public:
return refProblem.Model();
}
- [[nodiscard]] bool IsStatic() const override {
- return refProblem.IsStatic();
- }
-
+ [[nodiscard]] bool IsStatic() const override { return refProblem.IsStatic(); }
[[nodiscard]] string Name() const override {
return "Initial Prestress Problem";
@@ -231,6 +291,7 @@ public:
std::unique_ptr<ElasticityProblem> GetElasticityProblem();
-std::unique_ptr<ElasticityProblem> GetElasticityProblem(const std::string &problemName);
+std::unique_ptr<ElasticityProblem>
+GetElasticityProblem(const std::string &problemName);
-#endif //ELASTICITYPROBLEM_HPP
+#endif // ELASTICITYPROBLEM_HPP
diff --git a/src/elasticity/solvers/IterativePressureSolver.cpp b/src/elasticity/solvers/IterativePressureSolver.cpp
index a7c2514129367e6b908e19212e7f52af2654bca9..5b0cd9d576314c00f19a222a802dab2b55b77a3a 100644
--- a/src/elasticity/solvers/IterativePressureSolver.cpp
+++ b/src/elasticity/solvers/IterativePressureSolver.cpp
@@ -4,7 +4,6 @@
#include <chrono>
#include "LagrangeTransfer.hpp"
-
void IterativePressureSolver::Initialize(IElasticity &assemble, Vector &u) {
problem.Scale(0.0);
assemble.SetInitialValue(u);
@@ -24,11 +23,8 @@ bool IterativePressureSolver::Method(IElasticity &assemble, Vector &u) {
Vector uNew(u);
Initialize(assemble, uNew);
-
while (!iteration.IsFinished()) {
- //for loop overlevels:
- //{ u0 u1,....
- bool converged = Step(assemble, uNew);
+ bool converged = Step(assemble, uNew);
if(converged) {
//mapping von corse to finer
@@ -46,7 +42,6 @@ bool IterativePressureSolver::Method(IElasticity &assemble, Vector &u) {
return Method(assemble, u);
return false;
}
- //}
}
mout.EndBlock();
return true;
diff --git a/src/electrophysiology/MainMonodomain.cpp b/src/electrophysiology/MainMonodomain.cpp
index 402509dfc6369aecd1edb969345768709329593b..3213cc0e07b241ece5bd1032d57afd59980ff6cd 100644
--- a/src/electrophysiology/MainMonodomain.cpp
+++ b/src/electrophysiology/MainMonodomain.cpp
@@ -41,7 +41,7 @@ void MainMonodomain::Initialize() {
elphyA->ResetTime(times[0], times[1], times[2]);
// Initialize Vectors
- potential = std::make_unique<Vector>(0.0, elphyA->GetDisc());
+ potential = std::make_unique<Vector>(0.0, elphyA->GetSharedDisc());
problem->SetEvaluationCells(*potential);
bool printMeshOnly = false;
diff --git a/src/electrophysiology/MainOneCell.cpp b/src/electrophysiology/MainOneCell.cpp
index 58a5811d80c0a55c9fe6238d3dc37152e31d9158..700ede9d1265dcf2017eb67a7c7f6a5865970617 100644
--- a/src/electrophysiology/MainOneCell.cpp
+++ b/src/electrophysiology/MainOneCell.cpp
@@ -23,7 +23,7 @@ Vector &MainOneCell::Run() {
auto cellModel = GetElphyModel();
elphyA = CreateElphyAssemble(*problem, *cellModel, 0);
elphyA->ResetTime(times[0], times[1], times[2]);
- potential = std::make_unique<Vector>(cellModel->InitialValue(ELPHY), elphyA->GetDisc());
+ potential = std::make_unique<Vector>(cellModel->InitialValue(ELPHY), elphyA->GetSharedDisc());
if (problem->Name() == "EigenvalueProblem") {
MElphySolver solver(*cellModel, elphyA->ExcitationData(), "Eigenvalue");
solver.SetVerbose(15);
diff --git a/src/electrophysiology/assemble/IElphyAssemble.cpp b/src/electrophysiology/assemble/IElphyAssemble.cpp
index aee71d85046bfd963898d17e77065432b5238403..e49e36a2affa8b38edb70610c43dee5ff2ea8758 100644
--- a/src/electrophysiology/assemble/IElphyAssemble.cpp
+++ b/src/electrophysiology/assemble/IElphyAssemble.cpp
@@ -13,9 +13,9 @@
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),
- dataDisc(elphyProblem.GetMeshes(), degree, exactUpTo < 0 ? 2 * degree : exactUpTo, DomainPart,
- 3),
+ disc(std::make_shared<const MultiPartDiscretization>(elphyProblem.GetMeshes(), degree, exactUpTo < 0 ? 2 * degree : exactUpTo, DomainPart, 1)),
+ dataDisc(std::make_shared<const MultiPartDiscretization>(elphyProblem.GetMeshes(), degree, exactUpTo < 0 ? 2 * degree : exactUpTo, DomainPart,
+ 3)),
excData(std::move(InterpolateMeshData(dataDisc))) {
Config::Get("ElphyVerbose", verbose);
Config::Get("IsExternalCurrentSmooth", iextContinuous);
diff --git a/src/electrophysiology/assemble/IElphyAssemble.hpp b/src/electrophysiology/assemble/IElphyAssemble.hpp
index b2862a77f0dee1eb3825ac5ea6e0ec0624f2a45d..68425c47874c1211e3444705e46afc2b3cce66eb 100644
--- a/src/electrophysiology/assemble/IElphyAssemble.hpp
+++ b/src/electrophysiology/assemble/IElphyAssemble.hpp
@@ -25,8 +25,8 @@ protected:
//Discretization
int degree{1};
- const MultiPartDiscretization disc;
- const MultiPartDiscretization dataDisc;
+ std::shared_ptr<const MultiPartDiscretization> disc;
+ std::shared_ptr<const MultiPartDiscretization> dataDisc;
Vector excData;
int printVSteps{1};
@@ -36,12 +36,10 @@ public:
IElphyAssemble(ElphyProblem &elphyProblem, int degree, int exactUpTo = -1);
const IDiscretization &GetDisc() const {
- return disc;
+ return *disc;
}
- std::unique_ptr<IDiscretization> CreateDisc(int size) const {
- return std::make_unique<MultiPartDiscretization>(disc, size);
- }
+ std::shared_ptr<const IDiscretization> GetSharedDisc() const { return disc; }
virtual void Initialize(Vector &u) const {};
diff --git a/src/electrophysiology/assemble/Monodomain.cpp b/src/electrophysiology/assemble/Monodomain.cpp
index 948244a0600e799faa0c05645866200f32baf17e..3c2bb6a1a71e0005685175fbf69c5f26484f938a 100644
--- a/src/electrophysiology/assemble/Monodomain.cpp
+++ b/src/electrophysiology/assemble/Monodomain.cpp
@@ -196,8 +196,8 @@ void Monodomain::RHS(const cell &c, const Vectors &VCW, Vector &r) const {
MonodomainLINodes::MonodomainLINodes(ElphyProblem &elphyProblem, MCellModel &cellModel,
int degree, int exactUpTo) : Monodomain(elphyProblem, cellModel,
degree, exactUpTo) {
- IionVEC = std::make_unique<Vector>(0.0, GetDisc());
- DVIionVEC = std::make_unique<Vector>(0.0, GetDisc());
+ IionVEC = std::make_unique<Vector>(0.0, GetSharedDisc());
+ DVIionVEC = std::make_unique<Vector>(0.0, GetSharedDisc());
}
void MonodomainLINodes::updateIionVecs(const Vectors &VCW) {
@@ -340,8 +340,8 @@ MonodomainQuadrature::MonodomainQuadrature(ElphyProblem &elphyProblem, MCellMode
nJ = E.size();
}
- CWdisc = std::make_unique<DoFDiscretization<CellDoF>>(u.GetDisc().GetMeshes(), nQ * nJ);
- CWcells = std::make_unique<Vector>(*CWdisc);
+ CWdisc = std::make_shared<DoFDiscretization<CellDoF>>(u.GetDisc().GetMeshes(), nQ * nJ);
+ CWcells = std::make_unique<Vector>(CWdisc);
VCWQuadrature = std::make_unique<Vectors>(elphyModel.ElphySize(), *CWcells);
elphyModel.Initialize(*VCWQuadrature);
}
@@ -420,7 +420,7 @@ void MonodomainQuadrature::RHS(const cell &c, const Vectors &VCW, Vector &r) con
MonodomainFullNewton::MonodomainFullNewton(ElphyProblem &elphyProblem, MCellModel &cellModel,
int degree, int exactUpTo)
: Monodomain(elphyProblem, cellModel, degree, exactUpTo) {
- VACAW = std::make_unique<Vectors>(elphyModel.ElphySize(), GetDisc());
+ VACAW = std::make_unique<Vectors>(elphyModel.ElphySize(), GetSharedDisc());
elphyModel.Initialize(*VACAW);
/*if(MK==nullptr){
diff --git a/src/electrophysiology/assemble/Monodomain.hpp b/src/electrophysiology/assemble/Monodomain.hpp
index 561ac224ce00ce38dc1c3e28f3a59393f429da2d..3cd2c453b7e8bccc0ed8798e85ce0c52387a53d3 100644
--- a/src/electrophysiology/assemble/Monodomain.hpp
+++ b/src/electrophysiology/assemble/Monodomain.hpp
@@ -137,7 +137,7 @@ public:
class MonodomainQuadrature : public Monodomain {
std::unique_ptr<Vectors> VCWQuadrature;
std::unique_ptr<Vector> CWcells;
- std::unique_ptr<DoFDiscretization<CellDoF>> CWdisc;
+ std::shared_ptr<DoFDiscretization<CellDoF>> CWdisc;
public:
MonodomainQuadrature(ElphyProblem &elphyProblem, MCellModel &cellModel, const Vector &u,
int degree, int exactUpTo = -1);
@@ -199,7 +199,7 @@ public:
MonodomainSemiImplicitNodes(ElphyProblem &elphyProblem, MCellModel &cellModel, int degree,
int exactUpTo = -1)
: Monodomain(elphyProblem, cellModel, degree, exactUpTo) {
- IionVEC = std::make_unique<Vector>(0.0, GetDisc());
+ IionVEC = std::make_unique<Vector>(0.0, GetSharedDisc());
};
@@ -210,17 +210,15 @@ public:
};
class MonodomainSemiImplicitCells : public Monodomain {
protected:
- std::unique_ptr<LagrangeDiscretization> disc0;
+ std::shared_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);
-
+ disc0 = std::make_shared<LagrangeDiscretization >(mesh, 0, 1);
+ IionCell = std::make_unique<Vector>(0.0, disc0);;
};
@@ -230,7 +228,6 @@ public:
virtual void RHS(const cell &c, const Vector &u, Vector &r) const override;
void updateIionVecs(const Vectors &VCW) override;
-
};
class MonodomainSemiImplicitCellsIextPerCell : public MonodomainSemiImplicitCells {
std::unique_ptr<Vectors> externalCurrentOnCells;;
@@ -239,7 +236,7 @@ public:
MonodomainSemiImplicitCellsIextPerCell(ElphyProblem &elphyProblem, MCellModel &cellModel, int degree,
int exactUpTo = -1)
: MonodomainSemiImplicitCells(elphyProblem, cellModel, degree, exactUpTo) {
-
+
};
void RHS(const cell &c, const Vector &u, Vector &r) const override;
diff --git a/src/electrophysiology/problem/ElphyProblem.hpp b/src/electrophysiology/problem/ElphyProblem.hpp
index 3d85d757de0f639f94aea6adcc447f69bf5a069b..0c4bc23c4b28292a16c79d59ab13654a8bfa2a04 100644
--- a/src/electrophysiology/problem/ElphyProblem.hpp
+++ b/src/electrophysiology/problem/ElphyProblem.hpp
@@ -3,7 +3,7 @@
#include "utility/Config.hpp"
-#include "IProblem.hpp"
+#include "CardMechIProblem.hpp"
#include "CardiacProblem.hpp"
#include "CardiacData.hpp"
@@ -154,17 +154,17 @@ public:
const EvaluationPointData &getCellAt(int i) const { return cellsBesideMeshPoints.at(i); };
};
-class ElphyProblem : public IElphyProblem, public IProblem {
+class ElphyProblem : public IElphyProblem, public CardMechIProblem {
protected:
bool meshFromConfig() const { return meshesName != ""; }
public:
- ElphyProblem() : IProblem(), IElphyProblem() {}
+ ElphyProblem() : CardMechIProblem(), IElphyProblem() {}
- explicit ElphyProblem(bool at) : IProblem(), IElphyProblem(at) {}
+ explicit ElphyProblem(bool at) : CardMechIProblem(), IElphyProblem(at) {}
ElphyProblem(bool at, std::array<double, 4> sigmaParams)
- : IProblem(), IElphyProblem(at, sigmaParams) {}
+ : CardMechIProblem(), IElphyProblem(at, sigmaParams) {}
[[nodiscard]] const string &GetMeshName() const { return meshesName; }
diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.cpp b/src/electrophysiology/solvers/LinearImplicitSolver.cpp
index 6bfb193378799b3c18db2aa4c314a08008aa98dd..da62aad48f23645769864528160fb94a35d9a78c 100644
--- a/src/electrophysiology/solvers/LinearImplicitSolver.cpp
+++ b/src/electrophysiology/solvers/LinearImplicitSolver.cpp
@@ -31,8 +31,8 @@ void LinearImplicitSolver::Initialize(IElphyAssemble &A, Vector &V) {
if ((vcw_c) == nullptr)
{
const Meshes &mesh = (*potential).GetMeshes();
- disc0=std::make_unique<LagrangeDiscretization >(mesh, 0, 1);
- cellSizeV=std::make_unique<Vector>(*disc0);
+ disc0=std::make_shared<LagrangeDiscretization >(mesh, 0, 1);
+ cellSizeV=std::make_unique<Vector>(disc0);
}
vcw_c = std::make_unique<Vectors>(M.Size(), *cellSizeV);
//gating_c = std::make_unique<Vectors>(M.Size(), gamma_f_c);
@@ -401,7 +401,7 @@ void SemiImplicitSolverOnCells:: Initialize(IElphyAssemble &A, Vector &V){
{
const Meshes &mesh = (*potential).GetMeshes();
disc0=std::make_unique<LagrangeDiscretization >(mesh, 0, 1);
- cellSizeV=std::make_unique<Vector>(*disc0);
+ cellSizeV=std::make_unique<Vector>(disc0);
}
vcw_c = std::make_unique<Vectors>(M.Size(), *cellSizeV);
diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.hpp b/src/electrophysiology/solvers/LinearImplicitSolver.hpp
index 966d86250409e45d463e7d26181b8705cef6049b..02557a89f824e9eb624caaf1959fe9ba2c1f4e5b 100644
--- a/src/electrophysiology/solvers/LinearImplicitSolver.hpp
+++ b/src/electrophysiology/solvers/LinearImplicitSolver.hpp
@@ -8,7 +8,6 @@
#include <MCellModel.hpp>
#include "Newton.hpp"
-
void QuadraticIonStep(MCellModel &cellModel, double t, double dt, std::vector<double> &VW,
std::vector<std::vector<double>> &g);
@@ -49,7 +48,7 @@ protected:
int ionSteps{1};
bool plotCa{false};
std::unique_ptr<Vector> cellSizeV;
- std::unique_ptr<LagrangeDiscretization> disc0;
+ std::shared_ptr<LagrangeDiscretization> disc0;
std::unique_ptr<Vectors> externalCurrentOnCells;
void selectIonScheme(const string &solvingIonScheme);
@@ -123,7 +122,6 @@ public:
virtual void updateValues(Vectors &values, Vector &gamma_f_c){Exit("only defined for semi-implicit solver");};
void GatingVectorAti( Vector &values, int i) const override;
-
};
/*class LiniearImplicitSolverQuad: public LinearImplicitSolver{
@@ -159,11 +157,8 @@ public:
};
class SemiImplicitSolverOnCells : public LinearImplicitSolver {
- //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){};
diff --git a/test/cellmodels/TestElphyModels.cpp b/test/cellmodels/TestElphyModels.cpp
index 163adbe37f6e1d682b0a30e10c90dc7ca697d1d7..89276b673dc0a55e4b905a0e9b24b43d8e315d5c 100644
--- a/test/cellmodels/TestElphyModels.cpp
+++ b/test/cellmodels/TestElphyModels.cpp
@@ -82,7 +82,7 @@ public:
WithCellData(DataContainer({1, 0, 0, 0, 1, 0, 0, 0, 1})).
WithVertexData(DataContainer({0.0, 0.002, 1.0})).
CreateUnique();
- auto disc = LagrangeDiscretization(*meshes, 1);
+ auto disc = std::make_shared<const LagrangeDiscretization>(*meshes, 1);
dummyData = std::make_unique<Vector>(0.0, disc);
}
diff --git a/test/cellmodels/TestSolvers.cpp b/test/cellmodels/TestSolvers.cpp
index a9ace41ce27507c5b3743037dae49769ae92964a..6da3b67415ad46f732ad68fc60959d7577f75d9c 100644
--- a/test/cellmodels/TestSolvers.cpp
+++ b/test/cellmodels/TestSolvers.cpp
@@ -39,7 +39,7 @@ public:
WithCellData(DataContainer({1, 0, 0, 0, 1, 0, 0, 0, 1})).
WithVertexData(DataContainer({0.0, 0.002, 1.0})).
CreateUnique();
- auto disc = LagrangeDiscretization(*meshes, 0);
+ auto disc = std::make_shared<const LagrangeDiscretization>(*meshes, 0);
dummyData = std::make_unique<Vector>(0.0, disc);
}
diff --git a/test/cellmodels/TestTensionModel.cpp b/test/cellmodels/TestTensionModel.cpp
index eba0ad96bdc3d6b3b55b6ee6ec2261e6e11bc5e0..1152c7e165f1ebb0c4a1ae53b474bb7df395235b 100644
--- a/test/cellmodels/TestTensionModel.cpp
+++ b/test/cellmodels/TestTensionModel.cpp
@@ -26,7 +26,7 @@ public:
WithCellData(DataContainer({1, 0, 0, 0, 1, 0, 0, 0, 1})).
WithVertexData(DataContainer({0.0, 0.002, 1.0})).
CreateUnique();
- auto disc = LagrangeDiscretization(*meshes, 1);
+ auto disc = std::make_shared<const LagrangeDiscretization>(*meshes, 1);
dummyData = std::make_unique<Vector>(0.0, disc);
}
diff --git a/test/coupled/CMakeLists.txt b/test/coupled/CMakeLists.txt
index fd7ca81cf58b0eb21762883687a89c87839c88e4..b1121a9e04d8bda8a4e301bc4b0429a75072ab94 100644
--- a/test/coupled/CMakeLists.txt
+++ b/test/coupled/CMakeLists.txt
@@ -1,5 +1,3 @@
-add_mpp_test(TestCardiacTransfers COUPLED)
-add_mpp_test(TestDeformedDiffusion COUPLED)
-
add_mpi_test(TestCardiacTransfers COUPLED)
-add_mpi_test(TestCoupledConsistency COUPLED)
\ No newline at end of file
+add_mpi_test(TestCoupledConsistency COUPLED)
+add_mpi_test(TestDeformedDiffusion COUPLED)
diff --git a/test/coupled/TestCardiacTransfers.cpp b/test/coupled/TestCardiacTransfers.cpp
index 06312517bd07081438b147e8596d696267e61dfe..9ebcc6d54965528162d23d6fd32a30cdf2b90696 100644
--- a/test/coupled/TestCardiacTransfers.cpp
+++ b/test/coupled/TestCardiacTransfers.cpp
@@ -20,8 +20,8 @@ class TestCardiacTransfer : public TestWithParam<CTransferParam> {
protected:
std::unique_ptr<Meshes> coarseMeshes;
std::unique_ptr<Meshes> fineMeshes;
- std::unique_ptr<LagrangeDiscretization> coarseDisc;
- std::unique_ptr<MultiPartDiscretization> fineDisc;
+ std::shared_ptr<LagrangeDiscretization> coarseDisc;
+ std::shared_ptr<MultiPartDiscretization> fineDisc;
std::unique_ptr<Vector> coarse;
std::unique_ptr<Vector> fine;
@@ -37,12 +37,12 @@ protected:
WithIncludeVector(getIncludedChambers(GetParam().coarseChambers)).CreateUnique();
- coarseDisc = std::make_unique<LagrangeDiscretization>(*coarseMeshes, GetParam().coarseDegree,
+ coarseDisc = std::make_shared<LagrangeDiscretization>(*coarseMeshes, GetParam().coarseDegree,
GetParam().dim);
- coarse = std::make_unique<Vector>(*coarseDisc, GetParam().coarseLevel);
- fineDisc = std::make_unique<MultiPartDiscretization>(*fineMeshes, GetParam().fineDegree,
+ coarse = std::make_unique<Vector>(coarseDisc, GetParam().coarseLevel);
+ fineDisc = std::make_shared<MultiPartDiscretization>(*fineMeshes, GetParam().fineDegree,
DomainPart, GetParam().dim);
- fine = std::make_unique<Vector>(*fineDisc, GetParam().fineLevel);
+ fine = std::make_unique<Vector>(fineDisc, GetParam().fineLevel);
}
void SetUp() override {
diff --git a/test/coupled/TestCoupledConsistency.cpp b/test/coupled/TestCoupledConsistency.cpp
index 2e388f0af1e55b64fe20758f742be9866044df10..ccef0111663b9c83c60bb48eedf05af9d211d4ad 100644
--- a/test/coupled/TestCoupledConsistency.cpp
+++ b/test/coupled/TestCoupledConsistency.cpp
@@ -97,8 +97,8 @@ protected:
TEST_P(CoupledConsistencyTest, CompareCellModels) {
- Vector elphyPot(elphyA->GetDisc());
- Vector tensionPot(coupledA->GetDisc());
+ Vector elphyPot(elphyA->GetSharedDisc());
+ Vector tensionPot(coupledA->GetSharedDisc());
auto elphySolver = GetElphySolver(GetParam().ElphyModel, cellModel->ElphyModel(),
@@ -130,9 +130,9 @@ TEST_P(CoupledConsistencyTest, CompareCellModels) {
TEST_P(CoupledConsistencyTest, CompareSolvers) {
- Vector elphyPot(0.0, elphyA->GetDisc());
- Vector coupledPot(0.0, coupledA->GetDisc());
- Vector displacement(0.0, mechA->GetDisc());
+ Vector elphyPot(0.0, elphyA->GetSharedDisc());
+ Vector coupledPot(0.0, coupledA->GetSharedDisc());
+ Vector displacement(0.0, mechA->GetSharedDisc());
auto elphySolver = GetElphySolver(GetParam().ElphyModel, *cellModel,
diff --git a/test/coupled/TestDeformedDiffusion.cpp b/test/coupled/TestDeformedDiffusion.cpp
index 5f20982be8f5410add4bad1be20ad03493cec4fc..4438690203eff62073675778445e5ec8cd3544d7 100644
--- a/test/coupled/TestDeformedDiffusion.cpp
+++ b/test/coupled/TestDeformedDiffusion.cpp
@@ -78,9 +78,9 @@ protected:
};
TEST_F(DeformedDiffusionTest, CheckReference) {
- Vector elphyVec(0.0, elphyA->GetDisc());
- Vector mechVec(1.0, mechA->GetDisc());
- Vector vel(0.0, mechA->GetDisc());
+ Vector elphyVec(0.0, elphyA->GetSharedDisc());
+ Vector mechVec(1.0, mechA->GetSharedDisc());
+ Vector vel(0.0, mechA->GetSharedDisc());
std::unordered_map<Point, Tensor> conductivities{};
for (cell c = elphyVec.cells(); c != elphyVec.cells_end(); ++c) {
@@ -94,9 +94,9 @@ TEST_F(DeformedDiffusionTest, CheckReference) {
}
TEST_F(DeformedDiffusionTest, NoDeformation) {
- Vector elphyVec(0.0, elphyA->GetDisc());
- Vector mechVec(0.0, mechA->GetDisc());
- Vector vel(0.0, mechA->GetDisc());
+ Vector elphyVec(0.0, elphyA->GetSharedDisc());
+ Vector mechVec(0.0, mechA->GetSharedDisc());
+ Vector vel(0.0, mechA->GetSharedDisc());
deformedProblem->UpdateDeformation(mechVec, vel, elphyVec);
@@ -109,9 +109,9 @@ TEST_F(DeformedDiffusionTest, NoDeformation) {
}
TEST_F(DeformedDiffusionTest, WithDeformation) {
- Vector elphyVec(0.0, elphyA->GetDisc());
- Vector mechVec(0.0, mechA->GetDisc());
- Vector vel(0.0, mechA->GetDisc());
+ Vector elphyVec(0.0, elphyA->GetSharedDisc());
+ Vector mechVec(0.0, mechA->GetSharedDisc());
+ Vector vel(0.0, mechA->GetSharedDisc());
for (row r = mechVec.rows(); r != mechVec.rows_end(); ++r) {
auto x = r();
mechVec(r, 0) = 2.0 * x[0];
@@ -129,9 +129,9 @@ TEST_F(DeformedDiffusionTest, WithDeformation) {
TEST_F(DeformedDiffusionTest, StretchAlongFiber) {
- Vector elphyVec(0.0, elphyA->GetDisc());
- Vector mechVec(0.0, mechA->GetDisc());
- Vector vel(0.0, mechA->GetDisc());
+ Vector elphyVec(0.0, elphyA->GetSharedDisc());
+ Vector mechVec(0.0, mechA->GetSharedDisc());
+ Vector vel(0.0, mechA->GetSharedDisc());
for (row r = mechVec.rows(); r != mechVec.rows_end(); ++r) {
mechVec(r, 0) = 2.0 * r()[0];
}
diff --git a/test/elasticity/TestAssembleConsistency.cpp b/test/elasticity/TestAssembleConsistency.cpp
index f0d78609262bcae47d557ba42dc3110fcd949c92..5555b877dbe61587500bdf495cb1409cf2064dc8 100644
--- a/test/elasticity/TestAssembleConsistency.cpp
+++ b/test/elasticity/TestAssembleConsistency.cpp
@@ -62,11 +62,11 @@ protected:
Config::PrintInfo();
elastProblem = GetElasticityProblem();//CardiacBeam
- elastProblem->CreateMeshes();
+ elastProblem->CreateMeshes(MeshesCreator());
assemble = CreateFiniteElasticityAssemble(*elastProblem, model,
GetParam().degree, false);
- u = std::make_unique<Vector>(0.0, assemble->GetDisc());
+ u = std::make_unique<Vector>(0.0, assemble->GetSharedDisc());
assemble->SetInitialValue(*u);
assemble->Initialize(*u);