Merge branch '152-add-memory-efficient-projection' into 'feature'

Resolve "Add memory efficient projection"

Closes #152

See merge request !166

cardmech/conf/elasticity_benchmarks/beam.conf

@@ -9,7 +9,7 @@ h_min=1e-6
 
 Model = ActiveStrainElasticity
 MechModel = Static
-MechDiscretization = DG
+#MechDiscretization = DG
 MechDiscretization = Conforming
 
 CalculateReference = 1

cardmech/conf/m++.conf

@@ -11,14 +11,14 @@ GeoPath = ../geo/
 # = Choose Configuration for Problems  =
 # ======================================
 # === Electrophysiology ================
-loadconf = elphy.conf
+#loadconf = elphy.conf
 #loadconf = monodomain.conf;
 # === Cardiac Mechanics ================
 #loadconf = coupled.conf
 #loadconf = klotz.conf
 #loadconf = elasticity_benchmarks/test.conf
 #loadconf = elasticity_benchmarks/linear.conf
-#loadconf = elasticity_benchmarks/beam.conf
+loadconf = elasticity_benchmarks/beam.conf
 #loadconf = elasticity_benchmarks/ellipsoid.conf
 #loadconf = elasticity.conf
 # ======================================

cardmech/src/elasticity/MainElasticity.cpp

 #include <DGDiscretization.hpp>
 #include <EGDiscretization.hpp>
 #include "VectorFieldElement.hpp"
-#include "VectorFieldFaceElement.hpp"
 
 #include "assemble/LagrangeElasticity.hpp"
 #include "assemble/DGElasticity.hpp"
@@ -117,7 +116,9 @@ void MainElasticity::initAssemble() {
     int degree = 1;
     bool isStatic = false;
     config.get("MechPolynomialDegree",  degree);
-    mechA = std::make_unique<LagrangeElasticity>(*elasticityProblem, *pressureProblem, *meshes, degree, isStatic);
+    mechA = CreateFiniteElasticityAssemble(*elasticityProblem, *pressureProblem, *meshes, degree, isStatic);
+
+
 }
 
 double MainElasticity::EvaluateQuantity(const Vector &solution, const std::string &quantity) const {

cardmech/src/elasticity/assemble/DGElasticity.cpp

@@ -2,7 +2,6 @@
 #include "problem/PressureProblem.hpp"
 #include "DGDoF.hpp"
 #include "DGVectorFieldElement.hpp"
-#include "DGVectorFieldFaceElement.hpp"
 #include "VectorAccess.hpp"
 
 DGElasticity::DGElasticity(const ElasticityProblem &eP, const PressureProblem &pP, const Meshes &mesh, int degree, bool isStatic)
@@ -621,10 +620,8 @@ void DGElasticity::GetStress(const Vector &displacement, Vector &stress) const {
 
 
 
-Vector DGElasticity::Projection(const Vector &u, const Vector &reference) const {
-  Vector projection(0.0, u);
-
-  int dlevel = reference.SpaceLevel() - u.SpaceLevel();
+void DGElasticity:: Project( const Vector &fine, Vector &u) const  {
+  int dlevel = fine.SpaceLevel() - u.SpaceLevel();
 
   std::vector<Point> childrenVerteces{};
   std::vector<Point> childVerteces{};
@@ -642,20 +639,20 @@ Vector DGElasticity::Projection(const Vector &u, const Vector &reference) const
         }
       }
     }
-    auto r = projection.find_row(c());
+    auto r = u.find_row(c());
     int childCount = recursiveCells[dlevel].size();
-    DGVectorFieldElement E(projection, c);
+    DGVectorFieldElement E(u, c);
 
     for(int i = 0;i< E.NodalPoints(); ++i){
       Point z = E.NodalPoint(i);
       bool found=false;
       for (const auto &child: recursiveCells[dlevel]) {
-        DGVectorFieldElement CE(reference, reference.find_cell(child->Center()));
+        DGVectorFieldElement CE(fine, fine.find_cell(child->Center()));
         for(int j = 0;j< CE.NodalPoints(); ++j){
           if(CE.NodalPoint(j) == z){
-            auto rf = reference.find_row(child->Center());
+            auto rf = fine.find_row(child->Center());
             for(int k=0;k<E.Dim();++k) {
-              projection(r, k*E.shape_size() + i) = reference(rf, k*CE.shape_size() + j);
+              u(r, k*E.shape_size() + i) = fine(rf, k*CE.shape_size() + j);
             }
             found = true;
           }
@@ -667,19 +664,16 @@ Vector DGElasticity::Projection(const Vector &u, const Vector &reference) const
       }
     }
   }
-  return projection;
 }
 
 
