Made AHA evaluation possible

cardmech/conf/elasticity.conf

@@ -60,7 +60,7 @@ PrestressSteps=5
 #       Meshes       #
 ######################
 #Mesh = BenchmarkBeam2DTet
-#Mesh = FullEllipsoid
+Mesh = FullEllipsoid
 Mesh = QuarterEllipsoid
 #Mesh = OrientedEllipsoid
 #Mesh = KovachevaHeart

cardmech/src/elasticity/MainElasticity.cpp

@@ -181,6 +181,14 @@ std::vector<double> MainElasticity::Evaluate(const Vector &solution) {
                    PrintInfoEntry("H1 Error", mechA->H1Error(solution), 1),
                    PrintInfoEntry("Energy Error", mechA->EnergyError(solution), 1));
     mout << "=============================================================" << endl;
+
+    auto ahaStress = mechA->AHAStress(solution, Point(0.0, 0.0, -20.0));
+    for (int i = 0; i < ahaStress.size(); ++i) {
+      mout << "AHA Strain Segment " << i + 1 << ": " << ahaStress[i][0] << endl;
+    }
+    SerialVtuPlot lrocPlot(solution);
+    mechA->PlotRLoC(solution, lrocPlot, Point(0.0, 0.0, -20));
+    lrocPlot.PlotFile("LRoC");
   }
   return elasticityProblem->EvaluationResults(solution);
 

cardmech/src/elasticity/assemble/FiniteElasticity.hpp

@@ -78,7 +78,9 @@ public:
   virtual void GetInvariant(const Vector &u, Vector &iota4) const = 0;
 
 
-  virtual void PlotFibres(const Vector &u, VtuPlot &plot) const {  }
+  virtual void PlotFibres(const Vector &u, VtuPlot &plot) const {}
+
+  virtual void PlotRLoC(const Vector &u, VtuPlot &plot, const Point &apex) const {}
 
   void PlotDisplacement(const Vector &u, int step, const string &varname) const;
 
@@ -101,6 +103,9 @@ public:
 
   virtual double MaxStress(const Vector &u) const = 0;
 
+  virtual std::array<VectorField, 17>
+  AHAStress(const Vector &u, const Point &apex) const { return std::array<VectorField, 17>{}; }
+
   void SetAcceleration(bool toggle) {
     rho = (toggle ? base_rho : 0.0);
   }

cardmech/src/elasticity/assemble/LagrangeElasticity.cpp

@@ -614,9 +614,9 @@ void LagrangeElasticity::GetInvariant(const Vector &u, Vector &iota4) const {
 
     const VectorField fibre = Direction(c.GetData());
     for (int j = 0; j < E.size(); ++j) {
-      Tensor Du = E.VectorGradient(E[j](), u);
-      iota4(E[j](), 0) += mat::invariant(4, Du, fibre);
-      iota4(E[j](), 1) += 1.0;
+      Point xi = E[j]();
+      iota4(xi, 0) += mat::invariant(4, DeformationGradient(E, xi, u), fibre);
+      iota4(xi, 1) += 1.0;
     }
   }
   //iota4.Collect();
@@ -654,6 +654,61 @@ void LagrangeElasticity::PlotFibres(const Vector &u, VtuPlot &plot) const {
   plot.AddData("Normal", n);
 }
 
