From e299e3c689390aaf82669620feab1b6d1b4de2b0 Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Mon, 18 Oct 2021 09:18:09 +0200
Subject: [PATCH] Use correct second fundamental tensor in
nonplanarcosseratshellenergy and surfacecosseratenergy
---
dune/gfe/nonplanarcosseratshellenergy.hh | 46 +++++----
dune/gfe/surfacecosseratenergy.hh | 20 ++--
test/CMakeLists.txt | 4 +
test/test-cylinder.cc | 114 +++++++++++++++++++++++
4 files changed, 147 insertions(+), 37 deletions(-)
create mode 100644 test/test-cylinder.cc
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index 97a1b6bd..4a5daa42 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/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 "Reï¬ned 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 "Reï¬ned 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++)
{
diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index 7ece7fe2..a9d55824 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/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 "Reï¬ned 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
//////////////////////////////////////////////////////////
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 0491fac8..2f725374 100644
--- a/test/CMakeLists.txt
+++ b/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()
diff --git a/test/test-cylinder.cc b/test/test-cylinder.cc
new file mode 100644
index 00000000..76419928
--- /dev/null
+++ b/test/test-cylinder.cc
@@ -0,0 +1,114 @@
+#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;
+ }
+ }
+ }
+}
--
GitLab