-Vector DGElasticity::Interpolation(const Vector &start, const Vector &u) const {
-  Vector interpolation(0.0, u);
-
-  int dlevel = u.SpaceLevel() - start.SpaceLevel();
+void DGElasticity::Interpolate(const Vector &coarse, Vector &u) const {
+  int dlevel = u.SpaceLevel() - coarse.SpaceLevel();
 
   std::vector<Point> childrenVerteces{};
   std::vector<Point> childVerteces{};
   std::vector<Point> coarseNodalPoints{};
-  for (cell c = start.cells(); c != start.cells_end(); ++c) {
+  for (cell c = coarse.cells(); c != coarse.cells_end(); ++c) {
     std::vector<std::vector<std::unique_ptr<Cell>>> recursiveCells(dlevel + 1);
     recursiveCells[0].emplace_back(CreateCell(c.Type(), c.Subdomain(), c.AsVector(), false));
     for (int l = 1; l <= dlevel; ++l) {
@@ -694,22 +688,22 @@ Vector DGElasticity::Interpolation(const Vector &start, const Vector &u) const {
     }
 
     int childCount = recursiveCells[dlevel].size();
-    DGVectorFieldElement E(start, c);
+    DGVectorFieldElement E(coarse, c);
 
 
     for (const auto &child: recursiveCells[dlevel]) {
-      auto fineRow = interpolation.find_row(child->Center());
+      auto fineRow = u.find_row(child->Center());
       //mout << "E.shape_size() " << E.shape_size() << endl;
-      cell fineCell = interpolation.find_cell(child->Center());
-      DGVectorFieldElement CE(start, fineCell);
+      cell fineCell = u.find_cell(child->Center());
+      DGVectorFieldElement CE(coarse, fineCell);
 
       // TODO: test for (int i = 0; i < E.shape_size(); ++i)
       for (int i = 0; i < CE.shape_size(); ++i) {
           auto x = c.GlobalToLocal(CE.NodalPoint(i));
-          VectorField val = E.VectorValue(x, start);
+          VectorField val = E.VectorValue(x, coarse);
         //mout << "u ( " << c.GlobalToLocal(child->Corner(i)) << " ) = " << val << endl;
         for (int k = 0; k < E.Dim(); ++k) {
-          interpolation(fineRow, k * E.shape_size() + i) = val[k];
+          u(fineRow, k * E.shape_size() + i) = val[k];
         }
       }
     }
@@ -724,8 +718,7 @@ Vector DGElasticity::Interpolation(const Vector &start, const Vector &u) const {
     }
   }*/
 
-  interpolation.Accumulate();
-  return interpolation;
+  u.Accumulate();
 }
 
 
@@ -767,8 +760,7 @@ double DGElasticity::EnergyNorm(const Vector &u) const {
 }
 
 double DGElasticity::EnergyError(const Vector &u, const Vector &reference) const {
-    Vector projection(0.0,u);
-    projection = Projection(u, reference);
+    Vector projection = Projection(u, reference);
     projection -= u;
     return EnergyNorm(projection);
 }

cardmech/src/elasticity/assemble/DGElasticity.hpp

@@ -60,9 +60,10 @@ public:
 
     double PressureEnergy(const cell &c, int face, int bc, const Vector &U) const;
 
-  Vector Interpolation(const Vector &start, const Vector &u) const override;
 
-    Vector Projection(const Vector &u, const Vector &reference) const override;
+  void Project( const Vector &fine, Vector &u) const override;
+
+  void Interpolate(const Vector &coarse, Vector &u) const override;
 
     //void Project(Vector &coarse, const Vector &fine) const;
 

cardmech/src/elasticity/assemble/EGElasticity.cpp

@@ -11,7 +11,6 @@
 #include "EGDoF.hpp"
 #include "EGDiscretization.hpp"
 #include "EGVectorFieldElement.hpp"
-#include "EGVectorFieldFaceElement.hpp"
 
 EGElasticity::EGElasticity(const ElasticityProblem &eP, const PressureProblem &pP, const Meshes &mesh, int degree, bool isStatic)
     : FiniteElasticity(eP, pP) , disc(mesh, degree, size){

cardmech/src/elasticity/assemble/FiniteElasticity.cpp

-#include <VectorFieldFaceElement.hpp>
 //#include <LagrangeDiscretization.hpp>
 #include "Plotting.hpp"
 #include "FiniteElasticity.hpp"

cardmech/src/elasticity/assemble/FiniteElasticity.hpp

@@ -86,17 +86,28 @@ public:
 
   void PlotDisplacement(const Vector &u, int step, const string &varname) const;
 
-  virtual Vector Projection(const Vector &u, const Vector &reference) const {
+  Vector Projection(const Vector &target, const Vector &fine) const{
+    Vector v(target);
+    Project(fine, v);
+    return v;
+  }
+
+  virtual void Project(const Vector &fine, Vector &u) const {
     Warning("No Projection implemented")
-    return Vector(u);
   }
 
-  virtual Vector Interpolation(const Vector &start, const Vector &u) const {
+  Vector Interpolation(const Vector &coarse, const Vector &target) const{
+    Vector v(target);
+    Interpolate(coarse, v);
+    return v;
+  }
+
+  virtual void Interpolate(const Vector &coarse, Vector &u) const {
     Warning("No Interpolation implemented")
-    return Vector(u);
   }
 
 
+
   virtual double L2Error(const Vector &u) const { return 0.0; };
   virtual double L2Error(const Vector &u, const Vector &reference) const { return 0.0; };
 

cardmech/src/elasticity/assemble/LagrangeElasticity.cpp

@@ -785,30 +785,20 @@ double LagrangeElasticity::DeformedVolume(const Vector &U) const {
   return PPM->Sum(volume);
 }
 
-Vector LagrangeElasticity::Projection(const Vector &u, const Vector &reference) const {
-  Vector projection(0.0, u);
-  /*mout << "u " << u << endl;
-  mout << projection << endl;
-  mout << "reference before Transfer " << reference << endl;*/
-  LagrangeTransfer T(projection, reference);
-  //mout << "reference " << reference << endl;
-  T.Project(projection, reference);
-  /*mout << "new projection " << projection;
-  mout << L2Norm(projection-u)<<endl;*/
-  return projection;
+void LagrangeElasticity::Project(const Vector &fine, Vector &u) const {
+  LagrangeTransfer T(u, fine);
+  T.Project(u, fine);
 }
 
-Vector LagrangeElasticity::Interpolation(const Vector &start, const Vector &u) const {
-  Vector interpolation(0.0, u);
-  LagrangeTransfer T(start, u);
-  T.Prolongate(start, interpolation);
+void LagrangeElasticity::Interpolate(const Vector &coarse, Vector &u) const {
+  LagrangeTransfer T(coarse, u);
+  T.Prolongate(coarse, u);
 
   mout.PrintInfo("Interpolation comparison", 10,
-                 PrintInfoEntry("Coarse L2 Norm", L2Norm(start), 1),
-                 PrintInfoEntry("Interpolated L2 Norm", L2Norm(interpolation), 1),
-                 PrintInfoEntry("Coarse L2 Avg Norm", L2AvgNorm(start), 1),
-                 PrintInfoEntry("Interpolated L2 Avg Norm", L2AvgNorm(interpolation), 1));
-  return interpolation;
+                 PrintInfoEntry("Coarse L2 Norm", L2Norm(coarse), 1),
+                 PrintInfoEntry("Interpolated L2 Norm", L2Norm(u), 1),
+                 PrintInfoEntry("Coarse L2 Avg Norm", L2AvgNorm(coarse), 1),
+                 PrintInfoEntry("Interpolated L2 Avg Norm", L2AvgNorm(u), 1));
 }
 
 void LagrangeElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {

cardmech/src/elasticity/assemble/LagrangeElasticity.hpp

@@ -59,8 +59,6 @@ public:
 
   double PressureEnergy(const cell &c, int face, int bc, const Vector &U) const;
 
-  Vector Projection(const Vector &u, const Vector &reference) const override;
-
   double L2Error(const Vector &u, const Vector &reference) const override;;
 
   double H1Error(const Vector &u, const Vector &reference) const override;
@@ -105,7 +103,9 @@ public:
 
   double PressureBndNorm(const Vector &u) const override;
 
-  Vector Interpolation(const Vector &start, const Vector &u) const override;
+  void Project( const Vector &fine, Vector &u) const override;
+
+  void Interpolate(const Vector &coarse, Vector &u) const override;
 
   void SetDisplacement(Vector &u) const override;
 

cardmech/src/elasticity/solvers/IterativePressureSolver.cpp

@@ -76,7 +76,6 @@ void IterativePressureSolver::Solve(Vector &displacement, int pressureSteps) {
 bool IterativePressureSolver::SolveStep(Newton &newton, double t, double dt, Vector &uNew) {
   std::array<double, 4> mechVolumes = mechA.ChamberVolume(uNew);
   problem.Initialize(dt, t, mechVolumes);
-  newton.PrintIterationInfo();
   for (int i = 0; i < maxPressureIterations; ++i) {
     mout << "Update Pressure" << endl;
     problem.Update(dt, t, mechVolumes);
@@ -88,6 +87,7 @@ bool IterativePressureSolver::SolveStep(Newton &newton, double t, double dt, Vec
                    PrintInfoEntry("LA Pressure", problem.Pressure(t, Origin, PRESSURE_BC_LA), 3),
                    PrintInfoEntry("RA Pressure", problem.Pressure(t, Origin, PRESSURE_BC_RA), 3));
     newton(mechA, uNew);
+    newton.PrintIterationInfo();
 
     mechVolumes = mechA.ChamberVolume(uNew);
     if (problem.Converged(mechVolumes)) {