+
+void LagrangeElasticity::PlotRLoC(const Vector &u, VtuPlot &plot, const Point &apex) const {
+  LagrangeDiscretization fibreDisc(u.GetDisc().GetMeshes(), 0, SpaceDimension);
+  Vector l(0.0, fibreDisc);
+  Vector r(0.0, fibreDisc);
+  Vector t(0.0, fibreDisc);
+  Vector lroc(0.0, fibreDisc);
+
+  for (cell c = u.cells(); c != u.cells_end(); ++c) {
+    if (getCellChamber(c.Subdomain()) != 0) continue;
+    int cellSegment = getCellMaterial(c.Subdomain());
+    if (cellSegment > 17) continue;
+
+    VectorFieldElement E(u, c);
+    auto F = DeformationGradient(E, c(), u);
+    auto f = F * Direction(c.GetData());
+    auto s = F * Sheet(c.GetData());
+    auto n = F * Normal(c.GetData());
+
+    // Calculate R-Lo-C Coordinate system
+    // https://doi.org/10.1053/euje.2000.0031
+    VectorField dir = apex - c();
+    dir /= norm(dir);
+
+    VectorField rad = Sheet(c.GetData());
+    rad *= (rad * dir > 0 ? -1.0 : 1.0);
+    VectorField lon = -1.0 * (dir - (rad * dir) / (rad * rad) * rad);
+    lon /= norm(lon);
+    VectorField cir = rad ^ lon;
+    cir /= norm(cir);
+
+    Tensor Q(rad, lon, cir);
+
+    double strainf = sqrt(f * f) - 1.0;
+    double strains = sqrt(s * s) - 1.0;
+    double strainn = sqrt(n * n) - 1.0;
+
+    VectorField strain = (Q * (f / norm(f))) * strainf
+                         + (Q * (s / norm(s))) * strains
+                         + (Q * (n / norm(n))) * strainn;
+
+    for (int d = 0; d < SpaceDimension; ++d) {
+      r(c(), d) = rad[d];
+      l(c(), d) = lon[d];
+      t(c(), d) = cir[d];
+      lroc(c(), d) = strain[d];
+    }
+  }
+
+  plot.AddData("Longitdinal", l);
+  plot.AddData("Radial", r);
+  plot.AddData("Circumferential", t);
+  plot.AddData("LRoC-Strain", lroc);
+}
+
 double LagrangeElasticity::InitialVolume(const Mesh &M) const {
   double volume = 0.0;
   for (cell c = M.cells(); c != M.cells_end(); ++c) {
@@ -679,3 +734,52 @@ double LagrangeElasticity::DeformedVolume(const Vector &U) const {
   }
   return PPM->Sum(volume);
 }
+
+std::array<VectorField, 17>
+LagrangeElasticity::AHAStress(const Vector &u, const Point &apex) const {
+  std::array<VectorField, 17> ahaStress{};
+  std::array<int, 17> ahaCells{};
+
+  for (cell c = u.cells(); c != u.cells_end(); ++c) {
+    if (getCellChamber(c.Subdomain()) != 0) continue;
+    int cellSegment = getCellMaterial(c.Subdomain());
+    if (cellSegment > 17) continue;
+
+    VectorFieldElement E(u, c);
+    auto F = DeformationGradient(E, c(), u);
+    auto f = F * Direction(c.GetData());
+    auto s = F * Sheet(c.GetData());
+    auto n = F * Normal(c.GetData());
+
+    // Calculate R-Lo-C Coordinate system
+    // https://doi.org/10.1053/euje.2000.0031
+    VectorField dir = apex - c();
+    dir /= norm(dir);
+
+    VectorField r = Sheet(c.GetData());
+    r *= (r * dir > 0 ? -1.0 : 1.0);
+    VectorField l = -1.0 * (dir - (r * dir) / (r * r) * r);
+    l /= norm(l);
+    VectorField t = r ^ l;
+    t /= norm(t);
+
+
+    Tensor Q(r, l, t);
+
+    double strainf = sqrt(f * f) - 1.0;
+    double strains = sqrt(s * s) - 1.0;
+    double strainn = sqrt(n * n) - 1.0;
+
+    VectorField strain = (Q * (f / norm(f))) * strainf
+                         + (Q * (s / norm(s))) * strains
+                         + (Q * (n / norm(n))) * strainn;
+    ahaStress[cellSegment - 1] += strain;
+    ahaCells[cellSegment - 1] += 1;
+  }
+
+  for (int i = 0; i < ahaStress.size(); ++i) {
+    ahaStress[i] = PPM->Sum(ahaStress[i]) / PPM->Sum(ahaCells[i]);
+  }
+  return ahaStress;
+}
+

cardmech/src/elasticity/assemble/LagrangeElasticity.hpp

@@ -64,9 +64,12 @@ public:
 
   void PlotFibres(const Vector &u, VtuPlot &plot) const override;
 
+  void PlotRLoC(const Vector &u, VtuPlot &plot, const Point &apex) const override;
+
 
   double MaxStress(const Vector &u) const override;
 
+  std::array<VectorField, 17> AHAStress(const Vector &u, const Point &apex) const override;
 
   double NormPrestress();