[spaces.L2.DG] allow componentwise global mapping

dune/gdt/spaces/l2/discontinuous-lagrange.hh

@@ -75,10 +75,11 @@ private:
   using GlobalBasisImplementation = DefaultGlobalBasis<GridViewType, r, 1, R>;
 
 public:
-  DiscontinuousLagrangeSpace(GridViewType grd_vw, const int order = 1)
+  DiscontinuousLagrangeSpace(GridViewType grd_vw, const int order = 1, const bool dimws_glbl_mppng = false)
     : BaseType()
     , grid_view_(grd_vw)
     , order_(order)
+    , dimwise_global_mapping((r == 1) ? false : dimws_glbl_mppng) // does not make sense in the scalar case
     , local_finite_elements_(std::make_unique<const LocalLagrangeFiniteElementFamily<D, d, R, r>>())
     , mapper_(nullptr)
     , basis_(nullptr)
@@ -90,6 +91,7 @@ public:
     : BaseType(other)
     , grid_view_(other.grid_view_)
     , order_(other.order_)
+    , dimwise_global_mapping(other.dimwise_global_mapping)
     , local_finite_elements_(std::make_unique<const LocalLagrangeFiniteElementFamily<D, d, R, r>>())
     , mapper_(nullptr)
     , basis_(nullptr)
@@ -166,7 +168,8 @@ public:
     if (mapper_)
       mapper_->update_after_adapt();
     else
-      mapper_ = std::make_unique<MapperImplementation>(grid_view_, *local_finite_elements_, order_);
+      mapper_ =
+          std::make_unique<MapperImplementation>(grid_view_, *local_finite_elements_, order_, dimwise_global_mapping);
     // ... and basis
     if (basis_)
       basis_->update_after_adapt();
@@ -178,6 +181,11 @@ public:
 private:
   const GridViewType grid_view_;
   const int order_;
+
+public:
+  const bool dimwise_global_mapping;
+
+private:
   std::unique_ptr<const LocalLagrangeFiniteElementFamily<D, d, R, r>> local_finite_elements_;
   std::unique_ptr<MapperImplementation> mapper_;
   std::unique_ptr<GlobalBasisImplementation> basis_;
@@ -188,9 +196,10 @@ private:
  * \sa DiscontinuousLagrangeSpace
  */
 template <size_t r, class GV, class R = double>
-DiscontinuousLagrangeSpace<GV, r, R> make_discontinuous_lagrange_space(GV grid_view, const int order)
+DiscontinuousLagrangeSpace<GV, r, R>
+make_discontinuous_lagrange_space(GV grid_view, const int order, const bool dimwise_global_mapping = false)
 {
-  return DiscontinuousLagrangeSpace<GV, r, R>(grid_view, order);
+  return DiscontinuousLagrangeSpace<GV, r, R>(grid_view, order, dimwise_global_mapping);
 }
 
 
@@ -198,9 +207,10 @@ DiscontinuousLagrangeSpace<GV, r, R> make_discontinuous_lagrange_space(GV grid_v
  * \sa DiscontinuousLagrangeSpace
  */
 template <class GV, class R = double>
-DiscontinuousLagrangeSpace<GV, 1, R> make_discontinuous_lagrange_space(GV grid_view, const int order)
+DiscontinuousLagrangeSpace<GV, 1, R>
+make_discontinuous_lagrange_space(GV grid_view, const int order, const bool dimwise_global_mapping = false)
 {
-  return DiscontinuousLagrangeSpace<GV, 1, R>(grid_view, order);
+  return DiscontinuousLagrangeSpace<GV, 1, R>(grid_view, order, dimwise_global_mapping);
 }
 
 

dune/gdt/spaces/mapper/discontinuous.hh

@@ -39,6 +39,7 @@ public:
   using typename BaseType::D;
 
 private:
+  static constexpr size_t r = LocalFiniteElementFamily::r;
   using Implementation = MultipleCodimMultipleGeomTypeMapper<GV>;
 
 public:
@@ -47,13 +48,20 @@ public:
 
   DiscontinuousMapper(const GridViewType& grd_vw,
                       const LocalFiniteElementFamily& local_finite_elements,
-                      const int order)
+                      const int order,
+                      const bool dimwise_global_numbering = false)
     : grid_view_(grd_vw)
     , local_finite_elements_(local_finite_elements)
     , order_(order)
+    , dimwise_global_numbering_((r == 1) ? false : dimwise_global_numbering) // does not make sense in the scalar case
     , mapper_(grid_view_, mcmgElementLayout())
-    , offset_(mapper_.size())
+    , offset_()
   {
+    if (!dimwise_global_numbering_) // only use first component
+      offset_[0].resize(mapper_.size());
+    else
+      for (size_t s = 0; s < r; ++s)
+        offset_[s].resize(mapper_.size());
     this->update_after_adapt();
   }
 
@@ -76,7 +84,7 @@ public:
 
   size_t size() const override final
   {
-    return size_;
+    return std::accumulate(size_.begin(), size_.end(), 0);
   }
 
   size_t max_local_size() const override final
