initial version of curvedsurfacegrid based on curvedgeometry

dune.module

@@ -7,5 +7,5 @@ Module: dune-curvedsurfacegrid
 Version: 0.1.0
 Maintainer: florian.stenger@tu-dresden.de
 #depending on
-Depends: dune-common dune-geometry dune-grid dune-localfunctions dune-curvilineargeometry
+Depends: dune-common dune-geometry dune-grid dune-localfunctions dune-curvedgeometry
 Suggests: dune-alugrid dune-foamgrid

dune/curvedsurfacegrid/coordprovider.hh

@@ -5,6 +5,9 @@
 
 #include <array>
 
+#include <dune/localfunctions/lagrange/equidistantpoints.hh>
+#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
+
 namespace Dune
 {
 
@@ -220,26 +223,37 @@ namespace Dune
 
       typedef typename HostGrid::template Codim< codimension >::Entity HostEntity;
 
+      using LocalFECache = LagrangeLocalFiniteElementCache<ctype, ctype, mydimension, Grid::order>;
+      using LocalFiniteElement = typename LocalFECache::FiniteElementType;
+
     public:
       CoordProvider ( const HostEntity &hostEntity,
                     const CoordFunction &coordFunction )
         : coordFunction_(coordFunction),
-          hostEntity_(hostEntity)
+          hostEntity_(hostEntity),
+          localFECache_(),
+          localFE_(localFECache_.get(hostEntity.type()))
       {}
 
       void calculate ( std::vector< Coordinate > &vertices ) const
       {
-        const std::size_t numCorners = hostEntity_.geometry().corners();
-        std::array<Coordinate, (1 << mydim)> corners;
-        for( std::size_t i = 0; i < numCorners; ++i )
-          corners[i] = hostEntity_.geometry().corner(i);
-        InterpolatoryVerticesGenerator<Coordinate, mydim, Grid::order>::generate(corners, vertices);
-        for(auto &v : vertices) coordFunction_.evaluate(v, v);
+        vertices.resize(localFE_.size());
+
+        auto f = [this,geo = hostEntity_.geometry()](auto const& local) -> Coordinate {
+          Coordinate res;
+          this->coordFunction_.evaluate(geo.global(local), res);
+          return res;
+        };
+
+        localFE_.localInterpolation().interpolate(f, vertices);
       }
 
     private:
       const CoordFunction &coordFunction_;
       const HostEntity &hostEntity_;
+
+      LocalFECache localFECache_;
+      LocalFiniteElement const& localFE_;
     };
 
     //CoordProvider (fake)
@@ -273,6 +287,7 @@ namespace Dune
 
       void calculate ( std::vector< Coordinate > &vertices ) const
       {
+        throw "not implemented";
         const GeometryType type = hostElement_.geometry().type();
         auto refElement = referenceElement< ctype, dimension >( type );
         const std::size_t numCorners = refElement.size( subEntity_, codimension, dimension );
@@ -328,6 +343,7 @@ namespace Dune
 
       void calculate ( std::vector< Coordinate > &vertices ) const
       {
+        throw "not implemented";
         const std::size_t numCorners = hostLocalGeometry_.corners();
         std::array<Coordinate, (1 << mydimension)> corners;
         // $flo: coordFunction is applied here which is in contrast to GeomentryGrid's behaviour!

dune/curvedsurfacegrid/geometry.hh

@@ -7,7 +7,7 @@
 
 #include <dune/common/typetraits.hh>
 
-#include <dune/curvilineargeometry/curvilineargeometry.hh>
+#include <dune/curvedgeometry/curvedgeometry.hh>
 #include <dune/geometry/referenceelements.hh>
 #include <dune/geometry/type.hh>
 
@@ -92,7 +92,7 @@ namespace Dune
       static const int codimension = dimension - mydimension;
 
     protected:
-      typedef CachedCurvilinearGeometry<ctype, mydimension, coorddimension> BasicMapping;
+      typedef CurvedGeometry<mydimension, coorddimension, Grid> BasicMapping;
 
       struct Mapping
         : public BasicMapping
@@ -101,7 +101,7 @@ namespace Dune
 
         template< class CoordVector >
         Mapping ( const GeometryType &type, const CoordVector &coords )
-          : BasicMapping( type, coords, Grid::order ),
+          : BasicMapping( type, coords ),
             refCount_( 0 )
         {}
 
@@ -181,24 +181,21 @@ namespace Dune
       bool affine () const { return mapping_->affine(); }
       GeometryType type () const { return mapping_->type(); }
 
-      std::vector<GlobalCoordinate> interpolatoryVertices() const { return mapping_->vertexSet(); }
-      int corners () const { return mapping_->nCorner(); }
+      std::vector<GlobalCoordinate> interpolatoryVertices() const { return mapping_->vertices(); }
+      int corners () const { return mapping_->corners(); }
       GlobalCoordinate corner ( const int i ) const { return mapping_->corner( i ); }
       GlobalCoordinate center () const { return mapping_->center(); }
 
       GlobalCoordinate global ( const LocalCoordinate &local ) const { return mapping_->global( local ); }
       LocalCoordinate local ( const GlobalCoordinate &global ) const
       {
-        LocalCoordinate l;
-        DUNE_UNUSED bool b = mapping_->local( global, l );
-        assert(b);
-        return l;
+        return mapping_->local(global);
       }
 
       GlobalCoordinate normal ( const LocalCoordinate &local ) const { return mapping_->normal( local ); }
 
       ctype integrationElement ( const LocalCoordinate &local ) const { return mapping_->integrationElement( local ); }
-      ctype volume () const { return mapping_->volume(std::numeric_limits<ctype>::epsilon()); }
+      ctype volume () const { return mapping_->volume(); }
 
       JacobianTransposed jacobianTransposed ( const LocalCoordinate &local ) const { return mapping_->jacobianTransposed( local ); }
       JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate &local ) const { return mapping_->jacobianInverseTransposed( local ); }

dune/curvedsurfacegrid/test/test_curvedsurfacegrid.cc

@@ -16,7 +16,7 @@
 #define STR(s) STR_HELPER(s)
 #define STR_HELPER(s) #s
 
-const int order = 3;
+const int order = 2;
 const int quad_order = order+5;
 const int num_levels = 5;
 

example/example.cc

@@ -37,10 +37,10 @@ int main(int argc, char** argv)
   //load host-grid (ALUGrid)
   using HostGridType = Dune::ALUGrid<2, 3, simplex, conforming>;
   GridFactory<HostGridType> gridFactory;
-  GmshReader<HostGridType>::read(gridFactory, "dune-curvedsurfacegrid/example/sphere_coarse.msh");
+  GmshReader<HostGridType>::read(gridFactory, "example/sphere_coarse.msh");
   std::shared_ptr<HostGridType> hostGrid(gridFactory.createGrid());
 
-  const int ip_order = 5; //interpolatory order for curved elements
+  const int ip_order = 3; //interpolatory order for curved elements
 
   //case1: projection to actual sphere
   {
@@ -58,7 +58,7 @@ int main(int argc, char** argv)
 
   //case2: projection to a surface-mesh of a sphere
   {
-    SurfaceDistanceCoordFunction surfdistCF("dune-curvedsurfacegrid/example/sphere_fine.vtu");
+    SurfaceDistanceCoordFunction surfdistCF("example/sphere_fine.vtu");
     using GridType = CurvedSurfaceGrid<HostGridType, SurfaceDistanceCoordFunction, ip_order>;
     std::shared_ptr<GridType> grid(new GridType(*hostGrid, surfdistCF));