Implement general surface normal fields for shells with nonplanar reference geometry

Previously, the NonplanarCosseratShellEnergy had the normal field of the unit sphere
hard-wired in the code.  With this patch, using general triangulated surface as
reference geometries is possible.

dune/gfe/nonplanarcosseratshellenergy.hh

@@ -181,10 +181,12 @@ public:
    * \param parameters The material parameters
    */
   NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
+                               const std::vector<UnitVector<double,3> >& vertexNormals,
                                const BoundaryPatch<GridView>* neumannBoundary,
                                const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction,
                                const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
-  : neumannBoundary_(neumannBoundary),
+  : vertexNormals_(vertexNormals),
+    neumannBoundary_(neumannBoundary),
     neumannFunction_(neumannFunction),
     volumeLoad_(volumeLoad)
   {
@@ -248,6 +250,9 @@ public:
   /** \brief Curvature parameters */
   double b1_, b2_, b3_;
 
+  /** \brief The normal vectors at the grid vertices.  This are used to compute the reference surface curvature. */
+  std::vector<UnitVector<double,3> > vertexNormals_;
+
   /** \brief The Neumann boundary */
   const BoundaryPatch<GridView>* neumannBoundary_;
 
@@ -274,14 +279,12 @@ energy(const typename Basis::LocalView& localView,
   ////////////////////////////////////////////////////////////////////////////////////
   //  Construct a linear (i.e., non-constant!) normal field on the surface
   ////////////////////////////////////////////////////////////////////////////////////
+  auto gridView = localView.globalBasis().gridView();
 
+  assert(vertexNormals_.size() == gridView.indexSet().size(gridDim));
   std::vector<UnitVector<double,3> > cornerNormals(element.subEntities(gridDim));
   for (size_t i=0; i<cornerNormals.size(); i++)
-  {
-    // TODO: WARNING WARNING WARNING:  Shape is hard-coded here!
-    // The sphere:
-    cornerNormals[i] = geometry.corner(i) / geometry.corner(i).two_norm();
-  }
+    cornerNormals[i] = vertexNormals_[gridView.indexSet().subIndex(element,i,2)];
 
   typedef typename Dune::PQkLocalFiniteElementCache<DT, double, gridDim, 1> P1FiniteElementCache;
   typedef typename P1FiniteElementCache::FiniteElementType P1LocalFiniteElement;

dune/gfe/vertexnormals.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GFE_VERTEXNORMALS_HH
+#define DUNE_GFE_VERTEXNORMALS_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/gfe/unitvector.hh>
+
+/** \brief Compute averaged vertex normals for a 2d-in-3d grid
+ */
+template <class GridView>
+std::vector<UnitVector<typename GridView::ctype,3> > computeVertexNormals(const GridView& gridView)
+{
+  const auto& indexSet = gridView.indexSet();
+
+  assert(GridView::dimension == 2);
+
+  std::vector<Dune::FieldVector<typename GridView::ctype,3> > unscaledNormals(indexSet.size(2));
+  std::fill(unscaledNormals.begin(), unscaledNormals.end(), Dune::FieldVector<typename GridView::ctype,3>(0.0));
+
+  for (const auto& element : elements(gridView))
+  {
+    for (int i=0; i<element.subEntities(2); i++)
+    {
+      auto cornerPos = Dune::ReferenceElements<double,2>::general(element.type()).position(i,2);
+      auto tangent = element.geometry().jacobianTransposed(cornerPos);
+      auto cornerNormal = Arithmetic::crossProduct(tangent[0], tangent[1]);
+      cornerNormal /= cornerNormal.two_norm();
+
+      unscaledNormals[indexSet.subIndex(element,i,2)] += cornerNormal;
+    }
+  }
+
+  // Normalize
+  for (auto& normal : unscaledNormals)
+    normal /= normal.two_norm();
+
+  // Convert array of vectors to array of UnitVectors
+  std::vector<UnitVector<typename GridView::ctype,3> > result(indexSet.size(2));
+
+  for (std::size_t i=0; i<unscaledNormals.size(); i++)
+    result[i] = unscaledNormals[i];
+
+  return result;
+}
+
+
+#endif
\ No newline at end of file

src/cosserat-continuum.cc

@@ -43,6 +43,7 @@
 #include <dune/gfe/vtkreader.hh>
 #include <dune/gfe/geodesicfeassembler.hh>
 #include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/vertexnormals.hh>
 
 // grid dimension
 const int dim = 2;
@@ -311,7 +312,9 @@ int main (int argc, char *argv[]) try
     }
     else
     {
+      std::vector<UnitVector<double,3> > vertexNormals = computeVertexNormals(gridView);
       cosseratEnergyADOLCLocalStiffness = std::make_shared<NonplanarCosseratShellEnergy<FEBasis,3,adouble> >(materialParameters,
+                                                                                                             std::move(vertexNormals),
                                                                                                              &neumannBoundary,
                                                                                                              neumannFunction,
                                                                                                              volumeLoad);