added example and restructured geometry implementation

dune.module

@@ -9,4 +9,4 @@ Maintainer: simon.praetorius@tu-dresden.de
 # Required build dependencies
 Depends: dune-geometry dune-localfunctions
 # Optional build dependencies
-#Suggests:
+Suggests: dune-grid

dune/curvedgeometry/curvedgeometry.hh

@@ -131,19 +131,17 @@ namespace Dune
    *  \tparam  ct      coordinate type
    *  \tparam  mydim   geometry dimension
    *  \tparam  cdim    coordinate dimension
-   *  \tparam  Traits  traits allowing to tweak some implementation details
-   *                   (optional)
+   *  \tparam  polynomial_order
    *
    *  The requirements on the traits are documented along with their default,
    *  CurvedGeometryTraits.
    */
-  template <int mydim, int cdim, class GridImpl>
+  template <class ct, int mydim, int cdim, int polynomial_order>
   class CurvedGeometry
-      : public GeometryDefaultImplementation<mydim, cdim, GridImpl, CurvedGeometry>
   {
   public:
     /// coordinate type
-    using ctype = typename GridImpl::ctype;
+    using ctype = ct;
 
     using Traits = CurvedGeometryTraits<ctype>;
 
@@ -153,7 +151,7 @@ namespace Dune
     static const int coorddimension = cdim;
 
     /// Polynomial order of geometry
-    static const int order = GridImpl::order;
+    static const int order = polynomial_order;
 
     /// type of local coordinates
     using LocalCoordinate = FieldVector<ctype, mydimension>;

src/dune-curvedgeometry.cc

@@ -4,27 +4,55 @@
 #ifdef HAVE_CONFIG_H
 # include "config.h"
 #endif
+
+#include <cmath>
 #include <iostream>
-#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
-#include <dune/common/exceptions.hh> // We use exceptions
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/curvedgeometry/curvedgeometry.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
+
+const int order = 3;
 
 int main(int argc, char** argv)
 {
-  try{
-    // Maybe initialize MPI
-    Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
-    std::cout << "Hello World! This is dune-curvedgeometry." << std::endl;
-    if(Dune::MPIHelper::isFake)
-      std::cout<< "This is a sequential program." << std::endl;
-    else
-      std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
-        <<" processes!"<<std::endl;
-    return 0;
-  }
-  catch (Dune::Exception &e){
-    std::cerr << "Dune reported error: " << e << std::endl;
-  }
-  catch (...){
-    std::cerr << "Unknown exception thrown!" << std::endl;
+  using namespace Dune;
+  MPIHelper& helper = MPIHelper::instance(argc, argv);
+
+  YaspGrid<2> grid({1.0,1.0}, {4,4});
+  using LocalCoordinate  = FieldVector<double,2>;
+  using GlobalCoordinate = FieldVector<double,2>;
+  using WorldCoordinate  = FieldVector<double,3>;
+
+  using LocalFECache = LagrangeLocalFiniteElementCache<double, double, 2, order>;
+  LocalFECache localFeCache;
+
+  // coordinate projection
+  auto project = [](GlobalCoordinate const& global) -> WorldCoordinate
+  {
+    return { global[0],
+             global[1],
+             std::sin(global[0])*std::cos(global[1]) };
+  };
+
+  for (auto const& e : elements(grid.leafGridView()))
+  {
+    // projection from local coordinates
+    auto X = [project,geo=e.geometry()](LocalCoordinate const& local) -> WorldCoordinate { return project(geo.global(local)); };
+
+    // construct lagrange vertices
+    std::vector<WorldCoordinate> vertices;
+    localFeCache.get(e.type()).localInterpolation().interpolate(X, vertices);
+
+    // create a curved geometry
+    CurvedGeometry<double,2,3,order> geometry(e.type(), vertices);
+
+    auto const& quadRule = QuadratureRules<double,2>::rule(e.type(), 4);
+    for (auto const& qp : quadRule) {
+      auto Jt = geometry.jacobianTransposed(qp.position());
+      auto Jtinv = geometry.jacobianInverseTransposed(qp.position());
+    }
   }
 }