From d7c156ff9c47cbce98c6e2cd90ee5597882bfb0f Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 2 Jan 2025 12:23:37 +0100
Subject: [PATCH] CosseratShellDensity: Depend on grid element or intersection

Currently, densities only get the type used for coordinates
of the integration domain.  However, this is not enough:
The complete domain must be known, and this domain can be either
a grid element or a grid intersection.

As a first small step, this patch introduces ElementOrIntersection
as a template parameter of the CosseratShellDensity class.
---
 dune/gfe/assemblers/nonplanarcosseratshellenergy.hh | 8 +++-----
 dune/gfe/assemblers/surfacecosseratenergy.hh        | 6 ++++--
 dune/gfe/densities/cosseratshelldensity.hh          | 9 +++++++--
 src/cosserat-continuum.cc                           | 4 ++--
 src/film-on-substrate.cc                            | 4 ++--
 test/filmonsubstratetest.cc                         | 4 ++--
 test/nonplanarcosseratenergytest.cc                 | 4 ++--
 7 files changed, 22 insertions(+), 17 deletions(-)

diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index 0dc07d7d..aba5964a 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -79,21 +79,19 @@ class NonplanarCosseratShellEnergy
   typedef typename GridView::ctype DT;
   typedef Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > TargetSpace;
   typedef typename TargetSpace::ctype RT;
-  typedef typename GridView::template Codim<0>::Entity Entity;
+  typedef typename GridView::template Codim<0>::Entity Element;
 
   // some other sizes
   constexpr static int gridDim = GridView::dimension;
   constexpr static int dimworld = GridView::dimensionworld;
 
-  using Position = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
-
 public:
 
   /** \brief Constructor with a set of material parameters
    * \param parameters                  The material parameters
    * \param stressFreeStateGridFunction Pointer to a parametrization representing the Cosserat shell in a stress-free state
    */
-  NonplanarCosseratShellEnergy(const std::shared_ptr<Dune::GFE::CosseratShellDensity<Position,field_type> >& density,
+  NonplanarCosseratShellEnergy(const std::shared_ptr<Dune::GFE::CosseratShellDensity<Element,field_type> >& density,
                                const StressFreeStateGridFunction* stressFreeStateGridFunction)
     : stressFreeStateGridFunction_(stressFreeStateGridFunction),
     density_(density)
@@ -123,7 +121,7 @@ public:
   const StressFreeStateGridFunction* stressFreeStateGridFunction_;
 
   /** \brief The energy density of a Cosserat shell with nonplanar reference configuration */
-  const std::shared_ptr<Dune::GFE::CosseratShellDensity<Position,field_type> > density_;
+  const std::shared_ptr<Dune::GFE::CosseratShellDensity<Element,field_type> > density_;
 };
 
 template <class Basis, int dim, class field_type, class StressFreeStateGridFunction>
diff --git a/dune/gfe/assemblers/surfacecosseratenergy.hh b/dune/gfe/assemblers/surfacecosseratenergy.hh
index b9c57d73..bbd18f12 100644
--- a/dune/gfe/assemblers/surfacecosseratenergy.hh
+++ b/dune/gfe/assemblers/surfacecosseratenergy.hh
@@ -37,6 +37,8 @@ namespace Dune::GFE
     constexpr static int gridDim = GridView::dimension;
     static constexpr int boundaryDim = gridDim - 1;
 
+    using Intersection = typename GridView::Intersection;
+
   public:
 
     /** \brief Constructor with a set of material parameters
@@ -45,7 +47,7 @@ namespace Dune::GFE
      * \param curvedGeometryGridFunction The curvedGeometryGridFunction gives the geometry of the shell in stress-free state.
               When assembling, we deform the intersections using the curvedGeometryGridFunction and then use the deformed geometries.
      */
-    SurfaceCosseratEnergy(const std::shared_ptr<CosseratShellDensity<FieldVector<DT,3>,RT> >& density,
+    SurfaceCosseratEnergy(const std::shared_ptr<CosseratShellDensity<Intersection,RT> >& density,
                           const BoundaryPatch<GridView>* shellBoundary,
                           const CurvedGeometryGridFunction& curvedGeometryGridFunction)
       : shellBoundary_(shellBoundary),
@@ -199,7 +201,7 @@ namespace Dune::GFE
     const CurvedGeometryGridFunction curvedGeometryGridFunction_;
 
     /** \brief The density that is being integrated */
-    const std::shared_ptr<CosseratShellDensity<FieldVector<DT,3>,RT> > density_;
+    const std::shared_ptr<CosseratShellDensity<Intersection,RT> > density_;
 
   };
 }  // namespace Dune::GFE
