added datacollector for curved vtk-writer output

dune.module

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

dune/curvedgeometry/curvedgeometry.hh

@@ -53,13 +53,14 @@ namespace Dune
 
   /// \brief default traits class for CurvedGeometry
   /**
-   *  The CurvedGeometry (and CachedCurvedGeometry) allow tweaking
+   *  The CurvedGeometry allow tweaking
    *  some implementation details through a traits class.
    *
    *  This structure provides the default values.
    *
    *  \tparam  ct  coordinate type
-   *  \tparam  LFECache  A LocalFiniteElementVariantCache implementation, e.g. LagrangeLocalFiniteElementCache
+   *  \tparam  LFECache  A LocalFiniteElementVariantCache implementation, 
+   *                     e.g. LagrangeLocalFiniteElementCache
    */
   template <class ct, class LFECache>
   struct CurvedGeometryTraits
@@ -87,12 +88,12 @@ namespace Dune
    *  \tparam  ct      coordinate type
    *  \tparam  mydim   geometry dimension
    *  \tparam  cdim    coordinate dimension
-   *  \tparam  Traits  Parametrization of the geometry, see \ref CurvedGeometryTraits
+   *  \tparam  TraitsType  Parametrization of the geometry, see \ref CurvedGeometryTraits
    *
    *  The requirements on the traits are documented along with their default,
    *  CurvedGeometryTraits.
    */