@@ -91,47 +99,98 @@ public:
 
   size_t global_index(const ElementType& element, const size_t local_index) const override final
   {
-    if (local_index >= local_size(element))
-      DUNE_THROW(Exceptions::mapper_error,
-                 "local_size(element) = " << local_size(element) << "\n   local_index = " << local_index);
-    return offset_[mapper_.index(element)] + local_index;
-  }
+    const auto local_sz = local_size(element);
+    DUNE_THROW_IF(local_index >= local_sz,
+                  Exceptions::mapper_error,
+                  "local_size(element) = " << local_sz << "\n   local_index = " << local_index);
+    if (!dimwise_global_numbering_)
+      return offset_[0][mapper_.index(element)] + local_index;
+    else {
+      assert(local_sz % r == 0 && "This should not happen, see update_after_adapt()!");
+      const size_t num_local_indices_per_range_dim = local_sz / r;
+      const size_t dim_range_associated_with_local_index = int(local_index / num_local_indices_per_range_dim);
+      const size_t local_index_within_dim_range = local_index % num_local_indices_per_range_dim;
+      size_t base_offset = 0;
+      for (size_t s = 0; s < dim_range_associated_with_local_index; ++s)
+        base_offset += size_[s];
+      return base_offset + offset_[dim_range_associated_with_local_index][mapper_.index(element)]
+             + local_index_within_dim_range;
+    }
+  } // ... global_index(...)
 
   using BaseType::global_indices;
 
   void global_indices(const ElementType& element, DynamicVector<size_t>& indices) const override final
   {
-    const size_t offset = offset_[mapper_.index(element)];
     const auto local_sz = local_size(element);
     if (indices.size() < local_sz)
       indices.resize(local_sz, 0);
-    for (size_t ii = 0; ii < local_sz; ++ii)
-      indices[ii] = offset + ii;
+    if (!dimwise_global_numbering_) {
+      const size_t offset = offset_[0][mapper_.index(element)];
+      for (size_t ii = 0; ii < local_sz; ++ii)
+        indices[ii] = offset + ii;
+    } else {
+      assert(local_sz % r == 0 && "This should not happen, see update_after_adapt()!");
+      const size_t num_local_indices_per_range_dim = local_sz / r;
+      for (size_t local_index = 0; local_index < local_sz; ++local_index) {
+        const size_t dim_range_associated_with_local_index = int(local_index / num_local_indices_per_range_dim);
+        const size_t local_index_within_dim_range = local_index % num_local_indices_per_range_dim;
+        size_t offset = 0;
+        for (size_t s = 0; s < dim_range_associated_with_local_index; ++s)
+          offset += size_[s];
+        offset += offset_[dim_range_associated_with_local_index][mapper_.index(element)];
+        indices[local_index] = offset + local_index_within_dim_range;
+      }
+    }
   } // ... global_indices(...)
 
   void update_after_adapt() override final
   {
     mapper_.update(grid_view_);
-    size_ = 0;
     max_num_dofs_ = 0;
-    if (offset_.size() != mapper_.size())
-      offset_.resize(mapper_.size());
+    std::fill(size_.begin(), size_.end(), 0);
+    auto prepare = [&](const auto& i) {
+      if (offset_[i].size() != mapper_.size())
+        offset_[i].resize(mapper_.size());
+    };
+    if (!dimwise_global_numbering_) // only use first component
+      prepare(0);
+    else
+      for (size_t s = 0; s < r; ++s)
+        prepare(s);
     for (auto&& element : elements(grid_view_)) {
-      offset_[mapper_.index(element)] = size_;
       const auto local_sz = this->local_size(element);
-      size_ += local_sz;
+      if (!dimwise_global_numbering_) {
+        offset_[0][mapper_.index(element)] = size_[0];
+        size_[0] += local_sz;
+      } else {
+        DUNE_THROW_IF(!local_finite_elements_.get(element.geometry().type(), order_).powered(),
+                      Exceptions::mapper_error,
+                      "Using dimwise numbering only makes sense for powered FEs!");
+        DUNE_THROW_IF(local_sz % r != 0,
+                      Exceptions::mapper_error,
+                      "Not a power FE, as number of local keys associated with an element not a multiple of "
+                      "dim_range:\n   this->local_size(element) = "
+                          << local_sz << "\n   "
+                          << "dim_range = " << size_t(r));
+        for (size_t s = 0; s < r; ++s) {
+          offset_[s][mapper_.index(element)] = size_[s];
+          size_[s] += local_sz / r;
+        }
+      }
       max_num_dofs_ = std::max(max_num_dofs_, local_sz);
     }
-  }
+  } // ... update_after_adapt(...)
 
 private:
   const GridViewType& grid_view_;
   const LocalFiniteElementFamily& local_finite_elements_;
   const int order_;
+  const bool dimwise_global_numbering_;
   Implementation mapper_;
-  std::vector<size_t> offset_;
+  std::array<std::vector<size_t>, r> offset_;
+  std::array<size_t, r> size_;
   size_t max_num_dofs_;
-  size_t size_;
 }; // class DiscontinuousMapper