Merge branch 'fix-normalDerivative' into 'master'

Use correct second fundamental tensor in nonplanarcosseratshellenergy and surfacecosseratenergy

See merge request !88

.gitlab-ci.yml

@@ -51,7 +51,7 @@ dune:git clang:
   - *before
   script: duneci-standard-test
 
-dune:git dune-parmg dune-vtk dune-curvedgeometry gcc:
+dune:git dune-parmg dune-vtk dune-curvedgeometry dune-curvedgrid gcc:
   image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17
   before_script:
   - *patch-dune-common
@@ -59,9 +59,10 @@ dune:git dune-parmg dune-vtk dune-curvedgeometry gcc:
   - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
   - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-vtk.git
   - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
+  - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
   script: duneci-standard-test
 
-dune:git dune-parmg dune-vtk dune-curvedgeometry clang:
+dune:git dune-parmg dune-vtk dune-curvedgeometry dune-curvedgrid clang:
   image: registry.dune-project.org/docker/ci/dune:git-ubuntu-20.04-clang-10-20
   before_script:
   - *patch-dune-common
@@ -69,4 +70,5 @@ dune:git dune-parmg dune-vtk dune-curvedgeometry clang:
   - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
   - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-vtk.git
   - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
+  - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
   script: duneci-standard-test

dune/gfe/nonplanarcosseratshellenergy.hh

