diff --git a/dune/elasticity/assemblers/feassembler.hh b/dune/elasticity/assemblers/feassembler.hh
index 8e7c7fcb2c8506cab9f4d824724e0d6f56a17d76..eae94ae0c359f499a5c2d01a80222e32e1a06c43 100644
--- a/dune/elasticity/assemblers/feassembler.hh
+++ b/dune/elasticity/assemblers/feassembler.hh
@@ -13,6 +13,8 @@
#include <dune/functions/functionspacebases/concepts.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/grid/common/partitionset.hh>
+
#include "localfestiffness.hh"
diff --git a/dune/elasticity/materials/mooneyrivlindensity.hh b/dune/elasticity/materials/mooneyrivlindensity.hh
index 6871bc92b168156a50038d23d92057ba7d7f9339..3cf4d5d591a863bcacc5c300688066d516c60cd1 100644
--- a/dune/elasticity/materials/mooneyrivlindensity.hh
+++ b/dune/elasticity/materials/mooneyrivlindensity.hh
@@ -58,53 +58,62 @@ public:
for (int k=0; k<dim; k++)
C[i][j] += gradient[k][i] * gradient[k][j];
- //////////////////////////////////////////////////////////
- // Eigenvalues of C = F^TF
- //////////////////////////////////////////////////////////
-
-
- // eigenvalues of C, i.e., squared singular values \sigma_i of F
- Dune::FieldVector<field_type, dim> sigmaSquared;
- FMatrixHelp::eigenValues(C, sigmaSquared);
-
- field_type normFSquared = gradient.frobenius_norm2();
+ /* The Mooney-Rivlin-Density is given as a function of the eigenvalues of the right Cauchy-Green-Deformation tensor C = F^TF
+ or the right Cauchy-Green-Deformation tensor B = FF^T.
+ C = F^TF and B = FF^T have the same eigenvalues - we can use either one of them.
+
+ There are three Mooney-Rivlin-Variants:
+ ciarlet: W(F) = mooneyrivlin_a*(normF)^2 + mooneyrivlin_b*(normFinv)^2*detF + mooneyrivlin_c*(detF)^2 -
+ ((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))
+ 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)
+
+ For log and square: I1 and I2 are the first two invariants of C (or B), multiplied with a factor depending on det(F):
+ I1 = (det(F)^(-2/dim)) * [ first invariant of C ]
+ = (det(F)^(-2/dim)) * (sum of all eigenvalues of C)
+ = (det(F)^(-2/dim)) * trace(C)
+ = (det(F)^(-2/dim)) * trace(F^TF)
+ = (det(F)^(-2/dim)) * (frobenius_norm(F))^2
+ I2 = (det(F)^(-4/dim)) * [ second invariant of C ]
+ = (det(F)^(-4/dim)) * 0.5 * [ trace(C)^2 - tr(C^2) ]
+ = (det(F)^(-4/dim)) * [C_11C_22 + C_11C_33 + C_22C_33 - C_12C_21 - C_13C_31 - C_23C_32]
+ = (det(F)^(-4/dim)) * [C_11C_22 + C_11C_33 + C_22C_33 - C_12^2 - C_13^2 - C_23^2]
+ */
+
+ field_type frobeniusNormFsquared = gradient.frobenius_norm2();
field_type detF = gradient.determinant();
if (detF < 0)
DUNE_THROW( Dune::Exception, "det(F) < 0, so it is not possible to calculate the MooneyRivlinEnergy.");
- field_type normFinvSquared = 0;
-
- field_type c2Tilde = 0;
- for (int i = 0; i < dim; i++) {
- normFinvSquared += 1/sigmaSquared[i];
- // compute D, which is the sum of the squared eigenvalues
- for (int j = i+1; j < dim; j++)
- c2Tilde += sigmaSquared[j]*sigmaSquared[i];
- }
-
- field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
- // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
- c2Tilde /= pow(detF, 4.0/dim);
- field_type c2TildeMinus3 = c2Tilde - 3;
-
// Add the local energy density
field_type strainEnergy = 0;
if (mooneyrivlin_energy == "ciarlet")
{
+ //To calculate the squared frobeniusnorm of F^(-1), we actually invert F^(-1)
+ auto gradientInverse = gradient;
+ gradientInverse.invert();
+ field_type frobeinusNormFInverseSquared = gradientInverse.frobenius_norm2();
using std::log;
- return mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF);
- } else {
- strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
- mooneyrivlin_01 * c2TildeMinus3 +
- mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 +
- mooneyrivlin_02 * c2TildeMinus3 * c2TildeMinus3 +
- mooneyrivlin_11 * trCTildeMinus3 * c2TildeMinus3 +
- mooneyrivlin_30 * trCTildeMinus3 * trCTildeMinus3 * trCTildeMinus3 +
- mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3 +
- mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3 * c2TildeMinus3 +
- mooneyrivlin_03 * c2TildeMinus3 * c2TildeMinus3 * c2TildeMinus3;
+ return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF);
+ } else { // mooneyrivlin_energy is "log" or "square"
+ field_type a = pow(detF, 2.0/dim);
+ field_type invariant1Minus3 = frobeniusNormFsquared/a - 3;
+ field_type secondInvariantOfC = 0;
+ for (int i=0; i<dim-1; i++)
+ for (int j=i+1; j<dim; j++)
+ secondInvariantOfC += C[i][i]*C[j][j] - C[i][j]*C[i][j];
+ field_type invariant2Minus3 = secondInvariantOfC/(a*a) - 3;
+ strainEnergy = mooneyrivlin_10 * invariant1Minus3 +
+ mooneyrivlin_01 * invariant2Minus3 +
+ mooneyrivlin_20 * invariant1Minus3 * invariant1Minus3 +
+ mooneyrivlin_02 * invariant2Minus3 * invariant2Minus3 +
+ mooneyrivlin_11 * invariant1Minus3 * invariant2Minus3 +
+ mooneyrivlin_30 * invariant1Minus3 * invariant1Minus3 * invariant1Minus3 +
+ mooneyrivlin_21 * invariant1Minus3 * invariant1Minus3 * invariant2Minus3 +
+ mooneyrivlin_12 * invariant1Minus3 * invariant2Minus3 * invariant2Minus3 +
+ mooneyrivlin_03 * invariant2Minus3 * invariant2Minus3 * invariant2Minus3;
if (mooneyrivlin_energy == "log") {
using std::log;
field_type logDetF = log(detF);
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 65035164e33cd40989d38618a114079d84a4bb14..262d1474df7f0cc6dde019e24b4fee36f93fb729 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -8,5 +8,8 @@ dune_add_test(SOURCES materialtest
add_dune_ug_flags(materialtest)
add_dune_adolc_flags(materialtest)
+dune_add_test(SOURCES mooneyrivlintest
+ CMAKE_GUARD ADOLC_FOUND)
+
dune_add_test(SOURCES localhyperdualstiffnesstest
CMAKE_GUARD ADOLC_FOUND)
diff --git a/test/mooneyrivlintest.cc b/test/mooneyrivlintest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b69f5eb4151ad60a4236740d79004266a0ef1718
--- /dev/null
+++ b/test/mooneyrivlintest.cc
@@ -0,0 +1,176 @@
+#include <config.h>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/elasticity/assemblers/localadolcstiffness.hh>
+#include <dune/elasticity/assemblers/feassembler.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
+#include <dune/elasticity/materials/mooneyrivlindensity.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+// grid dimension
+const int dim = 3;
+
+//order of integration
+const int order = 1;
+
+using namespace Dune;
+
+//differentiation method: ADOL-C
+typedef adouble ValueType;
+
+using namespace Dune;
+
+typedef BlockVector<FieldVector<double,dim> > VectorType;
+typedef BCRSMatrix<FieldMatrix<double, dim, dim> > MatrixType;
+
+template <class Basis>
+int assembleAndCompare (const Basis basis, const ParameterTree parameters, const VectorType x,
+ double expectedEnergy, double expectedGradientTwoNorm, double expectedGradientInfinityNorm, double expectedMatrixFrobeniusNorm) {
+ using LocalView = typename Basis::LocalView;
+
+ std::string mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy");
+
+ auto elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(parameters);
+ auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, adouble>>(elasticDensity);
+ auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<LocalView>>(elasticEnergy);
+ Elasticity::FEAssembler<Basis,VectorType> assembler(basis, localADOLCStiffness);
+
+ //////////////////////////////////////////////////////////////
+ // Compute the energy, assemble and compare!
+ //////////////////////////////////////////////////////////////
+
+ VectorType gradient;
+ MatrixType hessianMatrix;
+ double energy = assembler.computeEnergy(x);
+
+ if ( std::abs(energy - expectedEnergy)/expectedEnergy > 1e-3)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "In " << mooneyrivlin_energy << ": Calculated energy is " << energy << " but '" << expectedEnergy << "' was expected!" << std::endl;
+ return 1;
+ }
+
+ assembler.assembleGradientAndHessian(x, gradient, hessianMatrix, true);
+
+ double gradientTwoNorm = gradient.two_norm();
+ double gradientInfinityNorm = gradient.infinity_norm();
+ double matrixFrobeniusNorm = hessianMatrix.frobenius_norm();
+ if ( std::abs(gradientTwoNorm - expectedGradientTwoNorm)/expectedGradientTwoNorm > 1e-4 || std::abs(gradientInfinityNorm - expectedGradientInfinityNorm)/expectedGradientInfinityNorm > 1e-8)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "In " << mooneyrivlin_energy << ": Calculated gradient max norm is " << gradientInfinityNorm << " but '" << expectedGradientInfinityNorm << "' was expected!" << std::endl;
+ std::cerr << "In " << mooneyrivlin_energy << ": Calculated gradient norm is " << gradientTwoNorm << " but '" << expectedGradientTwoNorm << "' was expected!" << std::endl;
+ return 1;
+ }
+
+ if ( std::abs(matrixFrobeniusNorm - expectedMatrixFrobeniusNorm)/expectedMatrixFrobeniusNorm > 1e-4)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "In " << mooneyrivlin_energy << ": Calculated matrix norm is " << matrixFrobeniusNorm << " but '" << expectedMatrixFrobeniusNorm << "' was expected!" << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
+
+int main (int argc, char *argv[])
+{
+
+ Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+
+ // ///////////////////////////////////////
+ // Create the grid
+ // ///////////////////////////////////////
+ using GridType = UGGrid<dim>;
+ auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0,0}, {1,1,1}, {3,3,3});
+ int refine = 3;
+ while(refine > 0) {
+ for (auto&& e : elements(grid->leafGridView())){
+ bool refine = false;
+ for (int i = 0; i < e.geometry().corners(); i++) {
+ refine = refine || (e.geometry().corner(i)[0] > 0.99);
+ }
+ grid->mark(refine ? 1 : 0, e);
+ }
+ grid->adapt();
+ refine--;
+ }
+ grid->loadBalance();
+
+ using GridView = GridType::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create a power basis
+ /////////////////////////////////////////////////////////////////////////
+
+ using namespace Functions::BasisFactory;
+ auto basis = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+
+ using PowerBasis = decltype(basis);
+
+ //////////////////////////////////////////////////////////////
+ // Prepare the iterate x where we want to assemble
+ //////////////////////////////////////////////////////////////
+ VectorType x(basis.size());
+ Functions::interpolate(basis, x, [](FieldVector<double,dim> x){
+ FieldVector<double,dim> y = {0.5*x[0]*x[0], 0.5*x[1] + 0.2*x[0], 4*std::abs(std::sin(x[2]))};
+ return y;
+ });
+
+ ParameterTree parameters;
+ parameters["mooneyrivlin_10"] = "1.67e+6";
+ parameters["mooneyrivlin_01"] = "1.94e+6";
+ parameters["mooneyrivlin_20"] = "2.42e+6";
+ parameters["mooneyrivlin_02"] = "6.52e+6";
+ parameters["mooneyrivlin_11"] = "-7.34e+6";
+ parameters["mooneyrivlin_30"] = "0";
+ parameters["mooneyrivlin_03"] = "0";
+ parameters["mooneyrivlin_21"] = "0";
+ parameters["mooneyrivlin_12"] = "0";
+ parameters["mooneyrivlin_k"] = "75e+6";
+
+ parameters["mooneyrivlin_energy"] = "square";
+ double expectedEnergy = 177758377;
+ double expectedGradientTwoNorm = 2.13791422e+09;
+ double expectedGradientInfinityNorm = 599188506;
+ double expectedMatrixFrobeniusNorm = 1.62529884e+11;
+ int testSquare = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+ parameters["mooneyrivlin_energy"] = "log";
+ expectedEnergy = 181180273;
+ expectedGradientTwoNorm = 2.29399362e+09;
+ expectedGradientInfinityNorm = 648551554;
+ expectedMatrixFrobeniusNorm = 1.67152642e+11;
+ int testLog = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+ parameters["mooneyrivlin_energy"] = "ciarlet";
+ parameters["mooneyrivlin_a"] = "2.42e+6";
+ parameters["mooneyrivlin_b"] = "6.52e+6";
+ parameters["mooneyrivlin_c"] = "-7.34e+6";
+ expectedEnergy = 83069427.2;
+ expectedGradientTwoNorm = 163210957;
+ expectedGradientInfinityNorm = 49284369.2;
+ expectedMatrixFrobeniusNorm = 5.92126562e+09;
+ int testCiarlet = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+ return testCiarlet + testLog + testSquare;
+}