Streamlined conductivity storage

src/coupled/problem/CoupledProblem.hpp

@@ -204,8 +204,8 @@ public:
     elphyProblem->InitializeEvaluationPoints();
   }
 
-  const Tensor GetConductivity(int domainpart, VectorField &fiber) const override {
-    return elphyProblem->GetConductivity(domainpart, fiber);
+  const Tensor Conductivity(int domainpart, VectorField &fiber) const override {
+    return elphyProblem->Conductivity(domainpart, fiber);
   }
 
   vector<double> EvaluationResults(const Vector &solution) const override {

src/electrophysiology/assemble/IElphyAssemble.cpp

@@ -224,7 +224,7 @@ void IElphyAssemble::PrintMesh(const Vector &u) const {
 
 Tensor IElphyAssemble::referenceConductivity(const cell &c) const {
   auto f = Direction(c.GetData());
-  const Tensor &Sigma = problem.GetConductivity(DomainPart(*c), f);
+  const Tensor &Sigma = problem.Conductivity(DomainPart(*c), f);
   // Deformed directions are not orthogonal
   return Sigma;
 }

src/electrophysiology/problem/DiffusionProblem.cpp

@@ -13,9 +13,6 @@ DiffusionProblem::DiffusionProblem(std::array<double, 9> entries) :
   if (!meshFromConfig()) {
     meshesName = "UnitBlock" + std::to_string(spaceDim) + "D" + (isTet ? "Tet" : "Quad");
   }
-
-  SigmaAtria = conductivity;
-  SigmaVentricle = conductivity;
 }
 
 
@@ -32,7 +29,8 @@ vector<double> DiffusionProblem::EvaluationResults(const Vector &solution) const
 
   return eval;
 }
-const Tensor DiffusionProblem::GetConductivity(int domainpart,VectorField &fiber) const{
+
+const Tensor DiffusionProblem::Conductivity(int domainpart, VectorField &fiber) const {
   return conductivity;
 }
 

src/electrophysiology/problem/DiffusionProblem.hpp

@@ -12,11 +12,13 @@ protected:
 public:
   explicit DiffusionProblem(std::array<double, 9> entries);
 
-  std::string Name() const  override {
+  std::string Name() const override {
     return "Diffusion";
   }
+
   vector<double> EvaluationResults(const Vector &solution) const override;
-  const Tensor GetConductivity(int domainpart,VectorField &fiber) const override;
+
+  const Tensor Conductivity(int domainpart, VectorField &fiber) const override;
 };
 
 class DiffusionExampleOne : public DiffusionProblem {

src/electrophysiology/problem/ElphyBiVentricleProblem.hpp

@@ -13,16 +13,14 @@ public:
     std::string spaceSmoothing;
     Config::Get("ExternalCurrentSmoothingInSpace", spaceSmoothing);
 
+    diffusionValues[0] = 0.001 * VectorField(0.2073, 0.0921, 0.023);
+    diffusionValues[1] = diffusionValues[0];
+
     if (spaceSmoothing == "discrete") {
       meshesName = "BiVentricle_unsmoothed";
     } else {
       meshesName = "BiVentricle_smoothed";
     }
-    Tensor SV(0.2073, 0.0, 0.0,
-              0.0, 0.0921, 0.0,
-              0.0, 0.0, 0.023);
-    //Divide for correct Units from 1/m to 1/mm
-    SigmaVentricle = (1.0 / 1000.0) * SV;
   }
 
   std::string Name() const override {

src/electrophysiology/problem/ElphyChamberCube2dProblem.hpp

@@ -8,8 +8,7 @@
 class ElphyChamberCube2dProblem : public ElphyProblem {
 public:
   explicit ElphyChamberCube2dProblem() : ElphyProblem(false) {
-    Tensor SV=SigmaVentricle;
-    SigmaAtria=SV;
+    diffusionValues[1] = diffusionValues[0];
 
     if (!meshFromConfig()) {
       meshesName = "ChamberCube0252DQuad";

src/electrophysiology/problem/ElphyFullHeartProblem.hpp

@@ -8,11 +8,6 @@
 class ElphyFullHeartProblem : public ElphyProblem {
 public:
   explicit ElphyFullHeartProblem() : ElphyProblem(false) {
-    Tensor SV(0.2073, 0.0, 0.0,
-              0.0, 0.0921, 0.0,
-              0.0, 0.0, 0.023);
-    //Divide for correct Units from 1/m to 1/mm
-    SigmaVentricle = (1.0 / 1000.0) * SV;
     if (!meshFromConfig()) {
       meshesName = "HeartT4";
     }
@@ -21,7 +16,6 @@ public:
   std::string Name() const override{
     return "FullHeart";
   }
-  //std::string MeshName() const override { return meshFromConfig()? meshName : "HeartT4"; }
 
   void InitializeEvaluationPoints() override;
 

src/electrophysiology/problem/ElphyProblem.hpp

@@ -48,9 +48,9 @@ protected:
   std::vector<Point> evaluationPoints{};
   int numberEvaluationPointsBesideMesh = 0;
   std::vector<EvaluationPointData> cellsBesideMeshPoints{};
-  std::array<double, 4> sigma;
-  Tensor SigmaVentricle;
-  Tensor SigmaAtria;
+  std::array<double, 4> sigma{};
+  /// Stores the diagonal entries of the diffusion tensor for the ventricles (0) and atria (1).
+  std::array<VectorField, 2> diffusionValues{};
 public:
   bool updateActivationTime{true};
 
@@ -60,67 +60,54 @@ public:
 
   IElphyProblem(bool at, std::array<double, 4> sigmaParams)
       : sigma(sigmaParams), updateActivationTime(at) {
-    InitializeConductivity();
-  }
-
-  void
-  SetNumberOfEvaluationPointsBesideMesh(int number) { numberEvaluationPointsBesideMesh = number; };
-
-  const int
-  getNumberEvaluationPointsBesideMesh() const { return numberEvaluationPointsBesideMesh; };
 
-  const int getSizeCellsEvaluationPoints() const { return cellsBesideMeshPoints.size(); };
-
-  void SetEvaluationCells(const Vector &V);
+    double sigma_l = (sigma[0] * sigma[2]) / (sigma[0] + sigma[2]);
+    double sigma_t = (sigma[1] * sigma[3]) / (sigma[1] + sigma[3]);
+    //Divide for correct Units from 1/m to 1/mm
+    diffusionValues[0] = (1.0 / 1000.0) * VectorField(sigma_l, sigma_t, sigma_t);
+    diffusionValues[1] = (2.0 / 1000.0) * VectorField(2.2188 * 0.251, 0.251,
+                                                      0.251); //TODO scaling mit 2 automatisieren
+  }
 
-  const EvaluationPointData &getCellAt(int i) const { return cellsBesideMeshPoints.at(i); };
 
   bool Includes(const cell &c, const Point P);
 
-  void InitializeConductivity() {
-    double sigma_l = (sigma[0] * sigma[2]) / (sigma[0] + sigma[2]);
-    double sigma_t = (sigma[1] * sigma[3]) / (sigma[1] + sigma[3]);
-    Tensor SV(sigma_l, 0.0, 0.0,
-              0.0, sigma_t, 0.0,
-              0.0, 0.0, sigma_t);
-    //Divide for correct Units from 1/m to 1/mm
-    SigmaVentricle = (1.0 / 1000.0) * SV;
-    Tensor SA(2.2188 * 0.251, 0, 0,
-              0, 0.251, 0,
-              0, 0, 0.251);
-    SigmaAtria = 2 * (1.0 / 1000.0) * SA; //TODO scaling mit 2 automatisiern
+  virtual double Potential(double time, const Point &x) const { return 0.0; }
 
-  }
+  virtual double IonCurrent(double time, const Point &x) const { return 0.0; }
 
-  virtual const Tensor GetConductivity(int domainpart, VectorField &fiber) const {
+  virtual const Tensor Conductivity(int domainpart, VectorField &fiber) const {
     auto S = Product(fiber, fiber);
-    if (domainpart == 1) {
-
-      return (SigmaAtria[0][0] * S + SigmaAtria[1][1] * (One - S));//SigmaAtria;
-    } else {
-      return (SigmaVentricle[0][0] * S + SigmaVentricle[1][1] * (One - S));//SigmaVentricle;
-    }
+    return diffusionValues[domainpart][1] * One
+           + (diffusionValues[domainpart][0] - diffusionValues[domainpart][1]) * S;
   }
 
   std::vector<double>
   getEvaluationAt(const Vector &solution, const std::vector<Point> points) const;
 
-  void setUpdateActivationTime(bool b) { updateActivationTime = b; }
-
   const int getSizeEvaluationPoints() const { return evaluationPoints.size(); };
 
   const Point &getEvaluationPointAt(int i) const { return evaluationPoints.at(i); }
 
-  virtual double Potential(double time, const Point &x) const { return 0.0; }
+  virtual void chamberSpecialization(const string &chamberName) const {}
 
-  virtual double IonCurrent(double time, const Point &x) const { return 0.0; }
+  void
+  SetNumberOfEvaluationPointsBesideMesh(int number) { numberEvaluationPointsBesideMesh = number; };
 
-  virtual void chamberSpecialization(const string &chamberName) const {}
+  const int
+  getNumberEvaluationPointsBesideMesh() const { return numberEvaluationPointsBesideMesh; };
+
+  const int getSizeCellsEvaluationPoints() const { return cellsBesideMeshPoints.size(); };
+
+  void SetEvaluationCells(const Vector &V);
+
+  const EvaluationPointData &getCellAt(int i) const { return cellsBesideMeshPoints.at(i); };
 };
 
 class ElphyProblem : public IElphyProblem, public IProblem {
 protected:
   bool meshFromConfig() const { return meshesName != ""; }
+
 public:
   ElphyProblem() : IProblem(), IElphyProblem() {}
 

src/electrophysiology/problem/VariableMeshProblem.hpp

@@ -8,9 +8,7 @@
 class VariableMeshProblem : public ElphyProblem {
 public:
   VariableMeshProblem() : ElphyProblem(false) {
-    Tensor SV=SigmaVentricle;
-    //Divide for correct Units from 1/m to 1/mm
-    SigmaAtria=SV;
+    diffusionValues[1] = diffusionValues[0];
 
 
     if (!meshFromConfig()) {
@@ -27,8 +25,6 @@ public:
   std::string Evaluate(const Vector &solution) const override;
 
   std::vector<double> EvaluationResults(const Vector &solution) const override;
-
-
 };