Skip to content
Snippets Groups Projects
Commit 97320b48 authored by Patrick Jaap's avatar Patrick Jaap
Browse files

Fix in mooneyrivlindensity.hh

Square the determinant in the Ciarlet variant
parent 603d9327
No related branches found
No related tags found
1 merge request!67Fix in mooneyrivlindensity
Pipeline #47848 passed
...@@ -66,7 +66,7 @@ public: ...@@ -66,7 +66,7 @@ public:
C = F^TF and B = FF^T have the same eigenvalues - we can use either one of them. C = F^TF and B = FF^T have the same eigenvalues - we can use either one of them.
There are three Mooney-Rivlin-Variants: There are three Mooney-Rivlin-Variants:
ciarlet: W(F) = mooneyrivlin_a*(normF)^2 + mooneyrivlin_b*(normFinv)^2*detF + mooneyrivlin_c*(detF)^2 - ciarlet: W(F) = mooneyrivlin_a*(normF)^2 + mooneyrivlin_b*(normFinv)^2*detF^2 + mooneyrivlin_c*(detF)^2 -
((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF) ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF)
log: W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * ln(det(F))^2 log: W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * ln(det(F))^2
square: W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * 0.5 * (det(F) - 1)^2 square: W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * 0.5 * (det(F) - 1)^2
...@@ -99,7 +99,7 @@ public: ...@@ -99,7 +99,7 @@ public:
gradientInverse.invert(); gradientInverse.invert();
field_type frobeinusNormFInverseSquared = gradientInverse.frobenius_norm2(); field_type frobeinusNormFInverseSquared = gradientInverse.frobenius_norm2();
using std::log; using std::log;
return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF); return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF);
} }
else else
{ {
......
...@@ -167,10 +167,10 @@ int main (int argc, char *argv[]) ...@@ -167,10 +167,10 @@ int main (int argc, char *argv[])
parameters["mooneyrivlin_a"] = "2.42e+6"; parameters["mooneyrivlin_a"] = "2.42e+6";
parameters["mooneyrivlin_b"] = "6.52e+6"; parameters["mooneyrivlin_b"] = "6.52e+6";
parameters["mooneyrivlin_c"] = "-7.34e+6"; parameters["mooneyrivlin_c"] = "-7.34e+6";
expectedEnergy = 83069427.2; expectedEnergy = 68753239.6;
expectedGradientTwoNorm = 163210957; expectedGradientTwoNorm = 31244643.1;
expectedGradientInfinityNorm = 49284369.2; expectedGradientInfinityNorm = 9673572.39;
expectedMatrixFrobeniusNorm = 5.92126562e+09; expectedMatrixFrobeniusNorm = 1.67660965e+09;
int testCiarlet = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm); int testCiarlet = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
return testCiarlet + testLog + testSquare; return testCiarlet + testLog + testSquare;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment