implementation of hex and tetra lagrange point sets

dune/vtk/datacollectors/lagrangedatacollector.hh

@@ -91,7 +91,7 @@ public:
 
         auto const& localKey = p.localKey();
         std::size_t idx = 3 * (localKey.codim() == dim
-            ? indexSet.subIndex(element, localKey.subEntity(), localKey.codim())
+            ? indexSet.subIndex(element, localKey.subEntity(), dim)
             : innerDOFs*indexSet.index(element) + (i - vertexDOFs) + shift);
 
         auto v = geometry.global(p.point());
@@ -138,7 +138,7 @@ public:
         auto const& p = pointSet[i];
         auto const& localKey = p.localKey();
         std::size_t idx = (localKey.codim() == dim
-            ? indexSet.subIndex(element, localKey.subEntity(), localKey.codim())
+            ? indexSet.subIndex(element, localKey.subEntity(), dim)
             : innerDOFs*indexSet.index(element) + (i - vertexDOFs) + shift);
         cells.connectivity.push_back(idx);
       }
@@ -172,7 +172,7 @@ public:
         auto const& p = pointSet[i];
         auto const& localKey = p.localKey();
         std::size_t idx = nComps * (localKey.codim() == dim
-            ? indexSet.subIndex(element, localKey.subEntity(), localKey.codim())
+            ? indexSet.subIndex(element, localKey.subEntity(), dim)
             : innerDOFs*indexSet.index(element) + (i - vertexDOFs) + shift);
 
         for (int comp = 0; comp < nComps; ++comp)

src/datacollector.cc

@@ -19,7 +19,12 @@
 #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>
 
@@ -46,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
@@ -70,8 +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);
-  if (GridView::dimension <= 2)
-    write_dc<LagrangeDataCollector<GridView,4>>(prefix + "_lagrange", 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>
@@ -81,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>(4);
+    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
 }