diff --git a/cardmech/conf/elasticity.conf b/cardmech/conf/elasticity.conf
index f85a8e43b0c711761072ce33a40979aade9e4d9d..cdec49e5eb35afe2e86f4e60c725e5cc2daeff9d 100644
--- a/cardmech/conf/elasticity.conf
+++ b/cardmech/conf/elasticity.conf
@@ -60,7 +60,7 @@ PrestressSteps=5
# Meshes #
######################
#Mesh = BenchmarkBeam2DTet
-#Mesh = FullEllipsoid
+Mesh = FullEllipsoid
Mesh = QuarterEllipsoid
#Mesh = OrientedEllipsoid
#Mesh = KovachevaHeart
diff --git a/cardmech/src/elasticity/MainElasticity.cpp b/cardmech/src/elasticity/MainElasticity.cpp
index 2f990806e0bfc7c0a972c54e30786e476b510e90..4342774df70bacad36f2abcceab21f5370bd5ca5 100644
--- a/cardmech/src/elasticity/MainElasticity.cpp
+++ b/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);
diff --git a/cardmech/src/elasticity/assemble/FiniteElasticity.hpp b/cardmech/src/elasticity/assemble/FiniteElasticity.hpp
index 6ab25c74945e65721425e70229cb7f169e46d160..8d23fc10c5f015980452558baf66c07ce12ff3e4 100644
--- a/cardmech/src/elasticity/assemble/FiniteElasticity.hpp
+++ b/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);
}
diff --git a/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp b/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp
index 7ac17bb81c1c3a75168de8f995aadb65e9d1e4b3..875c98649f497e867c4907b9a784ab06ba64529f 100644
--- a/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp
+++ b/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;
+}
+
diff --git a/cardmech/src/elasticity/assemble/LagrangeElasticity.hpp b/cardmech/src/elasticity/assemble/LagrangeElasticity.hpp
index 2a223b83e767b36df66820f439a045ff5512b31f..09f3ca9f0e07c6cf2df34e00ddddaac25b50682f 100644
--- a/cardmech/src/elasticity/assemble/LagrangeElasticity.hpp
+++ b/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();