-  template <class ct, int mydim, int cdim, class Traits>
+  template <class ct, int mydim, int cdim, class TraitsType>
   class CurvedGeometry
   {
   public:
@@ -124,6 +125,9 @@ namespace Dune
     /// type of reference element
     using ReferenceElement = typename ReferenceElements::ReferenceElement;
 
+    /// Parametrization of the geometry
+    using Traits = TraitsType;
+
   protected:
     using MatrixHelper = typename Traits::MatrixHelper;
 
@@ -141,7 +145,7 @@ namespace Dune
 
 
   public:
-    /// \brief constructor
+    /// \brief constructor from a vector of coefficients of the LocalBasis parametrizing the geometry
     /**
      *  \param[in]  refElement  reference element for the geometry
      *  \param[in]  vertices    vertices to store internally
@@ -157,6 +161,7 @@ namespace Dune
       assert(localFE_.size() == vertices_.size());
     }
 
+    /// \brief constructor from a local parametrization function, mapping local to (curved) global coordinates
     template <class Parametrization,
       std::enable_if_t<Std::is_callable<Parametrization(LocalCoordinate), GlobalCoordinate>::value, bool> = true>
     CurvedGeometry (const ReferenceElement& refElement, Parametrization&& param)
@@ -166,7 +171,7 @@ namespace Dune
       localInterpolation.interpolate(param, vertices_);
     }
 
-    /// \brief constructor
+    /// \brief constructor, forwarding to the other constructor that take a reference-element
     /**
      *  \param[in]  gt     geometry type
      *  \param[in]  param  either a vector of vertices, or a functor that can be used to construct the vertices
@@ -379,8 +384,8 @@ namespace Dune
       return out;
     }
 
-    /** \brief obtain the transposed of the Jacobian's inverse
-     *
+    /// \brief obtain the transposed of the Jacobian's inverse
+    /**
      *  The Jacobian's inverse is defined as a pseudo-inverse. If we denote
      *  the Jacobian by \f$J(x)\f$, the following condition holds:
      *  \f[J^{-1}(x) J(x) = I.\f]
@@ -392,12 +397,14 @@ namespace Dune
       return out;
     }
 
+    /// \brief obtain the reference-element related to this geometry
     friend ReferenceElement referenceElement (const CurvedGeometry& geometry)
     {
       return geometry.refElement();
     }
 
-    const std::vector<GlobalCoordinate>& vertices () const
+    /// \brief obtain the coefficients of the parametrization
+    const std::vector<GlobalCoordinate>& coefficients () const
     {
       return vertices_;
     }
@@ -440,6 +447,7 @@ namespace Dune
     std::vector<GlobalCoordinate> vertices_;
   };
 
+  /// \brief curved geometry parametrized with lagrange basis functions
   template <class ctype, int mydim, int cdim, int order>
   using LagrangeCurvedGeometry = CurvedGeometry<ctype,mydim,cdim,
     CurvedGeometryTraits<ctype, LagrangeLocalFiniteElementCache<ctype, ctype, mydim, order>> >;

dune/curvedgeometry/curvedgeometrydatacollector.hh

+#ifndef DUNE_CURVEDGEOMETRY_DATACOLLECTOR_HH
+#define DUNE_CURVEDGEOMETRY_DATACOLLECTOR_HH
+
+#if HAVE_DUNE_VTK
+
+#include <functional>
+#include <vector>
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/grid/common/partitionset.hh>
+#include <dune/vtk/forward.hh>
+#include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/datacollectors/quadraticdatacollector.hh>
+
+namespace Dune
+{
+  // Extract point coordinates from parametrized geometry
+  // and write nonlinear vtk elements
+  template <class GridView, class Geometry>
+  class CurvedGeometryDataCollector
+      : public UnstructuredDataCollectorInterface<GridView, CurvedGeometryDataCollector<GridView,Geometry>, Partitions::All>
+  {
+    using Self = CurvedGeometryDataCollector;
+    using Super = UnstructuredDataCollectorInterface<GridView, CurvedGeometryDataCollector<GridView,Geometry>, Partitions::All>;
+
+    using LocalCoord = typename GridView::template Codim<0>::Entity::Geometry::LocalCoordinate;
+    using GlobalCoord = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
+
+    using Mapping = std::function<typename Geometry::GlobalCoordinate(GlobalCoord)>;
+    using Parametrization = std::vector<std::vector<typename Geometry::GlobalCoordinate>>;
+
+  public:
+    using Super::dim;
+    using Super::partition;
+
+  public:
+    // Construct the curved geometry from a coordinate mapping
+    CurvedGeometryDataCollector (GridView const& gridView, Mapping mapping)
+      : Super(gridView)
+      , dataCollector_(gridView)
+      , mapping_(mapping)
+    {}
+
+    // construct the curved geometry from a stored parametrization, i.e. a vector of coefficient vectors
+    CurvedGeometryDataCollector (GridView const& gridView, Parametrization const& parametrization)
+      : Super(gridView)
+      , dataCollector_(gridView)
+      , parametrization_(&parametrization)
+    {}
+
+    // collect the points of the geometry. This uses the curved geometry instead of the element geometry
+    template <class T>
+    std::vector<T> pointsImpl () const
+    {
+      std::vector<T> data(Super::numPoints() * 3);
+      auto const& indexSet = gridView_.indexSet();
+      for (auto const& element : elements(gridView_, partition)) {
+        auto refElem = referenceElement<T,Geometry::mydimension>(element.type());
+        auto geometry = makeGeometry(indexSet, element);
+
+        // vertices
+        for (unsigned int i = 0; i < element.subEntities(dim); ++i) {
+          std::size_t idx = 3 * indexSet.subIndex(element, i, dim);
+          auto v = geometry.global(refElem.position(i,dim));
+          for (std::size_t j = 0; j < v.size(); ++j)
+            data[idx + j] = T(v[j]);
+          for (std::size_t j = v.size(); j < 3u; ++j)
+            data[idx + j] = T(0);
+        }
+        // edge centers
+        for (unsigned int i = 0; i < element.subEntities(dim-1); ++i) {
+          std::size_t idx = 3 * (indexSet.subIndex(element, i, dim-1) + gridView_.size(dim));
+          auto v = geometry.global(refElem.position(i,dim-1));
+          for (std::size_t j = 0; j < v.size(); ++j)
+            data[idx + j] = T(v[j]);
+          for (std::size_t j = v.size(); j < 3u; ++j)
+            data[idx + j] = T(0);
+        }
+      }
+      return data;
+    }
+
+  private:
+
+    // construct the geometry object, either from the mapping or the coefficient vector
+    template <class IndexSet, class Element>
+    Geometry makeGeometry(IndexSet const& indexSet, Element const& element) const
+    {
+      if (parametrization_ != nullptr) {
+        return Geometry{element.type(), (*parametrization_)[indexSet.index(element)]};
+      } else {
+        // mapping from local to global coordinates
+        auto X = [mapping=mapping_,geo=element.geometry()](LocalCoord const& local)
+          -> typename Geometry::GlobalCoordinate
+        {
+          return mapping(geo.global(local));
+        };
+
+        // create a curved geometry
+        return Geometry{element.type(), X};
+      }
+    }
+
+  public: // methods forwarded to quadratic data collector
+
+    void updateImpl ()
+    {
+      dataCollector_.update();
+    }
+
+    std::uint64_t numPointsImpl () const
+    {
+      return dataCollector_.numPoints();
+    }
+
+    std::vector<std::uint64_t> pointIdsImpl () const
+    {
+      return dataCollector_.pointIds();
+    }
+
+    std::uint64_t numCellsImpl () const
+    {
+      return dataCollector_.numCells();
+    }
+
+    Cells cellsImpl () const
+    {
+      return dataCollector_.cells();
+    }
+
+    template <class T, class GlobalFunction>
+    std::vector<T> pointDataImpl (GlobalFunction const& fct) const
+    {
+      return dataCollector_.template pointData<T>(fct);
+    }
+
+  private:
+    using Super::gridView_;
+    QuadraticDataCollector<GridView> dataCollector_;
+    Mapping mapping_;
+    Parametrization const* parametrization_ = nullptr;
+  };
+
+} // end namespace Dune
+
+#endif // HAVE_DUNE_VTK
+
+#endif // DUNE_CURVEDGEOMETRY_DATACOLLECTOR_HH

src/CMakeLists.txt

-add_executable("dune-curvedgeometry" dune-curvedgeometry.cc)
-target_link_dune_default_libraries("dune-curvedgeometry")
+add_executable("curvedgeometry" curvedgeometry.cc)
+target_link_dune_default_libraries("curvedgeometry")
+
+add_executable("parametrization" parametrization.cc)
+target_link_dune_default_libraries("parametrization")
+
+add_executable("lagrangesimplex" lagrangesimplex.cc)
+target_link_dune_default_libraries("lagrangesimplex")
+
+if (dune-vtk_FOUND)
+  add_executable("boundaryprojection" boundaryprojection.cc)
+  target_link_dune_default_libraries("boundaryprojection")
+
+  add_executable("vtkwriter" vtkwriter.cc)
+  target_link_dune_default_libraries("vtkwriter")
+endif (dune-vtk_FOUND)

src/boundaryprojection.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <cmath>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/curvedgeometry/curvedgeometry.hh>
+#include <dune/curvedgeometry/curvedgeometrydatacollector.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/vtk/vtkwriter.hh>
+
+const int order = 3;
+
+using namespace Dune;
+using LocalCoordinate  = FieldVector<double,2>;
+using GlobalCoordinate = FieldVector<double,2>;
+using WorldCoordinate  = FieldVector<double,2>;
+
+int main(int argc, char** argv)
+{
+  MPIHelper& helper = MPIHelper::instance(argc, argv);
+
+  using Grid = YaspGrid<2,EquidistantOffsetCoordinates<double,2>>;
+  Grid grid({-1.0,-1.0}, {1.0,1.0}, {20,20});
+  using GridView = typename Grid::LeafGridView;
+  GridView gridView = grid.leafGridView();
+
+  // coordinate projection
+  auto project = [](GlobalCoordinate const& global) -> WorldCoordinate
+  {
+    return global * (std::sqrt(2.0)/global.two_norm());
+  };
+
+  std::vector<std::vector<WorldCoordinate>> parametrization(gridView.size(0));
+
+  using Geometry = LagrangeCurvedGeometry<double, 2, 2, order>;
+  using LocalFECache = typename Geometry::Traits::LocalFiniteElementCache;
+  LocalFECache localFECache;
+  auto const& indexSet = gridView.indexSet();
+  for (auto const& e : elements(grid.leafGridView()))
+  {
+    // projection from local coordinates
+    auto id = [geo=e.geometry()](LocalCoordinate const& local) -> WorldCoordinate
+    {
+      return geo.global(local);
+    };
+
+    auto const& localFE = localFECache.get(e.type());
+    auto re = Dune::referenceElement<double,2>(e.type());
+
+    auto& localParametrization = parametrization[indexSet.index(e)];
+    localFE.localInterpolation().interpolate(id, localParametrization);
+    if (e.hasBoundaryIntersections())
+    {
+      std::vector<std::size_t> containedDOFs;
+
+      auto const& localCoeff = localFE.localCoefficients();
+      std::size_t localSize = localCoeff.size();
+      for (auto const& is : intersections(gridView, e)) {
+        if (!is.boundary())
+          continue;
+
+        // collect all DOFs on the boundary
+        std::size_t subEntityIndex = is.indexInInside();
+        std::size_t subEntityCodim = 1;
+        for (std::size_t i = 0; i < localSize; ++i)
+        {
+          auto localKey = localCoeff.localKey(i);
+          if (re.subEntities(subEntityIndex, subEntityCodim, localKey.codim()).contains(localKey.subEntity()))
+          {
+            containedDOFs.push_back(i);
+          }
+        }
+      }
+
+      // project only boundary DOFs
+      for (std::size_t i : containedDOFs)
+        localParametrization[i] = project(localParametrization[i]);
+    }
+  }
+
+  using DataCollector = CurvedGeometryDataCollector<GridView,Geometry>;
+  DataCollector dataCollector(gridView, parametrization);
+
+  using Writer = VtkUnstructuredGridWriter<GridView,DataCollector>;
+  Writer writer(stackobject_to_shared_ptr(dataCollector));
+
+  writer.write("boundaryprojection.vtu");
+}

src/dune-curvedgeometry.cc → src/curvedgeometry.cc

@@ -13,7 +13,6 @@
 #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;
 
@@ -44,7 +43,10 @@ int main(int argc, char** argv)
   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)); };
+    auto X = [project,geo=e.geometry()](LocalCoordinate const& local) -> WorldCoordinate
+    {
+      return project(geo.global(local));
+    };
 
     // create a curved geometry
     Geometry geometry(e.type(), X);
@@ -57,7 +59,8 @@ int main(int argc, char** argv)
       FieldMatrix<double, 2, 2> res;
       FMatrixHelp::multMatrix(Jt, Jtinv, res);
       res -= I;
-      test.check(res.frobenius_norm() < std::sqrt(std::numeric_limits<double>::epsilon()), "J^-1 * J == I");
+      test.check(res.frobenius_norm() < std::sqrt(std::numeric_limits<double>::epsilon()),
+                 "J^-1 * J == I");
     }
   }
 

src/lagrangesimplex.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+#include <dune/common/rangeutilities.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/fvector.hh>
+#include <dune/geometry/type.hh>
+#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
+#include <dune/localfunctions/lagrange/equidistantpoints.hh>
+#include <dune/localfunctions/lagrange/interpolation.hh>
+
+using namespace Dune;
+
+template <int dim, int order>
+void printBasis(LagrangeSimplexLocalFiniteElement<double, double, dim, order> const& localFE)
+{
+  auto id = [](auto const& local) { return local; };
+
+  std::vector<FieldVector<double,dim>> lagrangePoints;
+  localFE.localInterpolation().interpolate(id, lagrangePoints);
+
+  for (std::size_t i = 0; i < lagrangePoints.size(); ++i) {
+    std::cout << i << ") " << lagrangePoints[i] << std::endl;
+  }
+}
+
+template <unsigned int dim>
+void printBasis2(LocalLagrangeInterpolation<EquidistantPointSet, dim, double> const& localInterpolation)
+{
+  auto id = [](auto const& local) { return local; };
+
+  std::vector<FieldVector<double,dim>> lagrangePoints;
+  localInterpolation.interpolate(id, lagrangePoints);
+
+  for (std::size_t i = 0; i < lagrangePoints.size(); ++i) {
+    std::cout << "\\coordinate (V" << i << ") at (" << lagrangePoints[i][0] << "," << lagrangePoints[i][1] << ");" << std::endl;
+  }
+}
+
+int main(int argc, char** argv)
+{
+  auto range = StaticIntegralRange<std::size_t, 6, 1>{};
+  Hybrid::forEach(range, [](auto k) {
+    const int dim = 2;
+    LagrangeSimplexLocalFiniteElement<double, double, dim, decltype(k)::value> fe;
+    std::cout << "LagrangeSimplexLocalFiniteElement<dim=" << dim << ", k=" << k << ">:" << std::endl;
+    printBasis(fe);
+  });
+
+
+  for (int k = 1; k < 6; ++k) {
+    const int dim = 2;
+    using Factory = LagrangeInterpolationFactory<EquidistantPointSet, dim, double>;
+    using Interpolation = LocalLagrangeInterpolation<EquidistantPointSet, dim, double>;
+    const Interpolation* interpolation = Factory::create<typename Impl::SimplexTopology<dim>::type>(k);
+    std::cout << "LocalLagrangeInterpolation<dim=" << dim << ", k=" << k << ">:" << std::endl;
+    printBasis2(*interpolation);
+  }
+}

src/parametrization.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <cmath>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/curvedgeometry/curvedgeometry.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/grid/yaspgrid.hh>
+
+const int order = 3;
+
+using namespace Dune;
+using LocalCoordinate  = FieldVector<double,2>;
+using GlobalCoordinate = FieldVector<double,2>;
+using WorldCoordinate  = FieldVector<double,3>;
+
+int main(int argc, char** argv)
+{
+  MPIHelper& helper = MPIHelper::instance(argc, argv);
+
+  YaspGrid<2> grid({8.0,8.0}, {32,32});
+  auto gridView = grid.leafGridView();
+  auto const& indexSet = gridView.indexSet();
+
+  FieldMatrix<double,2,2> I{{1,0},{0,1}};
+
+  // coordinate projection
+  auto project = [](GlobalCoordinate const& global) -> WorldCoordinate
+  {
+    return { global[0],
+             global[1],
+             std::sin(global[0])*std::cos(global[1]) };
+  };
+
+  TestSuite test("curved geometry");
+
+  // 1. fill a global parametrization storage
+  using Geometry = LagrangeCurvedGeometry<double, 2, 3, order>;
+  using LocalFECache = typename Geometry::Traits::LocalFiniteElementCache;
+  LocalFECache localFECache;
+
+  std::vector<std::vector<WorldCoordinate>> parametrization(gridView.size(0));
+  for (auto const& e : elements(gridView))
+  {
+    // projection from local coordinates
+    auto X = [project,geo=e.geometry()](LocalCoordinate const& local) -> WorldCoordinate
+    {
+      return project(geo.global(local));
+    };
+
+    auto& localParametrization = parametrization[indexSet.index(e)];
+    auto const& localFE = localFECache.get(e.type());
+
+    // store coefficients of local parametrization in global storage
+    localFE.localInterpolation().interpolate(X, localParametrization);
+  }
+
+
+  // 2. use the stored parametrization to construct a geometry
+  for (auto const& e : elements(grid.leafGridView()))
+  {
+    // create a curved geometry
+    Geometry geometry(e.type(), parametrization[indexSet.index(e)]);
+
+    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());
+
+      FieldMatrix<double, 2, 2> res;
+      FMatrixHelp::multMatrix(Jt, Jtinv, res);
+      res -= I;
+      test.check(res.frobenius_norm() < std::sqrt(std::numeric_limits<double>::epsilon()),
+                 "J^-1 * J == I");
+    }
+  }
+
+  return test.report() ? 0 : 1;
+}

src/vtkwriter.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <cmath>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/curvedgeometry/curvedgeometry.hh>
+#include <dune/curvedgeometry/curvedgeometrydatacollector.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/vtk/vtkwriter.hh>
+
+const int order = 3;
+
+int main(int argc, char** argv)
+{
+  using namespace Dune;
+  MPIHelper& helper = MPIHelper::instance(argc, argv);
+
+  using Grid = YaspGrid<2>;
+  Grid grid({8.0,8.0}, {32,32});
+  using GridView = typename Grid::LeafGridView;
+  GridView gridView = grid.leafGridView();
+
+  using Geometry = LagrangeCurvedGeometry<double, 2, 3, order>;
+
+  // coordinate projection
+  using GlobalCoordinate = typename Grid::Codim<0>::Entity::Geometry::GlobalCoordinate;
+  using WorldCoordinate  = typename Geometry::GlobalCoordinate;
+  auto project = [](GlobalCoordinate const& global) -> WorldCoordinate
+  {
+    return { global[0],
+             global[1],
+             std::sin(global[0])*std::cos(global[1]) };
+  };
+
+  using DataCollector = CurvedGeometryDataCollector<GridView,Geometry>;
+  DataCollector dataCollector(gridView, project);
+
+  using Writer = VtkUnstructuredGridWriter<GridView,DataCollector>;
+  Writer writer(stackobject_to_shared_ptr(dataCollector));
+
+  writer.write("curvedgeometry.vtu");
+}