Added optimized interpolation and projection

cardmech/src/elasticity/assemble/DGElasticity.cpp

@@ -621,10 +621,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 +640,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 +665,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 +689,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 +719,7 @@ Vector DGElasticity::Interpolation(const Vector &start, const Vector &u) const {
     }
   }*/
 
-  interpolation.Accumulate();
-  return interpolation;
+  u.Accumulate();
 }
 
 
@@ -767,8 +761,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/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;