@@ -146,6 +146,9 @@ energy(const typename Basis::LocalView& localView,
   // The element geometry
   auto element = localView.element();
 
+  // The set of shape functions on this element
+  const auto& localFiniteElement = localView.tree().finiteElement();
+
 #if HAVE_DUNE_CURVEDGEOMETRY
   // Construct a curved geometry of this element of the Cosserat shell in stress-free state
   // When using element.geometry(), then the curvatures on the element are zero, when using a curved geometry, they are not
@@ -153,6 +156,7 @@ energy(const typename Basis::LocalView& localView,
   // this is used for the curved geometry approximation.
   // The variable local holds the local coordinates in the reference element
   // and localGeometry.global maps them to the world coordinates
+  auto curvedGeometryGridFunctionOrder = localFiniteElement.localBasis().order();
   Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim>>> geometry(referenceElement(element),
     [this,element](const auto& local) {
       if (not stressFreeStateGridFunction_) {
@@ -161,14 +165,11 @@ energy(const typename Basis::LocalView& localView,
       auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
       localGridFunction.bind(element);
       return localGridFunction(local);
-    }, 2); /*order*/
+    }, curvedGeometryGridFunctionOrder);
 #else
   auto geometry = element.geometry();
 #endif
 
-  // The set of shape functions on this element
-  const auto& localFiniteElement = localView.tree().finiteElement();
-
   ////////////////////////////////////////////////////////////////////////////////////
   //  Set up the local nonlinear finite element function
   ////////////////////////////////////////////////////////////////////////////////////
@@ -192,7 +193,7 @@ energy(const typename Basis::LocalView& localView,
     // The value of the local function
     RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
 
-    // The derivative of the local function
+    // The derivative of the local function w.r.t. the coordinate system of the tangent space
     auto derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
 
     //////////////////////////////////////////////////////////
@@ -204,6 +205,7 @@ energy(const typename Basis::LocalView& localView,
     value.q.matrix(R);
     auto RT = Dune::GFE::transpose(R);
 
+    //Derivative of the rotation w.r.t. the coordinate system of the tangent space
     Tensor3<field_type,3,3,gridDim> DR = value.quaternionTangentToMatrixTangent(derivative);
 
     //////////////////////////////////////////////////////////
@@ -252,29 +254,25 @@ energy(const typename Basis::LocalView& localView,
       for (int beta=0; beta<2; beta++)
         c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
 
+    Dune::FieldMatrix<double,3,3> b(0);
 #if HAVE_DUNE_CURVEDGEOMETRY
-    // Second fundamental form: The derivative of the normal field, on each quadrature point
-    auto normalDerivative = geometry.normalGradient(quad[pt].position());
-#else
-    //In case dune-curvedgeometry is not installed, the normal derivative is set to zero.
-    Dune::FieldMatrix<double,3,3> normalDerivative(0);
+    // In case dune-curvedgeometry is not installed, we assume the second fundamental form b
+    // ( (-1) * the derivative of the normal field, evaluated at each quadrature point) is zero!
+    auto grad_s_n = geometry.normalGradient(quad[pt].position());
+    grad_s_n *= (-1);
+    for (int i = 0; i < dimworld; i++)
+      for (int j = 0; j < dimworld; j++)
+        b[i][j] = grad_s_n[i][j];
 #endif
-
-    Dune::FieldMatrix<double,3,3> b(0);
-    for (int alpha=0; alpha<gridDim; alpha++)
-    {
-      Dune::FieldVector<double,3> vec;
-      for (int i=0; i<3; i++)
-        vec[i] = normalDerivative[i][alpha];
-      b -= Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
-    }
-
-    // Gauss curvature
-    auto K = b.determinant();
-
+    
     // Mean curvatue
     auto H = 0.5 * Dune::GFE::trace(b);
 
+    // Gauss curvature, calculated with the normalGradient in the euclidean coordinate system
+    // see e.g. formula (3.5) in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
+    auto bSquared = b*b;
+    auto K = 2*H*H - 0.5*Dune::GFE::trace(bSquared);
+
     //////////////////////////////////////////////////////////
     //  Strain tensors
     //////////////////////////////////////////////////////////
@@ -292,7 +290,7 @@ energy(const typename Basis::LocalView& localView,
 
     Ee = RT * grad_s_m - a;
 
-    // Elastic shell bending-curvature strain
+    // Elastic shell bending-curvature strain, e.g. formula (4.56) in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
     Dune::FieldMatrix<field_type,3,3> Ke(0);
     for (int alpha=0; alpha<gridDim; alpha++)
     {

dune/gfe/surfacecosseratenergy.hh

@@ -286,23 +286,17 @@ RT energy(const typename Basis::LocalView& localView,
           c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
 
       // Second fundamental form: The derivative of the normal field
-      auto normalDerivative = boundaryGeometry.normalGradient(quad[pt].position());
-
-      Dune::FieldMatrix<double,3,3> b(0);
-      for (int alpha=0; alpha<boundaryDim; alpha++)
-      {
-        Dune::FieldVector<double,3> vec;
-        for (int i=0; i<3; i++)
-          vec[i] = normalDerivative[i][alpha];
-        b -= Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
-      }
-
-      // Gauss curvature
-      auto K = b.determinant();
+      auto b = boundaryGeometry.normalGradient(quad[pt].position());
+      b *= (-1);
 
       // Mean curvatue
       auto H = 0.5 * Dune::GFE::trace(b);
 
+      // Gauss curvature, calculated with the normalGradient in the eucidean coordinate system
+      // see e.g. formula (3.5) from "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
+      auto bSquared = b*b;
+      auto K = 2*H*H - 0.5*Dune::GFE::trace(bSquared);
+
       //////////////////////////////////////////////////////////
       //  Strain tensors
       //////////////////////////////////////////////////////////

test/CMakeLists.txt

@@ -18,6 +18,10 @@ if (${DUNE_COMMON_VERSION} VERSION_GREATER_EQUAL 2.8)
   set(TESTS cosseratrodtest ${TESTS})
 endif()
 
+if (dune-curvedgrid_FOUND)
+  set(TESTS test-cylinder ${TESTS})
+endif()
+
 if (${DUNE_COMMON_VERSION} VERSION_GREATER_EQUAL 2.8)
   set(TESTS surfacecosseratstressassemblertest ${TESTS})
 endif()

test/test-cylinder.cc

+#include <config.h>
+
+#include <fenv.h>
+#include <array>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/curvedgeometry/curvedgeometry.hh>
+#include <dune/curvedgrid/grid.hh>
+#include <dune/curvedgrid/gridfunctions/analyticdiscretefunction.hh>
+#include <dune/curvedgrid/geometries/cylinder.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/gfe/linearalgebra.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+
+// grid dimension
+const int gridDim = 2;
+const int dimWorld = 3;
+const int approximationOrderGeometry = 4;
+const int approximationOrderAnalytics = 6;
+const int elementsCircle = 9;
+
+const auto pi = std::acos(-1.0);
+
+using namespace Dune;
+
+int main (int argc, char *argv[])
+{
+    MPIHelper::instance(argc,argv);
+    // Cylinder
+    static const double radius            = 4;
+    static const double height            = 2;
+    static const double cylinderFraction  = 0.75;
+
+    /////////////////////////////////////////////
+    //    Create the grid for the cylinder
+    /////////////////////////////////////////////
+
+    struct CylinderCreator
+            : public Dune :: AnalyticalCoordFunction< double, 2, 3, CylinderCreator >
+    {
+        FieldVector<double,3> operator() ( const FieldVector<double, 2> &x ) const
+        {
+            FieldVector<double,3> y;
+            y[0] = radius * std::cos(x[0]);
+            y[1] = radius * std::sin(x[0]);
+            y[2] = x[1];
+            return y;
+        }
+    };
+
+    using GridType = Dune::GeometryGrid< YaspGrid<gridDim>, CylinderCreator>;
+    std::shared_ptr<YaspGrid<gridDim>> hostGrid;
+    hostGrid = Dune::StructuredGridFactory<Dune::YaspGrid<2>>::createCubeGrid({0, 0}, {cylinderFraction*2*pi, height}, {elementsCircle,2});
+    auto cylinderCreator = std::make_shared<CylinderCreator>();
+    auto grid = std::make_shared<GridType>(hostGrid, cylinderCreator);
+    auto gridView = grid->leafGridView();
+
+    auto cylinder = CylinderProjection<double>{radius};
+    auto polynomialCylinderGF = analyticDiscreteFunction(cylinder, *grid, approximationOrderAnalytics);
+
+    auto analyticCylinderGF = cylinderGridFunction<GridType>(radius);
+    auto cylinderLocalFunction = localFunction(analyticCylinderGF);
+
+    auto quadOrder = 10;
+    for (const auto& element : elements(gridView, Dune::Partitions::interior)) {
+        using DT = decltype(gridView)::ctype;
+
+        Dune::CurvedGeometry<DT, gridDim, dimWorld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim>>>
+        polynomialGeometry(Dune::referenceElement(element.geometry()), [element, polynomialCylinderGF](const auto& local) {
+            auto localGridFunction = localFunction(polynomialCylinderGF);
+            localGridFunction.bind(element);
+            return localGridFunction(local);
+        }, approximationOrderGeometry);
+
+        cylinderLocalFunction.bind(element);
+        auto localGeometry = Dune::DefaultLocalGeometry<double,2,2>{};
+        auto analyticGeometry = Dune::LocalFunctionGeometry{element.type(), cylinderLocalFunction, localGeometry};
+        Dune::LagrangeLFECache<DT,DT,gridDim> cache(approximationOrderAnalytics);
+        auto refEle = referenceElement(element);
+        auto lFE = cache.get(refEle.type());
+
+        const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+        for (size_t pt=0; pt<quad.size(); pt++) {
+            // Check if mean curvature is correct
+            auto realMeanCurvature = std::abs(cylinder.mean_curvature(polynomialGeometry.global(quad[pt].position())));
+
+            auto normalDerivativeP = polynomialGeometry.normalGradient(quad[pt].position());
+            auto approximatedCurvatureP = 0.5*std::abs(Dune::GFE::trace(normalDerivativeP));
+            auto relativeDifferenceP = std::abs((realMeanCurvature - approximatedCurvatureP)/realMeanCurvature);
+
+            auto normalDerivativeA = analyticGeometry.normalGradient(quad[pt].position(), lFE);
+            auto approximatedCurvatureA = 0.5*std::abs(Dune::GFE::trace(normalDerivativeA));
+            std::cout << approximatedCurvatureA << std::endl;
+            auto relativeDifferenceA = std::abs((realMeanCurvature - approximatedCurvatureA)/realMeanCurvature);
+
+            if (relativeDifferenceP > 0.005){
+                std::cerr << "At point " << polynomialGeometry.global(quad[pt].position()) << " the curvature (approximated using a polynomial) is "
+                    << approximatedCurvatureP << " but " << realMeanCurvature << " was expected!" << std::endl;
+                return 1;
+            }
+            if (relativeDifferenceA > 0.005){
+                std::cerr << "At point " << polynomialGeometry.global(quad[pt].position()) << " the curvature (approximated using the real analytic cylinder function) is "
+                    << approximatedCurvatureA << " but " << realMeanCurvature << " was expected!" << std::endl;
+                return 1;
+            }
+        }
+    }
+}