From e7fdb0691070ac1184c85bfdd121aaaa7f7e612f Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Tue, 18 Feb 2020 23:48:08 +0100
Subject: [PATCH] added datacollector for curved vtk-writer output
---
dune.module | 2 +-
dune/curvedgeometry/curvedgeometry.hh | 26 +--
.../curvedgeometrydatacollector.hh | 148 ++++++++++++++++++
src/CMakeLists.txt | 18 ++-
src/boundaryprojection.cc | 96 ++++++++++++
...ne-curvedgeometry.cc => curvedgeometry.cc} | 9 +-
src/lagrangesimplex.cc | 67 ++++++++
src/parametrization.cc | 86 ++++++++++
src/vtkwriter.cc | 50 ++++++
9 files changed, 487 insertions(+), 15 deletions(-)
create mode 100644 dune/curvedgeometry/curvedgeometrydatacollector.hh
create mode 100644 src/boundaryprojection.cc
rename src/{dune-curvedgeometry.cc => curvedgeometry.cc} (91%)
create mode 100644 src/lagrangesimplex.cc
create mode 100644 src/parametrization.cc
create mode 100644 src/vtkwriter.cc
diff --git a/dune.module b/dune.module
index 008d78a..35f2124 100644
--- a/dune.module
+++ b/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
diff --git a/dune/curvedgeometry/curvedgeometry.hh b/dune/curvedgeometry/curvedgeometry.hh
index 1aa7c55..ce397fa 100644
--- a/dune/curvedgeometry/curvedgeometry.hh
+++ b/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>> >;
diff --git a/dune/curvedgeometry/curvedgeometrydatacollector.hh b/dune/curvedgeometry/curvedgeometrydatacollector.hh
new file mode 100644
index 0000000..67b7c4a
--- /dev/null
+++ b/dune/curvedgeometry/curvedgeometrydatacollector.hh
@@ -0,0 +1,148 @@
+#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_(¶metrization)
+ {}
+
+ // 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
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 1e26eeb..73c6ddd 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,2 +1,16 @@
-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)
diff --git a/src/boundaryprojection.cc b/src/boundaryprojection.cc
new file mode 100644
index 0000000..6173dc8
--- /dev/null
+++ b/src/boundaryprojection.cc
@@ -0,0 +1,96 @@
+// -*- 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");
+}
diff --git a/src/dune-curvedgeometry.cc b/src/curvedgeometry.cc
similarity index 91%
rename from src/dune-curvedgeometry.cc
rename to src/curvedgeometry.cc
index 8a7e697..cf5df5a 100644
--- a/src/dune-curvedgeometry.cc
+++ b/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");
}
}
diff --git a/src/lagrangesimplex.cc b/src/lagrangesimplex.cc
new file mode 100644
index 0000000..82bb4d0
--- /dev/null
+++ b/src/lagrangesimplex.cc
@@ -0,0 +1,67 @@
+// -*- 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);
+ }
+}
diff --git a/src/parametrization.cc b/src/parametrization.cc
new file mode 100644
index 0000000..104179d
--- /dev/null
+++ b/src/parametrization.cc
@@ -0,0 +1,86 @@
+// -*- 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;
+}
diff --git a/src/vtkwriter.cc b/src/vtkwriter.cc
new file mode 100644
index 0000000..7a8d14d
--- /dev/null
+++ b/src/vtkwriter.cc
@@ -0,0 +1,50 @@
+// -*- 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");
+}
--
GitLab