diff --git a/dune/gfe/densities/cosseratshelldensity.hh b/dune/gfe/densities/cosseratshelldensity.hh
index 447b5090..a8027211 100644
--- a/dune/gfe/densities/cosseratshelldensity.hh
+++ b/dune/gfe/densities/cosseratshelldensity.hh
@@ -24,16 +24,21 @@ namespace Dune::GFE
    *
    * (and earlier papers by Neff and Bîrsan).
    *
-   * \tparam Position Type for the points where the density will be evaluated
+   * \tparam ElementOrIntersection The domain of the density.
+   *   Can be either a grid element (i.e., a codimension-0 Entity)
+   *   or an Intersection.
    * \tparam field_type Type used for numbers
    */
-  template<class Position, class field_type>
+  template<class ElementOrIntersection, class field_type>
   class CosseratShellDensity
   {
     constexpr static int dimWorld = 3;
     constexpr static int domainDim = 2;
 
+    using Geometry = typename ElementOrIntersection::Geometry;
+
     // TODO: This is likely not a good long-term solution
+    using Position = Geometry::GlobalCoordinate;
     static_assert(Position::size()==dimWorld, "Position must be in world-coordinates.");
 
     using RigidBodyMotion = ProductManifold<RealTuple<field_type,dimWorld>, Rotation<field_type,dimWorld> >;
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index f656ae22..64c11c35 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -132,8 +132,8 @@ auto createCosseratEnergy(const ParameterTree& materialParameters,
   }
   else if constexpr (LocalCoordinate::size()==2 && dimworld==3)
   {
-    using GlobalCoordinate = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
-    auto density = std::make_shared<GFE::CosseratShellDensity<GlobalCoordinate, adouble> >(materialParameters);
+    using Element = typename GridView::template Codim<0>::Entity;
+    auto density = std::make_shared<GFE::CosseratShellDensity<Element, adouble> >(materialParameters);
 
     return std::make_shared<NonplanarCosseratShellEnergy<Basis, 3, adouble, decltype(creator)> >(density, &creator);
   }
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index e9684c2d..e822c73f 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -531,8 +531,8 @@ int main (int argc, char *argv[]) try
     auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunctionPtr);
 
     // The energy of the surface shell
-    auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<
-        FieldVector<double,3>, adouble> >(
+    using Intersection = typename GridView::Intersection;
+    auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<Intersection, adouble> >(
       materialParameters,
       fThickness,
       fLame);
diff --git a/test/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index d3f8b83c..80e88579 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -356,8 +356,8 @@ int main (int argc, char *argv[])
   auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensityWrapped);
   auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunction);
 
-  auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<
-      FieldVector<double,3>, adouble> >(
+  using Intersection = typename GridView::Intersection;
+  auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<Intersection, adouble> >(
     materialParameters,
     fThickness,
     fLame);
diff --git a/test/nonplanarcosseratenergytest.cc b/test/nonplanarcosseratenergytest.cc
index 8f3ab087..6566c170 100644
--- a/test/nonplanarcosseratenergytest.cc
+++ b/test/nonplanarcosseratenergytest.cc
@@ -116,8 +116,8 @@ double calculateEnergy(const FlatGridView& flatGridView,
   //  Construct the energy functional
   ///////////////////////////////////////////////////
 
-  using GlobalCoordinate = typename FlatGridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
-  auto density = std::make_shared<GFE::CosseratShellDensity<GlobalCoordinate, double> >(materialParameters);
+  using Element = typename FlatGridView::template Codim<0>::Entity;
+  auto density = std::make_shared<GFE::CosseratShellDensity<Element, double> >(materialParameters);
 
   using ShellEnergy = NonplanarCosseratShellEnergy<FlatFEBasis,
       3,                                               // Dimension of the target space
-- 
GitLab