removed IProblem, added CardMeshIProblem

mpp @ dac6a503

-Subproject commit ce3fae26092f94695226fe1051cd05e4cec86bd8
+Subproject commit dac6a50392fec76f0f325c36a446d928f871b77a

src/CardMechIProblem.hpp

+#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

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);

src/elasticity/problem/ElasticityProblem.hpp

 #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

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; }