diff --git a/mpp b/mpp
index ce3fae26092f94695226fe1051cd05e4cec86bd8..dac6a50392fec76f0f325c36a446d928f871b77a 160000
--- a/mpp
+++ b/mpp
@@ -1 +1 @@
-Subproject commit ce3fae26092f94695226fe1051cd05e4cec86bd8
+Subproject commit dac6a50392fec76f0f325c36a446d928f871b77a
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/elasticity/problem/ElasticityProblem.cpp b/src/elasticity/problem/ElasticityProblem.cpp
index 07f9d5c93e99299bd7140a4938e7d2b571c8a9b7..9be4e99f3be0888b2dcfb51f1eb90fa351875533 100644
--- a/src/elasticity/problem/ElasticityProblem.cpp
+++ b/src/elasticity/problem/ElasticityProblem.cpp
@@ -81,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"};
@@ -95,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 d4e7a1bab07bb9b7e6a688d2978bd135b21fc0f8..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,75 +41,109 @@ 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; }
+ 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"};
@@ -118,55 +153,64 @@ protected:
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; }
-
+ 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;
- }
+ 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()) {}
@@ -177,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);
}
@@ -203,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);
}
@@ -231,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";
@@ -243,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/electrophysiology/problem/ElphyProblem.hpp b/src/electrophysiology/problem/ElphyProblem.hpp
index 25aa156af5d328ecec665e31d6dd8cf8b363ecb6..d8c52d179568fef1afb3bd1f681a4a683656c73d 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"
@@ -123,17 +123,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; }