[spaces.h1.continuous-lagrange] allow for range dimension r > 1

.gitsuper

@@ -4,7 +4,7 @@ status = 1a3bcab04b011a5d6e44f9983cae6ff89fa695e8 bin (heads/master)
 	 28c9ce81c14c878a71e907ab05b9bb72df77883e config.opts (heads/master)
 	 f308c6637edd65dcb83c4c1a46feaf05b958130e dune-alugrid (v2.6.0-7-gf308c663)
 	 76d7f0c9886a061571cb8dc66dd45a4ef86e7a58 dune-common (v2.2.1-2269-g76d7f0c9)
-	+019022cfff5b13501f0386237c20483d4bb6610c dune-gdt (heads/saddlepoint)
+	+4b461b6c6e04d68f6b67677186ac9275a043a25f dune-gdt (heads/saddlepoint)
 	 5235397bc16d24c759a1672fed7b8cfde4852e52 dune-geometry (v2.2.0-834-g5235397)
 	 af5766f0df47e3d0b62ea486efb9cdbf8e1cfc52 dune-grid (v2.2.0-2671-gaf5766f0d)
 	 1369ae9329d0928480d6b18ed772fc77e1abf752 dune-grid-glue (v2.4.0-161-g1369ae9)
@@ -44,7 +44,7 @@ commit = 76d7f0c9886a061571cb8dc66dd45a4ef86e7a58
 [submodule.dune-gdt]
 remote = git@github.com:dune-community/dune-gdt.git
 status = c0b1735fab0ecbd4bb4f1eaa27cb65fe813e98f0 .vcsetup (heads/master)
-commit = 019022cfff5b13501f0386237c20483d4bb6610c
+commit = 4b461b6c6e04d68f6b67677186ac9275a043a25f
 
 [submodule.dune-geometry]
 remote = git@github.com:dune-community/dune-geometry.git

dune/gdt/spaces/h1/continuous-lagrange.hh

@@ -51,11 +51,11 @@ namespace GDT {
  *
  * \sa make_local_lagrange_finite_element
  */
-template <class GV, class R = double>
-class ContinuousLagrangeSpace : public SpaceInterface<GV, 1, 1, R>
+template <class GV, size_t r = 1, class R = double>
+class ContinuousLagrangeSpace : public SpaceInterface<GV, r, 1, R>
 {
-  using ThisType = ContinuousLagrangeSpace<GV, R>;
-  using BaseType = SpaceInterface<GV, 1, 1, R>;
+  using ThisType = ContinuousLagrangeSpace;
+  using BaseType = SpaceInterface<GV, r, 1, R>;
 
 public:
   using BaseType::d;
@@ -67,7 +67,7 @@ public:
 
 private:
   using MapperImplementation = ContinuousMapper<GridViewType, FiniteElementType>;
-  using GlobalBasisImplementation = DefaultGlobalBasis<GridViewType, 1, 1, R>;
+  using GlobalBasisImplementation = DefaultGlobalBasis<GridViewType, r, 1, R>;
 
 public:
   ContinuousLagrangeSpace(GridViewType grd_vw, const int order)
@@ -87,7 +87,7 @@ public:
     // create finite elements
     for (auto&& geometry_type : grid_view_.indexSet().types(0))
       finite_elements_->insert(
-          std::make_pair(geometry_type, make_local_lagrange_finite_element<D, d, R>(geometry_type, order)));
+          std::make_pair(geometry_type, make_local_lagrange_finite_element<D, d, R, r>(geometry_type, order)));
     // check
     if (d == 3 && finite_elements_->size() != 1)
       DUNE_THROW(Exceptions::space_error,
@@ -169,13 +169,23 @@ private:
 }; // class ContinuousLagrangeSpace
 
 
+/**
+ * \sa ContinuousLagrangeSpace
+ */
+template <size_t r, class GV, class R = double>
+ContinuousLagrangeSpace<GridView<GV>, r, R> make_continuous_lagrange_space(GridView<GV> grid_view, const int order)
+{
+  return ContinuousLagrangeSpace<GridView<GV>, r, R>(grid_view, order);
+}
+
+
 /**
  * \sa ContinuousLagrangeSpace
  */
 template <class GV, class R = double>
-ContinuousLagrangeSpace<GridView<GV>, R> make_continuous_lagrange_space(GridView<GV> grid_view, const int order)
+ContinuousLagrangeSpace<GridView<GV>, 1, R> make_continuous_lagrange_space(GridView<GV> grid_view, const int order)
 {
-  return ContinuousLagrangeSpace<GridView<GV>, R>(grid_view, order);
+  return ContinuousLagrangeSpace<GridView<GV>, 1, R>(grid_view, order);
 }
 
 

dune/gdt/spaces/mapper/continuous.hh

@@ -35,6 +35,7 @@ class ContinuousMapper : public MapperInterface<GV>
   static_assert(is_local_finite_element<FiniteElement>::value, "");
   using ThisType = ContinuousMapper<GV, FiniteElement>;
   using BaseType = MapperInterface<GV>;
+  static const size_t r = FiniteElement::r;
 
   template <int d>
   struct GeometryTypeLayout
@@ -86,9 +87,11 @@ public:
       const auto& coeffs = finite_element.coefficients();
       for (size_t ii = 0; ii < coeffs.size(); ++ii) {
         const auto& local_key = coeffs.local_key(ii);
-        DUNE_THROW_IF(local_key.index() != 0, // Would require twisting of DoFs and possibly more knowledge from the FE
+        // Currently only works if each subEntity has exactly r DoFs, if there is a variable number of DoFs per element
+        // we would need to do more complicated things in the global index mapping.
+        DUNE_THROW_IF(!(local_key.index() < r),
                       Exceptions::mapper_error,
-                      "This case is not covered yet, when we have more than one DoF per (sub)entity!");
+                      "This case is not covered yet, when we have a variable number of DoFs per (sub)entity!");
         // find the (sub)entity for this key
         const auto sub_entity = local_key.subEntity();
         const auto codim = local_key.codim();
@@ -149,8 +152,7 @@ public:
       DUNE_THROW(Exceptions::mapper_error,
                  "local_size(element) = " << coeffs.size() << "\n   local_index = " << local_index);
     const auto& local_key = coeffs.local_key(local_index);
-    // No need to assert local_key.index() == 0, has been checked in the ctor!
-    return mapper_->subIndex(element, local_key.subEntity(), local_key.codim());
+    return r * mapper_->subIndex(element, local_key.subEntity(), local_key.codim()) + local_key.index();
   } // ... mapToGlobal(...)
 
   using BaseType::global_indices;
@@ -164,7 +166,7 @@ public:
     for (size_t ii = 0; ii < local_sz; ++ii) {
       const auto& local_key = coeffs.local_key(ii);
       // No need to assert local_key.index() == 0, has been checked in the ctor!
-      indices[ii] = mapper_->subIndex(element, local_key.subEntity(), local_key.codim());
+      indices[ii] = r * mapper_->subIndex(element, local_key.subEntity(), local_key.codim()) + local_key.index();
     }
   } // ... globalIndices(...)