From f88922d99bbd2307b0c3104dd26f22b6fb3420b1 Mon Sep 17 00:00:00 2001
From: "Praetorius, Simon" <simon.praetorius@tu-dresden.de>
Date: Thu, 16 Apr 2020 12:48:41 +0200
Subject: [PATCH] higher order vtk reader
---
doc/triangles_2d_order3.vtu | 103 +++
doc/triangles_3d_order3.vtu | 140 ++++
dune.module | 2 +-
dune/vtk/datacollectors/CMakeLists.txt | 2 +
.../datacollectors/lagrangedatacollector.hh | 196 ++++++
dune/vtk/datacollectors/test/CMakeLists.txt | 2 +
.../test/test-lagrangedatacollector.cc | 119 ++++
dune/vtk/forward.hh | 3 +
dune/vtk/gridcreatorinterface.hh | 16 +-
dune/vtk/gridcreators/CMakeLists.txt | 1 +
.../vtk/gridcreators/continuousgridcreator.hh | 3 +-
.../gridcreators/discontinuousgridcreator.hh | 3 +-
dune/vtk/gridcreators/lagrangegridcreator.hh | 344 ++++++++++
dune/vtk/utility/CMakeLists.txt | 4 +
dune/vtk/utility/lagrangelfecache.hh | 35 +
dune/vtk/utility/lagrangepoints.hh | 605 ++++++++++++++++++
dune/vtk/utility/test/CMakeLists.txt | 13 +
dune/vtk/utility/test/test-lagrange.cc | 169 +++++
dune/vtk/utility/test/test-lfecache.cc | 56 ++
dune/vtk/vtkfunction.hh | 10 +-
dune/vtk/vtkreader.hh | 11 +-
dune/vtk/vtkreader.impl.hh | 22 +-
dune/vtk/vtktypes.cc | 41 ++
dune/vtk/vtktypes.hh | 15 +-
dune/vtk/vtkwriterinterface.hh | 8 +-
src/CMakeLists.txt | 11 +-
src/datacollector.cc | 27 +-
src/lagrangepoints.cc | 68 ++
src/lagrangereader.cc | 95 +++
29 files changed, 2097 insertions(+), 27 deletions(-)
create mode 100644 doc/triangles_2d_order3.vtu
create mode 100644 doc/triangles_3d_order3.vtu
create mode 100644 dune/vtk/datacollectors/lagrangedatacollector.hh
create mode 100644 dune/vtk/datacollectors/test/CMakeLists.txt
create mode 100644 dune/vtk/datacollectors/test/test-lagrangedatacollector.cc
create mode 100644 dune/vtk/gridcreators/lagrangegridcreator.hh
create mode 100644 dune/vtk/utility/lagrangelfecache.hh
create mode 100644 dune/vtk/utility/lagrangepoints.hh
create mode 100644 dune/vtk/utility/test/CMakeLists.txt
create mode 100644 dune/vtk/utility/test/test-lagrange.cc
create mode 100644 dune/vtk/utility/test/test-lfecache.cc
create mode 100644 src/lagrangepoints.cc
create mode 100644 src/lagrangereader.cc
diff --git a/doc/triangles_2d_order3.vtu b/doc/triangles_2d_order3.vtu
new file mode 100644
index 0000000..ff528b4
--- /dev/null
+++ b/doc/triangles_2d_order3.vtu
@@ -0,0 +1,103 @@
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+ <UnstructuredGrid>
+ <Piece NumberOfPoints="49" NumberOfCells="8">
+ <PointData>
+ <DataArray type="Float32" Name="DistanceToCenter" format="ascii" RangeMin="0" RangeMax="1.4142135381698608">
+ 1.4142135381698608 1 1.4142135381698608 1 0 1
+ 1.4142135381698608 1 1.4142135381698608 1.201850414276123 1.054092526435852 1.201850414276123
+ 0.9428090453147888 0.7453559637069702 0.6666666269302368 1.054092526435852 0.7453559637069702 0.471404492855072
+ 0.3333333134651184 0.6666666269302368 0.3333333134651184 1.054092526435852 1.201850414276123 0.745356023311615
+ 0.942808985710144 1.201850414276123 0.4714045226573944 0.7453559637069702 1.054092526435852 0.3333333730697632
+ 0.6666666269302368 1.054092526435852 0.745356023311615 0.4714045226573944 0.3333333730697632 1.201850414276123
+ 0.942808985710144 0.7453559637069702 0.6666666269302368 1.201850414276123 1.054092526435852 0.4714045822620392
+ 0.7453559637069702 1.054092526435852 0.7453559637069702 0.942808985710144 1.201850414276123 1.054092526435852
+ 1.201850414276123
+ </DataArray>
+ <DataArray type="Float32" Name="Polynomial" format="ascii" RangeMin="1" RangeMax="49">
+ 1 4 15 4 10 26
+ 15 26 49 1 1 1
+ 1 1 4 1 1 1
+ 4 4 4 5 8 5
+ 8 16 6 10 18 11
+ 15 5 5 6 11 8
+ 8 10 15 16 18 14
+ 19 29 19 25 36 29
+ 36
+ </DataArray>
+ </PointData>
+ <CellData>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="0" RangeMax="2.8284271247461903">
+ 0 0 0 1 0 0
+ 2 0 0 0 1 0
+ 1 1 0 2 1 0
+ 0 2 0 1 2 0
+ 2 2 0 0.3333333432674408 0 0
+ 0.6666666865348816 0 0 0 0.3333333432674408 0
+ 0.3333333432674408 0.3333333432674408 0 0.6666666865348816 0.3333333432674408 0
+ 1 0.3333333432674408 0 0 0.6666666865348816 0
+ 0.3333333432674408 0.6666666865348816 0 0.6666666865348816 0.6666666865348816 0
+ 1 0.6666666865348816 0 0.3333333432674408 1 0
+ 0.6666666865348816 1 0 1.3333333730697632 0 0
+ 1.6666666269302368 0 0 1.3333333730697632 0.3333333432674408 0
+ 1.6666666269302368 0.3333333432674408 0 2 0.3333333432674408 0
+ 1.3333333730697632 0.6666666865348816 0 1.6666666269302368 0.6666666865348816 0
+ 2 0.6666666865348816 0 1.3333333730697632 1 0
+ 1.6666666269302368 1 0 0 1.3333333730697632 0
+ 0.3333333432674408 1.3333333730697632 0 0.6666666865348816 1.3333333730697632 0
+ 1 1.3333333730697632 0 0 1.6666666269302368 0
+ 0.3333333432674408 1.6666666269302368 0 0.6666666865348816 1.6666666269302368 0
+ 1 1.6666666269302368 0 0.3333333432674408 2 0
+ 0.6666666865348816 2 0 1.3333333730697632 1.3333333730697632 0
+ 1.6666666269302368 1.3333333730697632 0 2 1.3333333730697632 0
+ 1.3333333730697632 1.6666666269302368 0 1.6666666269302368 1.6666666269302368 0
+ 2 1.6666666269302368 0 1.3333333730697632 2 0
+ 1.6666666269302368 2 0
+ <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
+ <Value index="0">
+ 0
+ </Value>
+ <Value index="1">
+ 2.8284271247
+ </Value>
+ </InformationKey>
+ <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
+ <Value index="0">
+ 0
+ </Value>
+ <Value index="1">
+ 2.8284271247
+ </Value>
+ </InformationKey>
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="48">
+ 0 1 3 9 10 13
+ 16 15 11 12 4 3
+ 1 20 19 16 13 14
+ 18 17 1 2 4 21
+ 22 24 26 18 14 23
+ 5 4 2 30 29 26
+ 24 25 28 27 3 4
+ 6 19 20 33 36 35
+ 31 32 7 6 4 40
+ 39 36 33 34 38 37
+ 4 5 7 29 30 42
+ 44 38 34 41 8 7
+ 5 48 47 44 42 43
+ 46 45
+ </DataArray>
+ <DataArray type="Int64" Name="offsets" format="ascii" RangeMin="10" RangeMax="80">
+ 10 20 30 40 50 60
+ 70 80
+ </DataArray>
+ <DataArray type="UInt8" Name="types" format="ascii" RangeMin="69" RangeMax="69">
+ 69 69 69 69 69 69
+ 69 69
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
diff --git a/doc/triangles_3d_order3.vtu b/doc/triangles_3d_order3.vtu
new file mode 100644
index 0000000..ec93fc4
--- /dev/null
+++ b/doc/triangles_3d_order3.vtu
@@ -0,0 +1,140 @@
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+ <UnstructuredGrid>
+ <Piece NumberOfPoints="49" NumberOfCells="8">
+ <PointData Vectors="Result">
+ <DataArray type="Float32" Name="DistanceToCenter" format="ascii" RangeMin="0" RangeMax="1.4142135381698608">
+ 1.4142135381698608 1 1.4142135381698608 1 0 1
+ 1.4142135381698608 1 1.4142135381698608 1.201850414276123 1.054092526435852 1.201850414276123
+ 0.9428090453147888 0.7453559637069702 0.6666666269302368 1.054092526435852 0.7453559637069702 0.471404492855072
+ 0.3333333134651184 0.6666666269302368 0.3333333134651184 1.054092526435852 1.201850414276123 0.745356023311615
+ 0.942808985710144 1.201850414276123 0.4714045226573944 0.7453559637069702 1.054092526435852 0.3333333730697632
+ 0.6666666269302368 1.054092526435852 0.745356023311615 0.4714045226573944 0.3333333730697632 1.201850414276123
+ 0.942808985710144 0.7453559637069702 0.6666666269302368 1.201850414276123 1.054092526435852 0.4714045822620392
+ 0.7453559637069702 1.054092526435852 0.7453559637069702 0.942808985710144 1.201850414276123 1.054092526435852
+ 1.201850414276123
+ </DataArray>
+ <DataArray type="Float32" Name="Polynomial" format="ascii" RangeMin="1" RangeMax="49">
+ 1 4 15 4 10 26
+ 15 26 49 1 1 1
+ 1 1 4 1 1 1
+ 4 4 4 5 8 5
+ 8 16 6 10 18 11
+ 15 5 5 6 11 8
+ 8 10 15 16 18 14
+ 19 29 19 25 36 29
+ 36
+ </DataArray>
+ <DataArray type="Float64" Name="Result" NumberOfComponents="3" format="ascii" RangeMin="0.08895213036368632" RangeMax="1.9954079615554763">
+ 0 0 1 0 0 1.8414709848078965
+ 0 0 1.9092974268256817 0 0 0.5403023058681398
+ 0 0 1.3817732906760363 0 0 1.4495997326938215
+ 0 0 -0.4161468365471424 0 0 0.4253241482607541
+ 0 0 0.4931505902785393 0 0 1.3271947061834561
+ 0 0 1.618369818683914 0 0 0.9449569430643503
+ 0 0 1.2721516492478064 0 0 1.5633267617482642
+ 0 0 1.7864279278722468 0 0 0.7858872484910437
+ 0 0 1.1130819546744999 0 0 1.4042570671749577
+ 0 0 1.6273582332989402 0 0 0.8674970120515959
+ 0 0 1.1586721245520537 0 0 1.9719379107108135
+ 0 0 1.9954079615554763 0 0 1.9168948537751636
+ 0 0 1.9403649046198266 0 0 1.854254369890032
+ 0 0 1.757825159201857 0 0 1.78129521004652
+ 0 0 1.6951846753167255 0 0 1.5122402165789532
+ 0 0 1.535710267423616 0 0 0.23523753468164693
+ 0 0 0.562432240865103 0 0 0.8536073533655608
+ 0 0 1.0767085194895434 0 0 -0.09572350846041702
+ 0 0 0.23147119772303906 0 0 0.522646310223497
+ 0 0 0.7457474763474795 0 0 -0.08895213036368632
+ 0 0 0.20222298213677153 0 0 1.2071754453924604
+ 0 0 1.2306454962371232 0 0 1.1445349615073286
+ 0 0 0.8762144022503964 0 0 0.8996844530950593
+ 0 0 0.8135739183652647 0 0 0.5557910741636709
+ 0 0 0.5792611250083339
+ <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
+ <Value index="0">
+ 0.088952130364
+ </Value>
+ <Value index="1">
+ 1.9954079616
+ </Value>
+ </InformationKey>
+ <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
+ <Value index="0">
+ 0.088952130364
+ </Value>
+ <Value index="1">
+ 1.9954079616
+ </Value>
+ </InformationKey>
+ </DataArray>
+ </PointData>
+ <CellData>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="1" RangeMax="2.871096916901671">
+ 0 0 1 1 0 1.8414709568023682
+ 2 0 1.9092974662780762 0 1 0.5403022766113281
+ 1 1 1.3817732334136963 2 1 1.4495997428894043
+ 0 2 -0.416146844625473 1 2 0.42532414197921753
+ 2 2 0.49315059185028076 0.3333333432674408 0 1.3271946907043457
+ 0.6666666865348816 0 1.6183698177337646 0 0.3333333432674408 0.9449569582939148
+ 0.3333333432674408 0.3333333432674408 1.2721517086029053 0.6666666865348816 0.3333333432674408 1.5633267164230347
+ 1 0.3333333432674408 1.7864279747009277 0 0.6666666865348816 0.7858872413635254
+ 0.3333333432674408 0.6666666865348816 1.113081932067871 0.6666666865348816 0.6666666865348816 1.40425705909729
+ 1 0.6666666865348816 1.6273581981658936 0.3333333432674408 1 0.8674970269203186
+ 0.6666666865348816 1 1.1586720943450928 1.3333333730697632 0 1.971937894821167
+ 1.6666666269302368 0 1.9954079389572144 1.3333333730697632 0.3333333432674408 1.9168949127197266
+ 1.6666666269302368 0.3333333432674408 1.940364956855774 2 0.3333333432674408 1.8542543649673462
+ 1.3333333730697632 0.6666666865348816 1.7578251361846924 1.6666666269302368 0.6666666865348816 1.7812951803207397
+ 2 0.6666666865348816 1.6951847076416016 1.3333333730697632 1 1.5122401714324951
+ 1.6666666269302368 1 1.5357102155685425 0 1.3333333730697632 0.23523753881454468
+ 0.3333333432674408 1.3333333730697632 0.5624322295188904 0.6666666865348816 1.3333333730697632 0.8536073565483093
+ 1 1.3333333730697632 1.0767085552215576 0 1.6666666269302368 -0.09572350978851318
+ 0.3333333432674408 1.6666666269302368 0.2314711958169937 0.6666666865348816 1.6666666269302368 0.5226463079452515
+ 1 1.6666666269302368 0.745747447013855 0.3333333432674408 2 -0.08895213156938553
+ 0.6666666865348816 2 0.20222298800945282 1.3333333730697632 1.3333333730697632 1.2071754932403564
+ 1.6666666269302368 1.3333333730697632 1.2306455373764038 2 1.3333333730697632 1.144534945487976
+ 1.3333333730697632 1.6666666269302368 0.8762143850326538 1.6666666269302368 1.6666666269302368 0.8996844291687012
+ 2 1.6666666269302368 0.8135738968849182 1.3333333730697632 2 0.5557910799980164
+ 1.6666666269302368 2 0.5792611241340637
+ <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
+ <Value index="0">
+ 1
+ </Value>
+ <Value index="1">
+ 2.8710969169
+ </Value>
+ </InformationKey>
+ <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
+ <Value index="0">
+ 1
+ </Value>
+ <Value index="1">
+ 2.8710969169
+ </Value>
+ </InformationKey>
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="48">
+ 0 1 3 9 10 13 16 15 11 12
+ 4 3 1 20 19 16 13 14 18 17
+ 1 2 4 21 22 24 26 18 14 23
+ 5 4 2 30 29 26 24 25 28 27
+ 3 4 6 19 20 33 36 35 31 32
+ 7 6 4 40 39 36 33 34 38 37
+ 4 5 7 29 30 42 44 38 34 41
+ 8 7 5 48 47 44 42 43 46 45
+ </DataArray>
+ <DataArray type="Int64" Name="offsets" format="ascii" RangeMin="10" RangeMax="80">
+ 10 20 30 40 50 60
+ 70 80
+ </DataArray>
+ <DataArray type="UInt8" Name="types" format="ascii" RangeMin="69" RangeMax="69">
+ 69 69 69 69 69 69
+ 69 69
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
diff --git a/dune.module b/dune.module
index e318f46..3ac7b47 100644
--- a/dune.module
+++ b/dune.module
@@ -8,4 +8,4 @@ Version: 0.2
Maintainer: simon.praetorius@tu-dresden.de
#depending on
Depends: dune-common (>= 2.6) dune-geometry (>= 2.6) dune-grid dune-uggrid
-Suggests: dune-functions dune-spgrid dune-polygongrid dune-alugrid
+Suggests: dune-functions dune-spgrid dune-polygongrid dune-alugrid dune-foamgrid
diff --git a/dune/vtk/datacollectors/CMakeLists.txt b/dune/vtk/datacollectors/CMakeLists.txt
index 3ee908c..0b6e52a 100644
--- a/dune/vtk/datacollectors/CMakeLists.txt
+++ b/dune/vtk/datacollectors/CMakeLists.txt
@@ -8,3 +8,5 @@ install(FILES
unstructureddatacollector.hh
yaspdatacollector.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtkwriter/datacollectors)
+
+add_subdirectory(test)
\ No newline at end of file
diff --git a/dune/vtk/datacollectors/lagrangedatacollector.hh b/dune/vtk/datacollectors/lagrangedatacollector.hh
new file mode 100644
index 0000000..d38f341
--- /dev/null
+++ b/dune/vtk/datacollectors/lagrangedatacollector.hh
@@ -0,0 +1,196 @@
+#pragma once
+
+#include <cassert>
+#include <map>
+#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/utility/lagrangepoints.hh>
+
+#include "unstructureddatacollector.hh"
+
+namespace Dune {
+
+/// Implementation of \ref DataCollector for lagrange cells
+template <class GridView, int ORDER = -1>
+class LagrangeDataCollector
+ : public UnstructuredDataCollectorInterface<GridView, LagrangeDataCollector<GridView,ORDER>, Partitions::All>
+{
+ using Self = LagrangeDataCollector;
+ using Super = UnstructuredDataCollectorInterface<GridView, Self, Partitions::All>;
+
+public:
+ static_assert(ORDER != 0, "Order 0 not supported");
+ using Super::dim;
+ using Super::partition; // NOTE: lagrange data-collector currently implemented for the All partition only
+
+public:
+ LagrangeDataCollector (GridView const& gridView, int order = (ORDER < 0 ? 2 : ORDER))
+ : Super(gridView)
+ , order_(order)
+ {
+ assert(order > 0 && "Order 0 not supported");
+ assert(ORDER < 0 || order == ORDER);
+ }
+
+ /// Construct the point sets
+ void updateImpl ()
+ {
+ auto const& indexSet = gridView_.indexSet();
+
+ pointSets_.clear();
+ for (auto gt : indexSet.types(0))
+ pointSets_.emplace(gt, order_);
+
+ for (auto& pointSet : pointSets_)
+ pointSet.second.build(pointSet.first);
+
+ numPoints_ = indexSet.size(dim);
+ for (auto const& pointSet : pointSets_) {
+ auto gt = pointSet.first;
+ auto refElem = referenceElement<double,dim>(gt);
+ numPoints_ += (pointSet.second.size() - refElem.size(dim)) * indexSet.size(gt);
+ }
+ }
+
+ /// Return number of lagrange nodes
+ std::uint64_t numPointsImpl () const
+ {
+ return numPoints_;
+ }
+
+ /// Return a vector of point coordinates.
+ /**
+ * The vector of point coordinates is composed of vertex coordinates first and second
+ * edge center coordinates.
+ **/
+ template <class T>
+ std::vector<T> pointsImpl () const
+ {
+ std::vector<T> data(this->numPoints() * 3);
+ auto const& indexSet = gridView_.indexSet();
+
+ std::size_t shift = indexSet.size(dim);
+
+ for (auto const& element : elements(gridView_, partition)) {
+ auto geometry = element.geometry();
+ auto refElem = referenceElement<T,dim>(element.type());
+
+ auto const& pointSet = pointSets_.at(element.type());
+ unsigned int vertexDOFs = refElem.size(dim);
+ unsigned int innerDOFs = pointSet.size() - vertexDOFs;
+
+ for (std::size_t i = 0; i < pointSet.size(); ++i) {
+ auto const& p = pointSet[i];
+ if (i < vertexDOFs)
+ assert(p.localKey().codim() == dim);
+
+ auto const& localKey = p.localKey();
+ std::size_t idx = 3 * (localKey.codim() == dim
+ ? indexSet.subIndex(element, localKey.subEntity(), dim)
+ : innerDOFs*indexSet.index(element) + (i - vertexDOFs) + shift);
+
+ auto v = geometry.global(p.point());
+ 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;
+ }
+
+ /// Return number of grid cells
+ std::uint64_t numCellsImpl () const
+ {
+ return gridView_.size(0);
+ }
+
+ /// \brief Return cell types, offsets, and connectivity. \see Cells
+ /**
+ * The cell connectivity is composed of cell vertices first and second cell edges,
+ * where the indices are grouped [vertex-indices..., (#vertices)+edge-indices...]
+ **/
+ Cells cellsImpl () const
+ {
+ Cells cells;
+ cells.connectivity.reserve(this->numPoints());
+ cells.offsets.reserve(this->numCells());
+ cells.types.reserve(this->numCells());
+
+ auto const& indexSet = gridView_.indexSet();
+ std::size_t shift = indexSet.size(dim);
+
+ std::int64_t old_o = 0;
+ for (auto const& element : elements(gridView_, partition)) {
+ auto refElem = referenceElement<double,dim>(element.type());
+ Vtk::CellType cellType(element.type(), Vtk::LAGRANGE);
+
+ auto const& pointSet = pointSets_.at(element.type());
+ unsigned int vertexDOFs = refElem.size(dim);
+ unsigned int innerDOFs = pointSet.size() - vertexDOFs;
+
+ for (std::size_t i = 0; i < pointSet.size(); ++i) {
+ auto const& p = pointSet[i];
+ auto const& localKey = p.localKey();
+ std::size_t idx = (localKey.codim() == dim
+ ? indexSet.subIndex(element, localKey.subEntity(), dim)
+ : innerDOFs*indexSet.index(element) + (i - vertexDOFs) + shift);
+ cells.connectivity.push_back(idx);
+ }
+
+ cells.offsets.push_back(old_o += pointSet.size());
+ cells.types.push_back(cellType.type());
+ }
+ return cells;
+ }
+
+ /// Evaluate the `fct` at element vertices and edge centers in the same order as the point coords.
+ template <class T, class GlobalFunction>
+ std::vector<T> pointDataImpl (GlobalFunction const& fct) const
+ {
+ int nComps = fct.ncomps();
+ std::vector<T> data(this->numPoints() * nComps);
+ auto const& indexSet = gridView_.indexSet();
+
+ std::size_t shift = indexSet.size(dim);
+
+ auto localFct = localFunction(fct);
+ for (auto const& element : elements(gridView_, partition)) {
+ localFct.bind(element);
+ auto refElem = referenceElement<T,dim>(element.type());
+
+ auto const& pointSet = pointSets_.at(element.type());
+ unsigned int vertexDOFs = refElem.size(dim);
+ unsigned int innerDOFs = pointSet.size() - vertexDOFs;
+
+ for (std::size_t i = 0; i < pointSet.size(); ++i) {
+ auto const& p = pointSet[i];
+ auto const& localKey = p.localKey();
+ std::size_t idx = nComps * (localKey.codim() == dim
+ ? indexSet.subIndex(element, localKey.subEntity(), dim)
+ : innerDOFs*indexSet.index(element) + (i - vertexDOFs) + shift);
+
+ for (int comp = 0; comp < nComps; ++comp)
+ data[idx + comp] = T(localFct.evaluate(comp, p.point()));
+ }
+ localFct.unbind();
+ }
+ return data;
+ }
+
+protected:
+ using Super::gridView_;
+
+ unsigned int order_;
+ std::uint64_t numPoints_ = 0;
+
+ using PointSet = VtkLagrangePointSet<typename GridView::ctype, GridView::dimension>;
+ std::map<GeometryType, PointSet> pointSets_;
+};
+
+} // end namespace Dune
diff --git a/dune/vtk/datacollectors/test/CMakeLists.txt b/dune/vtk/datacollectors/test/CMakeLists.txt
new file mode 100644
index 0000000..64403b5
--- /dev/null
+++ b/dune/vtk/datacollectors/test/CMakeLists.txt
@@ -0,0 +1,2 @@
+dune_add_test(SOURCES test-lagrangedatacollector.cc
+ LINK_LIBRARIES dunevtk)
\ No newline at end of file
diff --git a/dune/vtk/datacollectors/test/test-lagrangedatacollector.cc b/dune/vtk/datacollectors/test/test-lagrangedatacollector.cc
new file mode 100644
index 0000000..99d22c5
--- /dev/null
+++ b/dune/vtk/datacollectors/test/test-lagrangedatacollector.cc
@@ -0,0 +1,119 @@
+// -*- 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 <dune/common/hybridutilities.hh>
+#include <dune/common/rangeutilities.hh>
+#include <dune/common/test/testsuite.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#if HAVE_UG
+ #include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/vtk/datacollectors/continuousdatacollector.hh>
+#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
+
+/**
+ * This test compares LagrangeDataCollector<1> against the ContinuousDataCollector, thus
+ * just compares the linear order implementations
+ **/
+
+template <class DataCollector1, class DataCollector2>
+void testDataCollector (std::string prefix, Dune::TestSuite& testSuite, DataCollector1& dataCollector1, DataCollector2& dataCollector2)
+{
+ using std::sqrt;
+ auto tol = sqrt(std::numeric_limits<double>::epsilon());
+
+ dataCollector1.update();
+ dataCollector2.update();
+
+ // check that number of points and cells are equal
+ testSuite.check(dataCollector1.numPoints() == dataCollector2.numPoints(), prefix + "_numPoints");
+ testSuite.check(dataCollector1.numCells() == dataCollector2.numCells(), prefix + "_numCells");
+
+ auto points1 = dataCollector1.template points<double>();
+ auto points2 = dataCollector2.template points<double>();
+
+ // check that point sizes are equal
+ testSuite.check(points1.size() == points2.size(), prefix + "_points.size");
+ testSuite.check(points1.size() == 3*dataCollector1.numPoints(), prefix + "_points.size/3");
+
+ // check that point coordinates are equal
+ using std::abs;
+ for (std::size_t i = 0; i < points1.size(); ++i)
+ testSuite.check(abs(points1[i] - points2[i]) < tol, prefix + "_points[" + std::to_string(i) + "]");
+
+ auto cells1 = dataCollector1.cells();
+ auto cells2 = dataCollector2.cells();
+
+ // check that cell sizes are equal
+ testSuite.check(cells1.types.size() == cells2.types.size(), prefix + "_cells.types.size");
+ testSuite.check(cells1.offsets.size() == cells2.offsets.size(), prefix + "_cells.offsets.size");
+ testSuite.check(cells1.connectivity.size() == cells2.connectivity.size(), prefix + "_cells.connectivity.size");
+
+ // NOTE: cells.types do not need to be equal, e.g. LINEAR != LAGRANGE_LINEAR
+
+ // check that offsets are equal
+ for (std::size_t i = 0; i < cells1.offsets.size(); ++i)
+ testSuite.check(cells1.offsets[i] == cells2.offsets[i], prefix + "_cells.offsets[" + std::to_string(i) + "]");
+
+ // check that connectivities are equal
+ for (std::size_t i = 0; i < cells1.connectivity.size(); ++i)
+ testSuite.check(cells1.connectivity[i] == cells2.connectivity[i], prefix + "_cells.connectivity[" + std::to_string(i) + "]");
+}
+
+template <class GridView>
+void testGridView (std::string prefix, Dune::TestSuite& testSuite, GridView const& gridView)
+{
+ // 1. test linear order lagrange data-collector
+ {
+ Dune::ContinuousDataCollector<GridView> linearDataCollector(gridView);
+
+ // data collector with template order
+ Dune::LagrangeDataCollector<GridView, 1> lagrangeDataCollector1a(gridView);
+ testDataCollector(prefix + "_linear_template", testSuite, lagrangeDataCollector1a, linearDataCollector);
+
+ // data collector with runtime order
+ Dune::LagrangeDataCollector<GridView> lagrangeDataCollector1b(gridView, 1);
+ testDataCollector(prefix + "_linear_runtime", testSuite, lagrangeDataCollector1b, linearDataCollector);
+ }
+}
+
+template <class Grid>
+void testGrid (std::string prefix, Dune::TestSuite& testSuite, Grid& grid)
+{
+ grid.globalRefine(1);
+
+ testGridView(prefix + "_level0", testSuite, grid.levelGridView(0));
+ testGridView(prefix + "_leaf", testSuite, grid.leafGridView());
+}
+
+int main(int argc, char** argv)
+{
+ using namespace Dune;
+ MPIHelper::instance(argc, argv);
+
+ TestSuite testSuite;
+
+ YaspGrid<1> yaspGrid1({1.0}, {1});
+ YaspGrid<2> yaspGrid2({1.0,1.0}, {1,1});
+ YaspGrid<3> yaspGrid3({1.0,1.0,1.0}, {1,1,1});
+
+ testGrid("yaspgrid_1d", testSuite, yaspGrid1);
+ testGrid("yaspgrid_2d", testSuite, yaspGrid2);
+ testGrid("yaspgrid_3d", testSuite, yaspGrid3);
+
+#if HAVE_UG
+ auto ugGrid2 = StructuredGridFactory<UGGrid<2>>::createSimplexGrid({0.0,0.0}, {1.0,1.0}, {1u,1u});
+ auto ugGrid3 = StructuredGridFactory<UGGrid<3>>::createSimplexGrid({0.0,0.0,0.0}, {1.0,1.0,1.0}, {1u,1u,1u});
+
+ testGrid("uggrid_2d", testSuite, *ugGrid2);
+ testGrid("uggrid_3d", testSuite, *ugGrid3);
+#endif
+
+ return testSuite.exit();
+}
diff --git a/dune/vtk/forward.hh b/dune/vtk/forward.hh
index 3b6ac80..f3946ea 100644
--- a/dune/vtk/forward.hh
+++ b/dune/vtk/forward.hh
@@ -65,6 +65,9 @@ namespace Dune
template <class Grid>
struct SerialGridCreator;
+
+ template <class Grid>
+ class LagrangeGridCreator;
// @} gridcreators
diff --git a/dune/vtk/gridcreatorinterface.hh b/dune/vtk/gridcreatorinterface.hh
index c8c8fb4..02def76 100644
--- a/dune/vtk/gridcreatorinterface.hh
+++ b/dune/vtk/gridcreatorinterface.hh
@@ -17,15 +17,17 @@ namespace Dune
/**
* Construct a grid from data read from VTK files.
*
- * \tparam GridView Model of Dune::GridView
- * \tparam Derived Implementation of a concrete GridCreator.
+ * \tparam GridType Model of Dune::Grid
+ * \tparam GlobalCoordType Type of the global coordinates.
+ * \tparam DerivedType Implementation of a concrete GridCreator.
**/
- template <class G, class Derived>
+ template <class GridType, class DerivedType>
class GridCreatorInterface
{
public:
- using Grid = G;
+ using Grid = GridType;
using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
+ using Derived = DerivedType;
public:
/// Constructor. Stores a reference to the passed GridFactory
@@ -60,6 +62,12 @@ namespace Dune
return *factory_;
}
+ /// Return the associated (const) GridFactory
+ GridFactory<Grid> const& factory () const
+ {
+ return *factory_;
+ }
+
/// Return the mpi collective communicator
auto comm () const
{
diff --git a/dune/vtk/gridcreators/CMakeLists.txt b/dune/vtk/gridcreators/CMakeLists.txt
index 299b4d3..b605ec6 100644
--- a/dune/vtk/gridcreators/CMakeLists.txt
+++ b/dune/vtk/gridcreators/CMakeLists.txt
@@ -4,6 +4,7 @@ install(FILES
continuousgridcreator.hh
derivedgridcreator.hh
discontinuousgridcreator.hh
+ lagrangegridcreator.hh
parallelgridcreator.hh
serialgridcreator.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtkwriter/gridcreators)
diff --git a/dune/vtk/gridcreators/continuousgridcreator.hh b/dune/vtk/gridcreators/continuousgridcreator.hh
index 03270a2..3db1d0c 100644
--- a/dune/vtk/gridcreators/continuousgridcreator.hh
+++ b/dune/vtk/gridcreators/continuousgridcreator.hh
@@ -20,7 +20,8 @@ namespace Dune
struct ContinuousGridCreator
: public GridCreatorInterface<Grid, ContinuousGridCreator<Grid>>
{
- using Super = GridCreatorInterface<Grid, ContinuousGridCreator<Grid>>;
+ using Self = ContinuousGridCreator;
+ using Super = GridCreatorInterface<Grid, ContinuousGridCreator>;
using GlobalCoordinate = typename Super::GlobalCoordinate;
ContinuousGridCreator (GridFactory<Grid>& factory)
diff --git a/dune/vtk/gridcreators/discontinuousgridcreator.hh b/dune/vtk/gridcreators/discontinuousgridcreator.hh
index 5160b00..c285ba5 100644
--- a/dune/vtk/gridcreators/discontinuousgridcreator.hh
+++ b/dune/vtk/gridcreators/discontinuousgridcreator.hh
@@ -19,7 +19,8 @@ namespace Dune
struct DiscontinuousGridCreator
: public GridCreatorInterface<Grid, DiscontinuousGridCreator<Grid>>
{
- using Super = GridCreatorInterface<Grid, DiscontinuousGridCreator<Grid>>;
+ using Self = DiscontinuousGridCreator;
+ using Super = GridCreatorInterface<Grid, DiscontinuousGridCreator>;
using GlobalCoordinate = typename Super::GlobalCoordinate;
struct CoordLess
diff --git a/dune/vtk/gridcreators/lagrangegridcreator.hh b/dune/vtk/gridcreators/lagrangegridcreator.hh
new file mode 100644
index 0000000..855b7c4
--- /dev/null
+++ b/dune/vtk/gridcreators/lagrangegridcreator.hh
@@ -0,0 +1,344 @@
+#pragma once
+
+#include <cassert>
+#include <cstdint>
+#include <limits>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/std/optional.hh>
+#include <dune/geometry/utility/typefromvertexcount.hh>
+#include <dune/localfunctions/lagrange.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/vtk/forward.hh>
+#include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/gridcreatorinterface.hh>
+#include <dune/vtk/utility/lagrangepoints.hh>
+
+namespace Dune
+{
+ // \brief Create a grid from data that represents higher (lagrange) cells.
+ /**
+ * The grid is created from the first nodes of a cell parametrization, representing
+ * the corner vertices. Thus a piecewise "flat" grid is constructed. The parametrization
+ * is 1. passed as a local element parametrization to the `insertElement()` function of a
+ * gridFactory to allow the grid itself to handle the parametrization and 2. is stored
+ * internally that can be accessed by using this GridCreator object as a grid function,
+ * or by extracting locally the parametrization on each existing grid element after
+ * creation of the grid.
+ *
+ * So, the LagrangeGridCreator models both, a `GridCreator` and a `GridFunction`.
+ **/
+ template <class Grid>
+ struct LagrangeGridCreator
+ : public GridCreatorInterface<Grid, LagrangeGridCreator<Grid>>
+ {
+ using Self = LagrangeGridCreator;
+ using Super = GridCreatorInterface<Grid, Self>;
+ using GlobalCoordinate = typename Super::GlobalCoordinate;
+
+ using Nodes = std::vector<GlobalCoordinate>;
+
+ struct ElementParametrization
+ {
+ GeometryType type; //< Geometry type of the element
+ std::vector<std::int64_t> nodes; //< Indices of the w.r.t. `nodes_` vector
+ std::vector<unsigned int> corners; //< Insertion-indices of the element corner nodes
+ };
+
+ using Parametrization = std::vector<ElementParametrization>;
+ using Element = typename Grid::template Codim<0>::Entity;
+ using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+
+ class LocalParametrization;
+ class LocalFunction;
+
+ public:
+ using Super::factory;
+
+ LagrangeGridCreator (GridFactory<Grid>& factory)
+ : Super(factory)
+ {}
+
+ /// Implementation of the interface function `insertVertices()`
+ void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
+ std::vector<std::uint64_t> const& /*point_ids*/)
+ {
+ // store point coordinates in member variable
+ nodes_ = points;
+ }
+
+ template <class F>
+ using HasParametrizedElements = decltype(std::declval<F>().insertElement(std::declval<GeometryType>(),
+ std::declval<std::vector<unsigned int> const&>(), std::declval<std::function<GlobalCoordinate(LocalCoordinate)>>()));
+
+ /// Implementation of the interface function `insertElements()`
+ void insertElementsImpl (std::vector<std::uint8_t> const& types,
+ std::vector<std::int64_t> const& offsets,
+ std::vector<std::int64_t> const& connectivity)
+ {
+ assert(nodes_.size() > 0);
+
+ // mapping of node index to element-vertex index
+ std::vector<std::int64_t> elementVertices(nodes_.size(), -1);
+ parametrization_.resize(types.size());
+
+ std::int64_t vertexIndex = 0;
+ for (std::size_t i = 0; i < types.size(); ++i) {
+ auto type = Vtk::to_geometry(types[i]);
+ Vtk::CellType cellType{type};
+ auto refElem = referenceElement<double,Grid::dimension>(type);
+
+ std::int64_t shift = (i == 0 ? 0 : offsets[i-1]);
+ int nNodes = offsets[i] - shift;
+ int nVertices = refElem.size(Grid::dimension);
+
+ // insert vertices into grid and construct element vertices
+ std::vector<unsigned int> element(nVertices);
+ for (int j = 0; j < nVertices; ++j) {
+ auto index = connectivity.at(shift + j);
+ auto& vertex = elementVertices.at(index);
+ if (vertex < 0) {
+ factory().insertVertex(nodes_.at(index));
+ vertex = vertexIndex++;
+ }
+ element[j] = vertex;
+ }
+
+ // permute element indices
+ if (!cellType.noPermutation()) {
+ // apply index permutation
+ std::vector<unsigned int> cell(element.size());
+ for (int j = 0; j < element.size(); ++j)
+ cell[j] = element[cellType.permutation(j)];
+ std::swap(element, cell);
+ }
+
+ // fill vector of element parametrizations
+ parametrization_[i].type = type;
+ parametrization_[i].nodes.resize(nNodes);
+ for (int j = 0; j < nNodes; ++j)
+ parametrization_[i].nodes[j] = connectivity.at(shift + j);
+ parametrization_[i].corners = element;
+
+ // try to create element with parametrization
+ if constexpr (Std::is_detected_v<HasParametrizedElements, GridFactory<Grid>>) {
+ try {
+ factory().insertElement(type, element, localParametrization(i));
+ } catch (Dune::GridError const& /* notImplemented */) {
+ factory().insertElement(type, element);
+ }
+ } else {
+ factory().insertElement(type, element);
+ }
+ }
+ }
+
+ /// \brief Construct an element parametrization
+ /**
+ * The returned LocalParametrization is a mapping `GlobalCoordinate(LocalCoordinate)`
+ * where `LocalCoordinate is w.r.t. the local coordinate system in an element with
+ * given `insertionIndex` (defined by the inserted corner vertices) and
+ * `GlobalCoordinate` a world coordinate in the parametrized grid.
+ **/
+ LocalParametrization localParametrization (unsigned int insertionIndex) const
+ {
+ assert(!nodes_.empty() && !parametrization_.empty());
+ auto const& localParam = parametrization_.at(insertionIndex);
+ return LocalParametrization{nodes_, localParam, order(localParam)};
+ }
+
+ /// \brief Construct an element parametrization
+ /**
+ * The returned LocalParametrization is a mapping `GlobalCoordinate(LocalCoordinate)`
+ * where `LocalCoordinate is w.r.t. the local coordinate system in the passed element
+ * and `GlobalCoordinate` a world coordinate in the parametrized grid.
+ *
+ * Note, when an element is passed, it might have a different local coordinate system
+ * than the coordinate system used to defined the element parametrization. Thus coordinate
+ * transform is internally chained to the evaluation of the local parametrization. This
+ * local geometry transform is obtained by figuring out the permutation of corners in the
+ * element corresponding to the inserted corner vertices.
+ **/
+ LocalParametrization localParametrization (Element const& element) const
+ {
+ assert(!nodes_.empty() && !parametrization_.empty());
+
+ unsigned int insertionIndex = factory().insertionIndex(element);
+ auto const& localParam = parametrization_.at(insertionIndex);
+ assert(element.type() == localParam.type);
+
+ // collect indices of vertices
+ std::vector<unsigned int> indices(element.subEntities(Grid::dimension));
+ for (unsigned int i = 0; i < element.subEntities(Grid::dimension); ++i)
+ indices[i] = factory().insertionIndex(element.template subEntity<Grid::dimension>(i));
+
+ // calculate permutation vector
+ std::vector<unsigned int> permutation(indices.size());
+ for (std::size_t i = 0; i < indices.size(); ++i) {
+ auto it = std::find(localParam.corners.begin(), localParam.corners.end(), indices[i]);
+ assert(it != localParam.corners.end());
+ permutation[i] = std::distance(localParam.corners.begin(), it);
+ }
+
+ return LocalParametrization{nodes_, localParam, order(localParam), permutation};
+ }
+
+ /// \brief Local function representing the parametrization of the grid.
+ /**
+ * The returned LocalFunction object models a Functions::Concept::LocalFunction
+ * and can thus be bound to an element of the created grid and evaluated in the
+ * local coordinates of the bound element.
+ *
+ * It is implemented in terms of the \ref LocalParametrization function returned by
+ * the method \ref localParametrization(element). See comments there for furhter
+ * details.
+ *
+ * Note, this methods requires the GridCreator to be bassed by lvalue-reference. This
+ * is necessary, since we want to guarantee that all internal storage is preserved while
+ * evaluating the local function.
+ **/
+ friend LocalFunction localFunction (LagrangeGridCreator& gridCreator)
+ {
+ return LocalFunction{gridCreator};
+ }
+
+ /// Determine lagrange order from number of points
+ template <class LocalParam>
+ unsigned int order (LocalParam const& localParam) const
+ {
+ GeometryType type = localParam.type;
+ std::size_t nNodes = localParam.nodes.size();
+ for (unsigned int o = 1; o <= nNodes; ++o)
+ if (numLagrangePoints(type.id(), type.dim(), o) == nNodes)
+ return o;
+
+ return 1;
+ }
+
+ /// Determine lagrange order from number of points from the first element parametrization
+ unsigned int order () const
+ {
+ assert(!parametrization_.empty());
+ return order(parametrization_.front());
+ }
+
+ private:
+ /// All point coordinates inclusing the higher-order lagrange points
+ Nodes nodes_;
+
+ /// Parametrization for all elements
+ Parametrization parametrization_;
+ };
+
+
+ template <class Grid>
+ class LagrangeGridCreator<Grid>::LocalParametrization
+ {
+ using ctype = typename Grid::ctype;
+
+ using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
+ using LocalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::LocalCoordinate;
+ using LocalGeometry = MultiLinearGeometry<ctype,Grid::dimension,Grid::dimension>;
+
+ using LocalFE = LagrangeLocalFiniteElement<VtkLagrangePointSet, Grid::dimension, ctype, ctype>;
+ using LocalBasis = typename LocalFE::Traits::LocalBasisType;
+ using LocalBasisTraits = typename LocalBasis::Traits;
+
+ public:
+ /// Construct a local element parametrization
+ template <class Nodes, class LocalParam>
+ LocalParametrization (Nodes const& nodes, LocalParam const& param, unsigned int order)
+ : localFE_(param.type, order)
+ , localNodes_(param.nodes.size())
+ {
+ for (std::size_t i = 0; i < localNodes_.size(); ++i)
+ localNodes_[i] = nodes[param.nodes[i]];
+ }
+
+ /// Construct a local element parametrization for elements with permuted corners
+ template <class Nodes, class LocalParam, class Permutation>
+ LocalParametrization (Nodes const& nodes, LocalParam const& param, unsigned int order, Permutation const& permutation)
+ : LocalParametrization(nodes, param, order)
+ {
+ auto refElem = referenceElement<ctype,Grid::dimension>(param.type);
+ std::vector<LocalCoordinate> corners(permutation.size());
+ for (std::size_t i = 0; i < permutation.size(); ++i)
+ corners[i] = refElem.position(permutation[i], Grid::dimension);
+
+ localGeometry_.emplace(param.type, corners);
+ }
+
+ /// Evaluate the local parametrization in local coordinates
+ template <class LocalCoordinate>
+ GlobalCoordinate operator() (LocalCoordinate const& local) const
+ {
+ // map coordinates if element corners are permuted
+ LocalCoordinate x = localGeometry_ ? localGeometry_->global(local) : local;
+
+ LocalBasis const& localBasis = localFE_.localBasis();
+ localBasis.evaluateFunction(x, shapeValues_);
+ assert(shapeValues_.size() == localNodes_.size());
+
+ GlobalCoordinate out(0);
+ for (std::size_t i = 0; i < shapeValues_.size(); ++i)
+ out.axpy(shapeValues_[i], localNodes_[i]);
+
+ return out;
+ }
+
+ private:
+ LocalFE localFE_;
+ std::vector<GlobalCoordinate> localNodes_;
+ Std::optional<LocalGeometry> localGeometry_;
+
+ mutable std::vector<typename LocalBasisTraits::RangeType> shapeValues_;
+ };
+
+
+ template <class Grid>
+ class LagrangeGridCreator<Grid>::LocalFunction
+ {
+ using ctype = typename Grid::ctype;
+ using LocalContext = typename Grid::template Codim<0>::Entity;
+ using GlobalCoordinate = typename LocalContext::Geometry::GlobalCoordinate;
+ using LocalCoordinate = typename LocalContext::Geometry::LocalCoordinate;
+ using LocalParametrization = typename LagrangeGridCreator::LocalParametrization;
+
+ public:
+ explicit LocalFunction (LagrangeGridCreator& gridCreator)
+ : gridCreator_(&gridCreator)
+ {}
+
+ /// Collect a local parametrization on the element
+ void bind (LocalContext const& element)
+ {
+ localContext_ = element;
+ localParametrization_.emplace(gridCreator_->localParametrization(element));
+ }
+
+ void unbind () { /* do nothing */ }
+
+ /// Evaluate the local parametrization in local coordinates
+ GlobalCoordinate operator() (LocalCoordinate const& local) const
+ {
+ assert(!!localParametrization_);
+ return (*localParametrization_)(local);
+ }
+
+ /// Return the bound element
+ LocalContext const& localContext () const
+ {
+ return localContext_;
+ }
+
+ private:
+ LagrangeGridCreator const* gridCreator_;
+
+ LocalContext localContext_;
+ Std::optional<LocalParametrization> localParametrization_;
+ };
+
+} // end namespace Dune
diff --git a/dune/vtk/utility/CMakeLists.txt b/dune/vtk/utility/CMakeLists.txt
index e376849..cdc7dce 100644
--- a/dune/vtk/utility/CMakeLists.txt
+++ b/dune/vtk/utility/CMakeLists.txt
@@ -5,6 +5,10 @@ dune_add_library("filesystem" OBJECT
install(FILES
enum.hh
filesystem.hh
+ lagrangelfecache.hh
+ lagrangepoints.hh
string.hh
uid.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtkwriter/utility)
+
+ add_subdirectory(test)
\ No newline at end of file
diff --git a/dune/vtk/utility/lagrangelfecache.hh b/dune/vtk/utility/lagrangelfecache.hh
new file mode 100644
index 0000000..ac2db68
--- /dev/null
+++ b/dune/vtk/utility/lagrangelfecache.hh
@@ -0,0 +1,35 @@
+#pragma once
+
+#include <map>
+
+#include <dune/geometry/type.hh>
+#include <dune/localfunctions/lagrange.hh>
+#include <dune/vtk/utility/lagrangepoints.hh>
+
+namespace Dune
+{
+ // LocalFiniteElement for the lagrange pointset compatible with VTK
+ template <class Domain, class Range, int dim, int ORDER=-1>
+ class LagrangeLFECache
+ {
+ public:
+ using FiniteElementType = LagrangeLocalFiniteElement<VtkLagrangePointSet, dim, Domain, Range>;
+
+ LagrangeLFECache (unsigned int order = ORDER)
+ : order_(order)
+ {}
+
+ FiniteElementType const& get (GeometryType type)
+ {
+ auto it = data_.find(type);
+ if (it == data_.end())
+ it = data_.emplace(type,FiniteElementType(type,order_)).first;
+ return it->second;
+ }
+
+ private:
+ unsigned int order_;
+ std::map<GeometryType, FiniteElementType> data_;
+ };
+
+} // end namespace Dune
diff --git a/dune/vtk/utility/lagrangepoints.hh b/dune/vtk/utility/lagrangepoints.hh
new file mode 100644
index 0000000..b598c89
--- /dev/null
+++ b/dune/vtk/utility/lagrangepoints.hh
@@ -0,0 +1,605 @@
+#pragma once
+
+#include <cassert>
+#include <array>
+
+#include <dune/common/exceptions.hh>
+#include <dune/geometry/type.hh>
+#include <dune/localfunctions/lagrange/equidistantpoints.hh>
+
+namespace Dune {
+
+namespace Impl {
+ // forward declaration
+ template <class K, unsigned int dim>
+ class VtkLagrangePointSetBuilder;
+}
+
+
+/// \brief A set of lagrange points compatible with the numbering of VTK
+/**
+ * \tparam K Field-type for the coordinates
+ * \tparam dim Dimension of the coordinates
+ **/
+template <class K, unsigned int dim>
+class VtkLagrangePointSet
+ : public EmptyPointSet<K, dim>
+{
+ using Super = EmptyPointSet<K, dim>;
+
+public:
+ static const unsigned int dimension = dim;
+
+ VtkLagrangePointSet (std::size_t order)
+ : Super(order)
+ {
+ assert(order > 0);
+ }
+
+ /// Fill the lagrange points for the given geometry type
+ void build (GeometryType gt)
+ {
+ assert(gt.dim() == dimension);
+ builder_(gt, order(), points_);
+ }
+
+ /// Fill the lagrange points for the given topology type `Topology`
+ template <class Topology>
+ bool build ()
+ {
+ build(GeometryType(Topology{}));
+ return true;
+ }
+
+ /// Returns whether the point set support the given topology type `Topology` and can
+ /// generate point for the given order.
+ template <class Topology>
+ static bool supports (std::size_t order)
+ {
+ return true;
+ }
+
+ using Super::order;
+
+private:
+ using Super::points_;
+ Impl::VtkLagrangePointSetBuilder<K,dim> builder_;
+};
+
+
+namespace Impl {
+
+// Build for lagrange point sets in different dimensions
+// Specialized for dim=1,2,3
+template <class K, unsigned int dim>
+class VtkLagrangePointSetBuilder
+{
+public:
+ template <class Points>
+ void operator()(GeometryType, unsigned int, Points& points) const
+ {
+ DUNE_THROW(Dune::NotImplemented,
+ "Lagrange points not yet implemented for this GeometryType.");
+ }
+};
+
+/**
+ * The implementation of the point set builder is directly derived from VTK.
+ * Modification are a change in data-types and naming scheme. Additionally
+ * a LocalKey is created to reflect the concept of a Dune PointSet.
+ *
+ * Included is the license of the BSD 3-clause License included in the VTK
+ * source code from 2020/04/13 in commit b90dad558ce28f6d321420e4a6b17e23f5227a1c
+ * of git repository https://gitlab.kitware.com/vtk/vtk.
+ *
+ Program: Visualization Toolkit
+ Module: Copyright.txt
+
+ Copyright (c) 1993-2015 Ken Martin, Will Schroeder, Bill Lorensen
+ All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice,
+ this list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright notice,
+ this list of conditions and the following disclaimer in the documentation
+ and/or other materials provided with the distribution.
+
+ * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
+ of any contributors may be used to endorse or promote products derived
+ from this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
+ AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
+ ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+ OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ **/
+
+// Lagrange points on point geometries
+template <class K>
+class VtkLagrangePointSetBuilder<K,0>
+{
+ static constexpr int dim = 0;
+ using LP = LagrangePoint<K,dim>;
+ using Vec = typename LP::Vector;
+ using Key = LocalKey;
+
+public:
+ template <class Points>
+ void operator()(GeometryType gt, int /*order*/, Points& points) const
+ {
+ assert(gt.isVertex());
+ points.push_back(LP{Vec{},Key{0,0,0}});
+ }
+};
+
+
+// Lagrange points on line geometries
+template <class K>
+class VtkLagrangePointSetBuilder<K,1>
+{
+ static constexpr int dim = 1;
+ using LP = LagrangePoint<K,dim>;
+ using Vec = typename LP::Vector;
+ using Key = LocalKey;
+
+public:
+ template <class Points>
+ void operator()(GeometryType gt, int order, Points& points) const
+ {
+ assert(gt.isLine());
+
+ // Vertex nodes
+ points.push_back(LP{Vec{0.0}, Key{0,dim,0}});
+ points.push_back(LP{Vec{1.0}, Key{1,dim,0}});
+
+ if (order > 1) {
+ // Inner nodes
+ Vec p{0.0};
+ for (unsigned int i = 0; i < order-1; i++)
+ {
+ p[0] += 1.0 / order;
+ points.push_back(LP{p,Key{0,dim-1,i}});
+ }
+ }
+ }
+};
+
+
+// Lagrange points on 2d geometries
+template <class K>
+class VtkLagrangePointSetBuilder<K,2>
+{
+ static constexpr int dim = 2;
+ using LP = LagrangePoint<K,dim>;
+ using Vec = typename LP::Vector;
+ using Key = LocalKey;
+
+ friend class VtkLagrangePointSetBuilder<K,3>;
+
+public:
+ template <class Points>
+ void operator()(GeometryType gt, int order, Points& points) const
+ {
+ std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
+
+ if (gt.isTriangle())
+ buildTriangle(nPoints, order, points);
+ else if (gt.isQuadrilateral())
+ buildQuad(nPoints, order, points);
+ else {
+ DUNE_THROW(Dune::NotImplemented,
+ "Lagrange points not yet implemented for this GeometryType.");
+ }
+
+ assert(points.size() == nPoints);
+ }
+
+private:
+ // Construct the point set in a triangle element.
+ // Loop from the outside to the inside
+ template <class Points>
+ void buildTriangle (std::size_t nPoints, int order, Points& points) const
+ {
+ points.reserve(nPoints);
+
+ const int nVertexDOFs = 3;
+ const int nEdgeDOFs = 3 * (order-1);
+
+ static const unsigned int vertexPerm[3] = {0,1,2};
+ static const unsigned int edgePerm[3] = {0,2,1};
+
+ auto calcKey = [&](int p) -> Key
+ {
+ if (p < nVertexDOFs) {
+ return Key{vertexPerm[p], dim, 0};
+ }
+ else if (p < nVertexDOFs+nEdgeDOFs) {
+ unsigned int entityIndex = (p - nVertexDOFs) / (order-1);
+ unsigned int index = (p - nVertexDOFs) % (order-1);
+ return Key{edgePerm[entityIndex], dim-1, index};
+ }
+ else {
+ unsigned int index = p - (nVertexDOFs + nEdgeDOFs);
+ return Key{0, dim-2, index};
+ }
+ };
+
+ std::array<int,3> bindex;
+
+ double order_d = double(order);
+ for (std::size_t p = 0; p < nPoints; ++p) {
+ barycentricIndex(p, bindex, order);
+ Vec point{bindex[0] / order_d, bindex[1] / order_d};
+ points.push_back(LP{point, calcKey(p)});
+ }
+ }
+
+ // "Barycentric index" is a triplet of integers, each running from 0 to
+ // <Order>. It is the index of a point on the triangle in barycentric
+ // coordinates.
+ static void barycentricIndex (int index, std::array<int,3>& bindex, int order)
+ {
+ int max = order;
+ int min = 0;
+
+ // scope into the correct triangle
+ while (index != 0 && index >= 3 * order)
+ {
+ index -= 3 * order;
+ max -= 2;
+ min++;
+ order -= 3;
+ }
+
+ // vertex DOFs
+ if (index < 3)
+ {
+ bindex[index] = bindex[(index + 1) % 3] = min;
+ bindex[(index + 2) % 3] = max;
+ }
+ // edge DOFs
+ else
+ {
+ index -= 3;
+ int dim = index / (order - 1);
+ int offset = (index - dim * (order - 1));
+ bindex[(dim + 1) % 3] = min;
+ bindex[(dim + 2) % 3] = (max - 1) - offset;
+ bindex[dim] = (min + 1) + offset;
+ }
+ }
+
+
+ // Construct the point set in the quad element
+ // 1. build equispaced points with index tuple (i,j)
+ // 2. map index tuple to DOF index and LocalKey
+ template <class Points>
+ void buildQuad(std::size_t nPoints, int order, Points& points) const
+ {
+ points.resize(nPoints);
+
+ std::array<int,2> orders{order, order};
+ std::array<Vec,4> nodes{Vec{0., 0.}, Vec{1., 0.}, Vec{1., 1.}, Vec{0., 1.}};
+
+ for (int n = 0; n <= orders[1]; ++n) {
+ for (int m = 0; m <= orders[0]; ++m) {
+ // int idx = pointIndexFromIJ(m,n,orders);
+
+ const double r = double(m) / orders[0];
+ const double s = double(n) / orders[1];
+ Vec p = (1.0 - r) * (nodes[3] * s + nodes[0] * (1.0 - s)) +
+ r * (nodes[2] * s + nodes[1] * (1.0 - s));
+
+ auto [idx,key] = calcQuadKey(m,n,orders);
+ points[idx] = LP{p, key};
+ // points[idx] = LP{p, calcQuadKey(n,m,orders)};
+ }
+ }
+ }
+
+ // Obtain the VTK DOF index of the node (i,j) in the quad element
+ // and construct a LocalKey
+ static std::pair<int,Key> calcQuadKey (int i, int j, std::array<int,2> order)
+ {
+ const bool ibdy = (i == 0 || i == order[0]);
+ const bool jbdy = (j == 0 || j == order[1]);
+ const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0); // How many boundaries do we lie on at once?
+
+ int dof = 0;
+ unsigned int entityIndex = 0;
+ unsigned int index = 0;
+
+ if (nbdy == 2) // Vertex DOF
+ {
+ dof = (i ? (j ? 2 : 1) : (j ? 3 : 0));
+ entityIndex = (j ? (i ? 3 : 2) : (i ? 1 : 0));
+ return std::make_pair(dof,Key{entityIndex, dim, 0});
+ }
+
+ int offset = 4;
+ if (nbdy == 1) // Edge DOF
+ {
+ if (!ibdy) {
+ dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + offset;
+ entityIndex = j ? 3 : 2;
+ index = i-1;
+ }
+ else if (!jbdy) {
+ dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + offset;
+ entityIndex = i ? 1 : 0;
+ index = j-1;
+ }
+ return std::make_pair(dof, Key{entityIndex, dim-1, index});
+ }
+
+ offset += 2 * (order[0]-1 + order[1]-1);
+
+ // nbdy == 0: Face DOF
+ dof = offset + (i - 1) + (order[0]-1) * ((j - 1));
+ Key innerKey = VtkLagrangePointSetBuilder<K,2>::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second;
+ return std::make_pair(dof, Key{0, dim-2, innerKey.index()});
+ }
+};
+
+
+// Lagrange points on 3d geometries
+template <class K>
+class VtkLagrangePointSetBuilder<K,3>
+{
+ static constexpr int dim = 3;
+ using LP = LagrangePoint<K,dim>;
+ using Vec = typename LP::Vector;
+ using Key = LocalKey;
+
+public:
+ template <class Points>
+ void operator() (GeometryType gt, unsigned int order, Points& points) const
+ {
+ std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
+
+ if (gt.isTetrahedron())
+ buildTetra(nPoints, order, points);
+ else if (gt.isHexahedron())
+ buildHex(nPoints, order, points);
+ else {
+ DUNE_THROW(Dune::NotImplemented,
+ "Lagrange points not yet implemented for this GeometryType.");
+ }
+
+ assert(points.size() == nPoints);
+ }
+
+private:
+ // Construct the point set in the tetrahedron element
+ // 1. construct barycentric (index) coordinates
+ // 2. obtains the DOF index, LocalKey and actual coordinate from barycentric index
+ template <class Points>
+ void buildTetra (std::size_t nPoints, int order, Points& points) const
+ {
+ points.reserve(nPoints);
+
+ const int nVertexDOFs = 4;
+ const int nEdgeDOFs = 6 * (order-1);
+ const int nFaceDOFs = 4 * (order-1)*(order-2)/2;
+
+ static const unsigned int vertexPerm[4] = {0,1,2,3};
+ static const unsigned int edgePerm[6] = {0,2,1,3,4,5};
+ static const unsigned int facePerm[4] = {1,2,0,3};
+
+ auto calcKey = [&](int p) -> Key
+ {
+ if (p < nVertexDOFs) {
+ return Key{vertexPerm[p], dim, 0};
+ }
+ else if (p < nVertexDOFs+nEdgeDOFs) {
+ unsigned int entityIndex = (p - nVertexDOFs) / (order-1);
+ unsigned int index = (p - nVertexDOFs) % (order-1);
+ return Key{edgePerm[entityIndex], dim-1, index};
+ }
+ else if (p < nVertexDOFs+nEdgeDOFs+nFaceDOFs) {
+ unsigned int index = (p - (nVertexDOFs + nEdgeDOFs)) % ((order-1)*(order-2)/2);
+ unsigned int entityIndex = (p - (nVertexDOFs + nEdgeDOFs)) / ((order-1)*(order-2)/2);
+ return Key{facePerm[entityIndex], dim-2, index};
+ }
+ else {
+ unsigned int index = p - (nVertexDOFs + nEdgeDOFs + nFaceDOFs);
+ return Key{0, dim-3, index};
+ }
+ };
+
+ std::array<int,4> bindex;
+
+ double order_d = double(order);
+ for (std::size_t p = 0; p < nPoints; ++p) {
+ barycentricIndex(p, bindex, order);
+ Vec point{bindex[0] / order_d, bindex[1] / order_d, bindex[2] / order_d};
+ points.push_back(LP{point, calcKey(p)});
+ }
+ }
+
+ // "Barycentric index" is a set of 4 integers, each running from 0 to
+ // <Order>. It is the index of a point in the tetrahedron in barycentric
+ // coordinates.
+ static void barycentricIndex (std::size_t p, std::array<int,4>& bindex, int order)
+ {
+ const int nVertexDOFs = 4;
+ const int nEdgeDOFs = 6 * (order-1);
+
+ static const int edgeVertices[6][2] = {{0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3}};
+ static const int linearVertices[4][4] = {{0,0,0,1}, {1,0,0,0}, {0,1,0,0}, {0,0,1,0}};
+ static const int vertexMaxCoords[4] = {3,0,1,2};
+ static const int faceBCoords[4][3] = {{0,2,3}, {2,0,1}, {2,1,3}, {1,0,3}};
+ static const int faceMinCoord[4] = {1,3,0,2};
+
+ int max = order;
+ int min = 0;
+
+ // scope into the correct tetra
+ while (p >= 2 * (order * order + 1) && p != 0 && order > 3)
+ {
+ p -= 2 * (order * order + 1);
+ max -= 3;
+ min++;
+ order -= 4;
+ }
+
+ // vertex DOFs
+ if (p < nVertexDOFs)
+ {
+ for (int coord = 0; coord < 4; ++coord)
+ bindex[coord] = (coord == vertexMaxCoords[p] ? max : min);
+ }
+ // edge DOFs
+ else if (p < nVertexDOFs+nEdgeDOFs)
+ {
+ int edgeId = (p - nVertexDOFs) / (order-1);
+ int vertexId = (p - nVertexDOFs) % (order-1);
+ for (int coord = 0; coord < 4; ++coord)
+ {
+ bindex[coord] = min +
+ (linearVertices[edgeVertices[edgeId][0]][coord] * (max - min - 1 - vertexId) +
+ linearVertices[edgeVertices[edgeId][1]][coord] * (1 + vertexId));
+ }
+ }
+ // face DOFs
+ else
+ {
+ int faceId = (p - (nVertexDOFs+nEdgeDOFs)) / ((order-2)*(order-1)/2);
+ int vertexId = (p - (nVertexDOFs+nEdgeDOFs)) % ((order-2)*(order-1)/2);
+
+ std::array<int,3> projectedBIndex;
+ if (order == 3)
+ projectedBIndex[0] = projectedBIndex[1] = projectedBIndex[2] = 0;
+ else
+ VtkLagrangePointSetBuilder<K,2>::barycentricIndex(vertexId, projectedBIndex, order-3);
+
+ for (int i = 0; i < 3; i++)
+ bindex[faceBCoords[faceId][i]] = (min + 1 + projectedBIndex[i]);
+
+ bindex[faceMinCoord[faceId]] = min;
+ }
+ }
+
+private:
+ // Construct the point set in the heyhedral element
+ // 1. build equispaced points with index tuple (i,j,k)
+ // 2. map index tuple to DOF index and LocalKey
+ template <class Points>
+ void buildHex (std::size_t nPoints, int order, Points& points) const
+ {
+ points.resize(nPoints);
+
+ std::array<int,3> orders{order, order, order};
+ std::array<Vec,8> nodes{Vec{0., 0., 0.}, Vec{1., 0., 0.}, Vec{1., 1., 0.}, Vec{0., 1., 0.},
+ Vec{0., 0., 1.}, Vec{1., 0., 1.}, Vec{1., 1., 1.}, Vec{0., 1., 1.}};
+
+ for (int p = 0; p <= orders[2]; ++p) {
+ for (int n = 0; n <= orders[1]; ++n) {
+ for (int m = 0; m <= orders[0]; ++m) {
+ const double r = double(m) / orders[0];
+ const double s = double(n) / orders[1];
+ const double t = double(p) / orders[2];
+ Vec point = (1.0-r) * ((nodes[3] * (1.0-t) + nodes[7] * t) * s + (nodes[0] * (1.0-t) + nodes[4] * t) * (1.0-s)) +
+ r * ((nodes[2] * (1.0-t) + nodes[6] * t) * s + (nodes[1] * (1.0-t) + nodes[5] * t) * (1.0-s));
+
+ auto [idx,key] = calcHexKey(m,n,p,orders);
+ points[idx] = LP{point, key};
+ }
+ }
+ }
+ }
+
+ // Obtain the VTK DOF index of the node (i,j,k) in the hexahedral element
+ static std::pair<int,Key> calcHexKey (int i, int j, int k, std::array<int,3> order)
+ {
+ const bool ibdy = (i == 0 || i == order[0]);
+ const bool jbdy = (j == 0 || j == order[1]);
+ const bool kbdy = (k == 0 || k == order[2]);
+ const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0); // How many boundaries do we lie on at once?
+
+ int dof = 0;
+ unsigned int entityIndex = 0;
+ unsigned int index = 0;
+
+ if (nbdy == 3) // Vertex DOF
+ {
+ dof = (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0);
+ entityIndex = (i ? 1 : 0) + (j ? 2 : 0) + (k ? 4 : 0);
+ return std::make_pair(dof, Key{entityIndex, dim, 0});
+ }
+
+ int offset = 8;
+ if (nbdy == 2) // Edge DOF
+ {
+ entityIndex = (k ? 8 : 4);
+ if (!ibdy)
+ { // On i axis
+ dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
+ index = i;
+ entityIndex += (i ? 1 : 0);
+ }
+ else if (!jbdy)
+ { // On j axis
+ dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
+ index = j;
+ entityIndex += (j ? 3 : 2);
+ }
+ else
+ { // !kbdy, On k axis
+ offset += 4 * (order[0]-1) + 4 * (order[1]-1);
+ dof = (k - 1) + (order[2]-1) * (i ? (j ? 3 : 1) : (j ? 2 : 0)) + offset;
+ index = k;
+ entityIndex = (i ? 1 : 0) + (j ? 2 : 0);
+ }
+ return std::make_pair(dof, Key{entityIndex, dim-1, index});
+ }
+
+ offset += 4 * (order[0]-1 + order[1]-1 + order[2]-1);
+ if (nbdy == 1) // Face DOF
+ {
+ Key faceKey;
+ if (ibdy) // On i-normal face
+ {
+ dof = (j - 1) + ((order[1]-1) * (k - 1)) + (i ? (order[1]-1) * (order[2]-1) : 0) + offset;
+ entityIndex = (i ? 1 : 0);
+ faceKey = VtkLagrangePointSetBuilder<K,2>::calcQuadKey(j-1,k-1,{order[1]-2, order[2]-2}).second;
+ }
+ else {
+ offset += 2 * (order[1] - 1) * (order[2] - 1);
+ if (jbdy) // On j-normal face
+ {
+ dof = (i - 1) + ((order[0]-1) * (k - 1)) + (j ? (order[2]-1) * (order[0]-1) : 0) + offset;
+ entityIndex = (j ? 3 : 2);
+ faceKey = VtkLagrangePointSetBuilder<K,2>::calcQuadKey(i-1,k-1,{order[0]-2, order[2]-2}).second;
+ }
+ else
+ { // kbdy, On k-normal face
+ offset += 2 * (order[2]-1) * (order[0]-1);
+ dof = (i - 1) + ((order[0]-1) * (j - 1)) + (k ? (order[0]-1) * (order[1]-1) : 0) + offset;
+ entityIndex = (k ? 5 : 4);
+ faceKey = VtkLagrangePointSetBuilder<K,2>::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second;
+ }
+ }
+ return std::make_pair(dof, Key{entityIndex, dim-2, faceKey.index()});
+ }
+
+ // nbdy == 0: Body DOF
+ offset += 2 * ((order[1]-1) * (order[2]-1) + (order[2]-1) * (order[0]-1) + (order[0]-1) * (order[1]-1));
+ dof = offset + (i - 1) + (order[0]-1) * ((j - 1) + (order[1]-1) * ((k - 1)));
+ Key innerKey = VtkLagrangePointSetBuilder<K,3>::calcHexKey(i-1,j-1,k-1,{order[0]-2, order[1]-2, order[2]-2}).second;
+ return std::make_pair(dof, Key{0, dim-3, innerKey.index()});
+ }
+};
+
+}} // end namespace Dune::Impl
diff --git a/dune/vtk/utility/test/CMakeLists.txt b/dune/vtk/utility/test/CMakeLists.txt
new file mode 100644
index 0000000..f752b18
--- /dev/null
+++ b/dune/vtk/utility/test/CMakeLists.txt
@@ -0,0 +1,13 @@
+dune_add_test(NAME test-lagrange1
+ SOURCES test-lagrange.cc
+ COMPILE_DEFINITIONS "CHECKDIM=1")
+
+dune_add_test(NAME test-lagrange2
+ SOURCES test-lagrange.cc
+ COMPILE_DEFINITIONS "CHECKDIM=2")
+
+dune_add_test(NAME test-lagrange3
+ SOURCES test-lagrange.cc
+ COMPILE_DEFINITIONS "CHECKDIM=3")
+
+dune_add_test(SOURCES test-lfecache.cc)
\ No newline at end of file
diff --git a/dune/vtk/utility/test/test-lagrange.cc b/dune/vtk/utility/test/test-lagrange.cc
new file mode 100644
index 0000000..62ed50a
--- /dev/null
+++ b/dune/vtk/utility/test/test-lagrange.cc
@@ -0,0 +1,169 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <config.h>
+
+#include <dune/localfunctions/utility/field.hh>
+#include <dune/localfunctions/utility/basisprint.hh>
+
+#include <dune/localfunctions/lagrange/lagrangebasis.hh>
+#include <dune/vtk/utility/lagrangepoints.hh>
+
+/**
+ * \file
+ * \brief Performs some tests for the generic Lagrange
+ * shape functions on simplices.
+ *
+ * The topology can be chosen at compile time by setting TOPOLOGY
+ * to a string like
+ * \code
+ * Pyramid<Pyramid<Point> > >
+ * \endcode
+ * which generates a 2d simplex. If TOPOLOGY is not set, all
+ * topologies up to 4d are tested. Note, this may lead to prolonged
+ * compiler runs.
+ *
+ * For debugging purpuse the functions and the derivatives can be
+ * printed. You have to define the macro TEST_OUTPUT_FUNCTIONS to
+ * activate this function.
+ */
+
+#if HAVE_GMP
+typedef Dune::GMPField< 128 > StorageField;
+typedef Dune::GMPField< 512 > ComputeField;
+#else
+typedef double StorageField;
+typedef double ComputeField;
+#endif
+
+template <class Basis,class Points>
+bool test(const Basis &basis, const Points &points, bool verbose)
+{
+ bool ret = true;
+ std::vector< Dune::FieldVector< double, 1 > > y( basis.size() );
+ for( unsigned int index = 0; index < points.size(); ++index )
+ {
+ if (verbose)
+ std::cout << index << " " << points[ index ].point() << " "
+ << points[ index ].localKey()
+ << std::endl;
+ basis.evaluate( points[ index ].point(), y );
+ bool first = true;
+ for( unsigned int i = 0; i < y.size(); ++i )
+ {
+ if( fabs( y[ i ] - double( i == index ) ) > 1e-10 )
+ {
+ if (first) {
+ std::cout << "ERROR: "
+ << index << " -> "
+ << "x = " << points[ index ].point()
+ << " (codim = " << points[ index ].localKey().codim() << ", "
+ << "subentity = " << points[ index ].localKey().subEntity() << ", "
+ << "index = " << points[ index ].localKey().index() << "):" << std::endl;
+ first = false;
+ }
+ if (1)
+ std::cout << " y[ " << i << " ] = " << y[ i ] << " "
+ << " error : " << fabs( y[ i ] - double( i == index ) )
+ << std::endl;
+ ret = false;
+ }
+ }
+ }
+ return ret;
+}
+
+template <class Topology>
+bool test(unsigned int order, bool verbose = false)
+{
+ typedef Dune::LagrangeBasisFactory<Dune::VtkLagrangePointSet,Topology::dimension,StorageField,ComputeField> BasisFactory;
+ typedef Dune::LagrangeCoefficientsFactory< Dune::VtkLagrangePointSet, Topology::dimension,double > LagrangeCoefficientsFactory;
+
+ bool ret = true;
+
+ for (unsigned int o = 0; o <= order; ++o)
+ {
+ const typename LagrangeCoefficientsFactory::Object *pointsPtr = LagrangeCoefficientsFactory::template create< Topology >( o );
+
+ if ( pointsPtr == 0)
+ continue;
+
+ std::cout << "# Testing " << Topology::name() << " in dimension " << Topology::dimension << " with order " << o << std::endl;
+
+ typename BasisFactory::Object &basis = *BasisFactory::template create<Topology>(o);
+
+ ret |= test(basis,*pointsPtr,verbose);
+
+ // define the macro TEST_OUTPUT_FUNCTIONS to output files containing functions and
+ // derivatives in a human readabible form (aka LaTeX source)
+#ifdef TEST_OUTPUT_FUNCTIONS
+ std::stringstream name;
+ name << "lagrange_" << Topology::name() << "_p" << o << ".basis";
+ std::ofstream out(name.str().c_str());
+ Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField>(out,basis);
+ Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField>(out,basis);
+#endif // TEST_OUTPUT_FUNCTIONS
+
+ LagrangeCoefficientsFactory::release( pointsPtr );
+ BasisFactory::release( &basis );
+ }
+
+ if (verbose)
+ std::cout << std::endl << std::endl << std::endl;
+ if (!ret) {
+ std::cout << " FAILED !" << std::endl;
+ }
+ return ret;
+}
+
+#ifdef CHECKDIM
+ #if CHECKDIM==1
+ #define CHECKDIM1
+ #elif CHECKDIM==2
+ #define CHECKDIM2
+ #elif CHECKDIM==3
+ #define CHECKDIM3
+ #endif
+#else
+ #define CHECKDIM1
+ #define CHECKDIM2
+ #define CHECKDIM3
+#endif
+
+
+
+int main ( int argc, char **argv )
+{
+ using namespace Dune;
+ using namespace Impl;
+
+ unsigned int order = (argc < 2) ? 5 : atoi(argv[1]);
+
+ if (argc < 2)
+ {
+ std::cerr << "Usage: " << argv[ 0 ] << " <p>" << std::endl
+ << "Using default order of " << order << std::endl;
+ }
+#ifdef TOPOLOGY
+ return (test<TOPOLOGY>(order) ? 0 : 1 );
+#else
+ bool tests = true;
+
+#ifdef CHECKDIM1
+ tests &= test<Prism<Point> > (order); // line
+ tests &= test<Pyramid<Point> > (order); // line
+#endif
+
+#ifdef CHECKDIM2
+ tests &= test<Prism<Prism<Point> > > (order); // quad
+ tests &= test<Pyramid<Pyramid<Point> > >(order); // triangle
+#endif
+
+#ifdef CHECKDIM3
+ tests &= test<Prism<Prism<Prism<Point> > > >(order); // hexahedron
+ // tests &= test<Prism<Pyramid<Pyramid<Point> > > >(order);
+ tests &= test<Pyramid<Pyramid<Pyramid<Point> > > >(order); // tetrahedron
+#endif
+
+ return (tests ? 0 : 1);
+#endif // TOPOLOGY
+}
diff --git a/dune/vtk/utility/test/test-lfecache.cc b/dune/vtk/utility/test/test-lfecache.cc
new file mode 100644
index 0000000..2782090
--- /dev/null
+++ b/dune/vtk/utility/test/test-lfecache.cc
@@ -0,0 +1,56 @@
+// -*- 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 <dune/common/hybridutilities.hh>
+#include <dune/common/rangeutilities.hh>
+#include <dune/common/std/utility.hh>
+#include <dune/common/test/testsuite.hh>
+
+#include <dune/geometry/type.hh>
+
+#include <dune/vtk/utility/lagrangelfecache.hh>
+
+template <class FiniteElementCache>
+int test (FiniteElementCache& cache, Dune::GeometryType type)
+{
+ using FiniteElement = typename FiniteElementCache::FiniteElementType;
+ const FiniteElement& finiteElement = cache.get(type);
+
+ if (!(finiteElement.type() == type))
+ DUNE_THROW(Dune::Exception, "Wrong GeometryType in constructed local finite-element");
+
+ return finiteElement.localBasis().order();
+}
+
+int main()
+{
+ using namespace Dune;
+ const int max_dim = 3;
+ const unsigned int max_order = 4;
+
+ TestSuite testSuite;
+ Hybrid::forEach(StaticIntegralRange<int,max_dim+1,1>{},[&](auto dim)
+ {
+ // template order parameter
+ Hybrid::forEach(StaticIntegralRange<int,max_order+1,1>{},[&](auto k) {
+ using FiniteElementCache = LagrangeLFECache<double, double, dim.value, k.value>;
+ FiniteElementCache cache{};
+ testSuite.check(test(cache, Dune::GeometryTypes::simplex(dim)) == k);
+ testSuite.check(test(cache, Dune::GeometryTypes::cube(dim)) == k);
+ });
+
+ // runtime order parameter
+ for (unsigned int k = 1; k <= max_order; ++k) {
+ using FiniteElementCache = LagrangeLFECache<double, double, dim>;
+ FiniteElementCache cache{k};
+ testSuite.check(test(cache, Dune::GeometryTypes::simplex(dim)) == k);
+ testSuite.check(test(cache, Dune::GeometryTypes::cube(dim)) == k);
+ }
+ });
+
+ return testSuite.exit();
+}
diff --git a/dune/vtk/vtkfunction.hh b/dune/vtk/vtkfunction.hh
index d745da6..48444ec 100644
--- a/dune/vtk/vtkfunction.hh
+++ b/dune/vtk/vtkfunction.hh
@@ -22,12 +22,12 @@ namespace Dune
class VtkFunction
{
template <class F>
- using HasLocalFunction = decltype(localFunction(std::declval<F>()));
+ using LocalFunction = decltype(localFunction(std::declval<F>()));
using Domain = typename GridView::template Codim<0>::Entity::Geometry::LocalCoordinate;
template <class F>
- using Range = std::decay_t<decltype(std::declval<F>()(std::declval<Domain>()))>;
+ using Range = std::decay_t<std::result_of_t<F(Domain)>>;
private:
@@ -71,14 +71,14 @@ namespace Dune
* GridFunction if not given.
**/
template <class F,
- class = void_t<HasLocalFunction<F>> >
+ class = void_t<LocalFunction<F>> >
VtkFunction (F&& fct, std::string name,
Std::optional<int> ncomps = {},
Std::optional<Vtk::DataTypes> type = {})
: localFct_(localFunction(std::forward<F>(fct)))
, name_(std::move(name))
{
- using R = Range<decltype(localFunction(std::forward<F>(fct)))>;
+ using R = Range<LocalFunction<F>>;
ncomps_ = ncomps ? *ncomps : sizeOf<R>();
type_ = type ? *type : Vtk::Map::type<R>();
@@ -86,7 +86,7 @@ namespace Dune
/// Constructor that forward the number of components and data type to the other constructor
template <class F,
- class = void_t<HasLocalFunction<F>> >
+ class = void_t<LocalFunction<F>> >
VtkFunction (F&& fct, Vtk::FieldInfo fieldInfo,
Std::optional<Vtk::DataTypes> type = {})
: VtkFunction(std::forward<F>(fct), fieldInfo.name(), fieldInfo.ncomps(), type)
diff --git a/dune/vtk/vtkreader.hh b/dune/vtk/vtkreader.hh
index c92bc82..294ced1 100644
--- a/dune/vtk/vtkreader.hh
+++ b/dune/vtk/vtkreader.hh
@@ -41,13 +41,12 @@ namespace Dune
std::uint64_t offset = 0;
};
- using Entity = typename Grid::template Codim<0>::Entity;
- using GlobalCoordinate = typename Entity::Geometry::GlobalCoordinate;
+ using GlobalCoordinate = typename GridCreator::GlobalCoordinate;
public:
/// Constructor. Creates a new GridCreator with the passed factory
template <class... Args,
- disableCopyMove<VtkReader, Args...> = 0>
+ std::enable_if_t<std::is_constructible<GridCreator, Args...>::value,int> = 0>
explicit VtkReader (Args&&... args)
: VtkReader(std::make_shared<GridCreator>(std::forward<Args>(args)...))
{}
@@ -76,6 +75,12 @@ namespace Dune
reader.readFromFile(filename);
}
+ /// Obtains the creator of the grid
+ GridCreator& gridCreator ()
+ {
+ return *creator_;
+ }
+
/// Advanced read methods
/// @{
diff --git a/dune/vtk/vtkreader.impl.hh b/dune/vtk/vtkreader.impl.hh
index 9427cde..7acc338 100644
--- a/dune/vtk/vtkreader.impl.hh
+++ b/dune/vtk/vtkreader.impl.hh
@@ -292,7 +292,7 @@ template <class IStream, class T, class Sections>
Sections readDataArray (IStream& input, std::vector<T>& values, std::size_t max_size,
Sections section, Sections parent_section)
{
- values.reserve(max_size);
+ values.reserve(max_size < std::size_t(-1) ? max_size : 0);
using S = std::conditional_t<(sizeof(T) <= 1), std::uint16_t, T>; // problem when reading chars as ints
std::size_t idx = 0;
@@ -300,11 +300,27 @@ Sections readDataArray (IStream& input, std::vector<T>& values, std::size_t max_
trim(line);
if (line.substr(1,10) == "/DataArray")
return parent_section;
+ if (line[0] == '<')
+ break;
std::istringstream stream(line);
S value;
for (; stream >> value; idx++)
values.push_back(T(value));
+ if (idx >= max_size)
+ break;
+ }
+
+ return section;
+}
+
+template <class IStream, class Sections>
+Sections skipRestOfDataArray (IStream& input, Sections section, Sections parent_section)
+{
+ for (std::string line; std::getline(input, line);) {
+ ltrim(line);
+ if (line.substr(1,10) == "/DataArray")
+ return parent_section;
}
return section;
@@ -324,6 +340,8 @@ VtkReader<Grid,Creator>::readPoints (std::ifstream& input, std::string name)
std::vector<T> point_values;
sec = readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS);
+ if (sec != POINTS)
+ sec = skipRestOfDataArray(input, POINTS_DATA_ARRAY, POINTS);
assert(sec == POINTS);
assert(point_values.size() == 3*numberOfPoints_);
@@ -385,7 +403,7 @@ VtkReader<Grid,Creator>::readCells (std::ifstream& input, std::string name)
max_size = vec_offsets.back();
else
max_size = numberOfCells_ * max_vertices;
- sec = readDataArray(input, vec_connectivity, max_size, CELLS_DATA_ARRAY, CELLS);
+ sec = readDataArray(input, vec_connectivity, std::size_t(-1), CELLS_DATA_ARRAY, CELLS);
} else if (name == "global_point_ids") {
sec = readDataArray(input, vec_point_ids, numberOfPoints_, CELLS_DATA_ARRAY, CELLS);
assert(vec_point_ids.size() == numberOfPoints_);
diff --git a/dune/vtk/vtktypes.cc b/dune/vtk/vtktypes.cc
index 3ea7b9a..90f993e 100644
--- a/dune/vtk/vtktypes.cc
+++ b/dune/vtk/vtktypes.cc
@@ -50,6 +50,21 @@ GeometryType to_geometry (std::uint8_t cell)
case HEXAHEDRON: return GeometryTypes::hexahedron;
case WEDGE: return GeometryTypes::prism;
case PYRAMID: return GeometryTypes::pyramid;
+
+ // Quadratic VTK cell types
+ case QUADRATIC_EDGE: return GeometryTypes::line;
+ case QUADRATIC_TRIANGLE: return GeometryTypes::triangle;
+ case QUADRATIC_QUAD: return GeometryTypes::quadrilateral;
+ case QUADRATIC_TETRA: return GeometryTypes::tetrahedron;
+ case QUADRATIC_HEXAHEDRON: return GeometryTypes::hexahedron;
+
+ // Arbitrary order Lagrange elements
+ case LAGRANGE_CURVE: return GeometryTypes::line;
+ case LAGRANGE_TRIANGLE: return GeometryTypes::triangle;
+ case LAGRANGE_QUADRILATERAL:return GeometryTypes::quadrilateral;
+ case LAGRANGE_TETRAHEDRON: return GeometryTypes::tetrahedron;
+ case LAGRANGE_HEXAHEDRON: return GeometryTypes::hexahedron;
+ case LAGRANGE_WEDGE: return GeometryTypes::prism;
default:
DUNE_THROW(RangeError, "CellType does not map to GeometryType.");
std::abort();
@@ -153,6 +168,32 @@ CellType::CellType (GeometryType const& t, CellParametrization parametrization)
std::cerr << "Geometry Type not supported by VTK!\n";
std::abort();
}
+ } else if (parametrization == LAGRANGE) {
+ if (t.isLine()) {
+ type_ = LAGRANGE_CURVE;
+ }
+ else if (t.isTriangle()) {
+ type_ = LAGRANGE_TRIANGLE;
+ }
+ else if (t.isQuadrilateral()) {
+ type_ = LAGRANGE_QUADRILATERAL;
+ }
+ else if (t.isTetrahedron()) {
+ type_ = LAGRANGE_TETRAHEDRON;
+ }
+ else if (t.isHexahedron()) {
+ type_ = LAGRANGE_HEXAHEDRON;
+ }
+ else if (t.isPrism()) {
+ type_ = LAGRANGE_WEDGE;
+ }
+ else if (t.isPyramid()) {
+ type_ = LAGRANGE_PYRAMID;
+ }
+ else {
+ std::cerr << "Geometry Type not supported by VTK!\n";
+ std::abort();
+ }
}
}
diff --git a/dune/vtk/vtktypes.hh b/dune/vtk/vtktypes.hh
index 3d41852..67f0c49 100644
--- a/dune/vtk/vtktypes.hh
+++ b/dune/vtk/vtktypes.hh
@@ -33,7 +33,8 @@ namespace Dune
enum CellParametrization {
LINEAR,
- QUADRATIC
+ QUADRATIC,
+ LAGRANGE
};
enum CellTypes : std::uint8_t {
@@ -57,7 +58,15 @@ namespace Dune
QUADRATIC_TRIANGLE = 22,
QUADRATIC_QUAD = 23,
QUADRATIC_TETRA = 24,
- QUADRATIC_HEXAHEDRON = 25
+ QUADRATIC_HEXAHEDRON = 25,
+ // Arbitrary order Lagrange elements
+ LAGRANGE_CURVE = 68,
+ LAGRANGE_TRIANGLE = 69,
+ LAGRANGE_QUADRILATERAL = 70,
+ LAGRANGE_TETRAHEDRON = 71,
+ LAGRANGE_HEXAHEDRON = 72,
+ LAGRANGE_WEDGE = 73,
+ LAGRANGE_PYRAMID = 74,
};
GeometryType to_geometry (std::uint8_t);
@@ -110,7 +119,7 @@ namespace Dune
private:
std::uint8_t type_;
std::vector<int> permutation_;
- bool noPermutation_;
+ bool noPermutation_ = true;
};
diff --git a/dune/vtk/vtkwriterinterface.hh b/dune/vtk/vtkwriterinterface.hh
index 877d760..33121d2 100644
--- a/dune/vtk/vtkwriterinterface.hh
+++ b/dune/vtk/vtkwriterinterface.hh
@@ -92,9 +92,9 @@ namespace Dune
* Constructor for possible arguments.
**/
template <class Function, class... Args>
- VtkWriterInterface& addPointData (Function const& fct, Args&&... args)
+ VtkWriterInterface& addPointData (Function&& fct, Args&&... args)
{
- pointData_.emplace_back(fct, std::forward<Args>(args)...);
+ pointData_.emplace_back(std::forward<Function>(fct), std::forward<Args>(args)...);
return *this;
}
@@ -106,9 +106,9 @@ namespace Dune
* See \ref VtkFunction Constructor for possible arguments.
**/
template <class Function, class... Args>
- VtkWriterInterface& addCellData (Function const& fct, Args&&... args)
+ VtkWriterInterface& addCellData (Function&& fct, Args&&... args)
{
- cellData_.emplace_back(fct, std::forward<Args>(args)...);
+ cellData_.emplace_back(std::forward<Function>(fct), std::forward<Args>(args)...);
return *this;
}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0ec1e3d..0c09ee2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,3 +1,5 @@
+set(GRID_PATH "${CMAKE_CURRENT_SOURCE_DIR}/../doc/")
+
dune_add_test(SOURCES vtkreader.cc
LINK_LIBRARIES dunevtk)
@@ -17,6 +19,14 @@ dune_add_test(SOURCES datacollector.cc
LINK_LIBRARIES dunevtk
CMAKE_GUARD dune-functions_FOUND)
+dune_add_test(SOURCES lagrangepoints.cc
+ LINK_LIBRARIES dunevtk)
+
+dune_add_test(SOURCES lagrangereader.cc
+ COMPILE_DEFINITIONS "GRID_PATH=\"${GRID_PATH}\""
+ LINK_LIBRARIES dunevtk
+ CMAKE_GUARD dune-foamgrid_FOUND)
+
dune_add_test(SOURCES structuredgridwriter.cc
LINK_LIBRARIES dunevtk
CMAKE_GUARD dune-functions_FOUND)
@@ -37,7 +47,6 @@ dune_add_test(SOURCES vectorwriter.cc
LINK_LIBRARIES dunevtk
CMAKE_GUARD dune-functions_FOUND)
-
if (dune-polygongrid_FOUND)
# CMAKE_GUARD can not be used, since a dummy target is created and linked against dunepolygongrid
dune_add_test(SOURCES polygongrid.cc
diff --git a/src/datacollector.cc b/src/datacollector.cc
index 607591e..9f126d2 100644
--- a/src/datacollector.cc
+++ b/src/datacollector.cc
@@ -19,13 +19,19 @@
#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
#include <dune/vtk/datacollectors/continuousdatacollector.hh>
#include <dune/vtk/datacollectors/discontinuousdatacollector.hh>
#include <dune/vtk/datacollectors/quadraticdatacollector.hh>
+#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
using namespace Dune;
using namespace Dune::Functions;
@@ -45,6 +51,7 @@ void write_dc (std::string prefix, GridView const& gridView, Fct1 const& fct1, F
template <class GridView>
void write (std::string prefix, GridView const& gridView)
{
+ std::cout << prefix << "..." << std::endl;
#if ! DUNE_VERSION_NEWER(DUNE_FUNCTIONS, 2, 6)
using namespace BasisBuilder;
#else
@@ -69,6 +76,10 @@ void write (std::string prefix, GridView const& gridView)
write_dc<ContinuousDataCollector<GridView>>(prefix + "_continuous", gridView, p1Interpol, p1Analytic);
write_dc<DiscontinuousDataCollector<GridView>>(prefix + "_discontinuous", gridView, p1Interpol, p1Analytic);
write_dc<QuadraticDataCollector<GridView>>(prefix + "_quadratic", gridView, p1Interpol, p1Analytic);
+
+ Hybrid::forEach(StaticIntegralRange<int,7,1>{}, [&](auto p) {
+ write_dc<LagrangeDataCollector<GridView,p>>(prefix + "_lagrange_p" + std::to_string(p), gridView, p1Interpol, p1Analytic);
+ });
}
template <int I>
@@ -78,12 +89,24 @@ int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
- Hybrid::forEach(std::make_tuple(int_<1>{}, int_<2>{}, int_<3>{}), [](auto dim)
+ Hybrid::forEach(StaticIntegralRange<int,4,1>{}, [](auto dim)
{
using GridType = YaspGrid<dim.value>;
FieldVector<double,dim.value> upperRight; upperRight = 1.0;
- auto numElements = filledArray<dim.value,int>(8);
+ auto numElements = filledArray<dim.value,int>(2);
GridType grid(upperRight, numElements, 0, 0);
write("datacollector_yasp", grid.leafGridView());
});
+
+#if HAVE_UG
+ Hybrid::forEach(StaticIntegralRange<int,4,2>{}, [](auto dim)
+ {
+ using GridType = UGGrid<dim.value>;
+ FieldVector<double,dim.value> lowerLeft; lowerLeft = 0.0;
+ FieldVector<double,dim.value> upperRight; upperRight = 1.0;
+ auto numElements = filledArray<dim.value,unsigned int>(4);
+ auto grid = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
+ write("datacollector_ug", grid->leafGridView());
+ });
+#endif
}
diff --git a/src/lagrangepoints.cc b/src/lagrangepoints.cc
new file mode 100644
index 0000000..a40f83e
--- /dev/null
+++ b/src/lagrangepoints.cc
@@ -0,0 +1,68 @@
+// -*- 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 <iostream>
+#include <vector>
+
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/indices.hh>
+
+#include <dune/vtk/utility/lagrangepoints.hh>
+#include <dune/localfunctions/lagrange/interpolation.hh>
+#include <dune/localfunctions/lagrange/lagrangebasis.hh>
+#include <dune/localfunctions/lagrange/lagrangecoefficients.hh>
+#include <dune/localfunctions/utility/localfiniteelement.hh>
+
+using namespace Dune;
+
+template <std::size_t dim>
+void write (std::string prefix, index_constant<dim>)
+{
+ for (int order = 1; order < 6; ++order) {
+ std::cout << "order: " << order << std::endl;
+ {
+ VtkLagrangePointSet<double, dim> pointSet(order);
+ pointSet.build(GeometryTypes::cube(dim));
+
+ std::size_t i = 0;
+ std::cout << "Cube:" << GeometryTypes::cube(dim) << std::endl;
+ for (auto const& p : pointSet) {
+ std::cout << i++ << ") p = " << p.point() << ", key = " << p.localKey() << std::endl;
+ }
+
+ using BasisF = LagrangeBasisFactory<VtkLagrangePointSet, dim, double, double>;
+ using CoefficientF = LagrangeCoefficientsFactory<VtkLagrangePointSet, dim, double>;
+ using InterpolationF = LagrangeInterpolationFactory<VtkLagrangePointSet, dim, double>;
+ GenericLocalFiniteElement<BasisF, CoefficientF, InterpolationF> localFE(GeometryTypes::cube(dim), order);
+
+ auto const& localBasis = localFE.localBasis();
+ auto const& localCoefficints = localFE.localCoefficients();
+ auto const& localInterpolation = localFE.localInterpolation();
+ }
+
+ {
+ VtkLagrangePointSet<double, dim> pointSet(order);
+ pointSet.build(GeometryTypes::simplex(dim));
+
+ std::size_t i = 0;
+ std::cout << "Simplex:" << GeometryTypes::simplex(dim) << std::endl;
+ for (auto const& p : pointSet) {
+ std::cout << i++ << ") p = " << p.point() << ", key = " << p.localKey() << std::endl;
+ }
+ }
+ }
+}
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ Hybrid::forEach(std::make_tuple(index_constant<2>{}), [](auto dim)
+ {
+ write("lagrangepoints", dim);
+ });
+}
diff --git a/src/lagrangereader.cc b/src/lagrangereader.cc
new file mode 100644
index 0000000..8a3543c
--- /dev/null
+++ b/src/lagrangereader.cc
@@ -0,0 +1,95 @@
+// -*- 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 <iostream>
+#include <vector>
+
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/exceptions.hh> // We use exceptions
+#include <dune/common/filledarray.hh>
+#include <dune/common/test/testsuite.hh>
+
+// #include <dune/grid/uggrid.hh>
+// #include <dune/alugrid/grid.hh>
+#include <dune/foamgrid/foamgrid.hh>
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/vtk/vtkreader.hh>
+#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
+#include <dune/vtk/gridcreators/lagrangegridcreator.hh>
+#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
+
+using namespace Dune;
+
+#ifndef GRID_PATH
+#define GRID_PATH
+#endif
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ const int dim = 2;
+ const int dow = 3;
+ const int order = 3;
+
+ // using GridType = UGGrid<dim>;
+ // using GridType = Dune::ALUGrid<dim,dow,Dune::simplex,Dune::conforming>;
+ using GridType = FoamGrid<dim,dow>;
+ using GridView = typename GridType::LeafGridView;
+ using DataCollector = LagrangeDataCollector<GridView, order>;
+ using GridCreator = LagrangeGridCreator<GridType>;
+
+ std::string filename = "triangles_" + std::to_string(dow) + "d_order" + std::to_string(order);
+
+ TestSuite testSuite;
+
+ { // Test using the (static) file-reader interface
+ std::cout << "Test 1..." << std::endl;
+ auto gridPtr = VtkReader<GridType,GridCreator>::read(GRID_PATH "/" + filename + ".vtu");
+ auto& grid = *gridPtr;
+ grid.globalRefine(2);
+
+ VtkUnstructuredGridWriter<GridView,DataCollector> vtkWriter(grid.leafGridView(), Vtk::ASCII);
+ vtkWriter.write(filename + "_out1.vtu");
+ }
+
+ { // Test using an instantiated reader
+ std::cout << "Test 2..." << std::endl;
+ GridFactory<GridType> factory;
+ VtkReader<GridType,GridCreator> reader(factory);
+ reader.readFromFile(GRID_PATH "/" + filename + ".vtu");
+ auto gridPtr = factory.createGrid();
+ auto& grid = *gridPtr;
+
+ auto&& creator = reader.gridCreator();
+ testSuite.check(creator.order() == order, "order");
+
+ VtkUnstructuredGridWriter<GridView,DataCollector> vtkWriter(grid.leafGridView(), Vtk::ASCII);
+ vtkWriter.addPointData(creator, "param", 3);
+ vtkWriter.write(filename + "_out2.vtu");
+ }
+
+ { // Test using an explicit grid-creator
+ std::cout << "Test 3..." << std::endl;
+ GridFactory<GridType> factory;
+ GridCreator creator(factory);
+ VtkReader<GridType,GridCreator> reader(creator);
+ reader.readFromFile(GRID_PATH "/" + filename + ".vtu");
+ auto gridPtr = factory.createGrid();
+ auto& grid = *gridPtr;
+
+ testSuite.check(creator.order() == order, "order");
+
+ VtkUnstructuredGridWriter<GridView,DataCollector> vtkWriter(grid.leafGridView(), Vtk::ASCII);
+ vtkWriter.addPointData(creator, "param", 3);
+ vtkWriter.write(filename + "_out3.vtu");
+ }
+
+ return testSuite.exit();
+}
--
GitLab