[spaces.mapper] add FixedOrderScalarDiscontinuousMapper

dune/gdt/spaces/mapper/default.hh

@@ -33,6 +33,9 @@ namespace GDT {
 template <class GL, class FiniteElementType>
 class FixedOrderMultipleCodimMultipleGeomTypeMapper;
 
+template <class GL, class FiniteElementType>
+class FixedOrderScalarDiscontinuousMapper;
+
 
 namespace internal {
 
@@ -71,6 +74,27 @@ private:
 }; // class FixedOrderMultipleCodimMultipleGeomTypeMapperTraits
 
 
+template <class GL, class FiniteElementType>
+class FixedOrderScalarDiscontinuousMapperTraits
+{
+  static_assert(XT::Grid::is_layer<GL>::value, "");
+
+  template <int dim_>
+  struct GeometryTypeLayout
+  {
+    bool contains(const GeometryType& gt) const
+    {
+      return gt.dim() == dim_;
+    }
+  };
+
+public:
+  using derived_type = FixedOrderScalarDiscontinuousMapper<GL, FiniteElementType>;
+  using BackendType = MultipleCodimMultipleGeomTypeMapper<GL, GeometryTypeLayout>;
+  using EntityType = XT::Grid::extract_entity_t<GL>;
+};
+
+
 } // namespace internal
 
 
@@ -218,6 +242,112 @@ private:
 }; // class FixedOrderMultipleCodimMultipleGeomTypeMapper
 
 
+template <class GL, class FiniteElementType>
+class FixedOrderScalarDiscontinuousMapper
+    : public MapperInterface<internal::FixedOrderScalarDiscontinuousMapperTraits<GL, FiniteElementType>>
+{
+public:
+  using Traits = internal::FixedOrderScalarDiscontinuousMapperTraits<GL, FiniteElementType>;
+
+private:
+  using ThisType = FixedOrderScalarDiscontinuousMapper<GL, FiniteElementType>;
+  using BaseType = MapperInterface<Traits>;
+  using D = typename GL::ctype;
+  static const constexpr size_t d = GL::dimension;
+
+public:
+  using typename BaseType::EntityType;
+  using typename BaseType::BackendType;
+
+  FixedOrderScalarDiscontinuousMapper(
+      const GL& grid_layer,
+      const std::shared_ptr<std::map<GeometryType, std::shared_ptr<FiniteElementType>>>& finite_elements)
+    : finite_elements_(finite_elements)
+    , mapper_(new BackendType(grid_layer))
+    , offset_(new std::vector<size_t>(mapper_->size()))
+    , max_num_dofs_(0)
+    , size_(0)
+  {
+    for (auto&& element : elements(grid_layer)) {
+      (*offset_)[mapper_->subIndex(element, 0, 0)] = size_;
+      const auto finite_element_search_result = finite_elements_->find(element.geometry().type());
+      if (finite_element_search_result == finite_elements_->end())
+        DUNE_THROW(mapper_error,
+                   "Missing finite element for the required geometry type " << element.geometry().type() << "!");
+      const auto& finite_element = *finite_element_search_result->second;
+      size_ += finite_element.size();
+      max_num_dofs_ = std::max(max_num_dofs_, XT::Common::numeric_cast<size_t>(finite_element.size()));
+    }
+  } // ... FixedOrderScalarDiscontinuousMapper(...)
+
+  FixedOrderScalarDiscontinuousMapper(const ThisType&) = default;
+  FixedOrderScalarDiscontinuousMapper(ThisType&&) = default;
+  FixedOrderScalarDiscontinuousMapper& operator=(const ThisType&) = delete;
+  FixedOrderScalarDiscontinuousMapper& operator=(ThisType&&) = delete;
+
+  const BackendType& backend() const
+  {
+    return *mapper_;
+  }
+
+  size_t size() const
+  {
+    return size_;
+  }
+
+  size_t maxNumDofs() const
+  {
+    return max_num_dofs_;
+  }
+
+  size_t numDofs(const EntityType& entity) const
+  {
+    const auto finite_element_search_result = finite_elements_->find(entity.geometry().type());
+    if (finite_element_search_result == finite_elements_->end())
+      DUNE_THROW(XT::Common::Exceptions::internal_error,
+                 "This must not happen after the checks in the ctor, the grid layer did not report all geometry types!"
+                     << "\n   entity.geometry().type() = "
+                     << entity.geometry().type());
+    const auto& finite_element = *finite_element_search_result->second;
+    return finite_element.size();
+  } // ... numDofs(...)
+
+  template <int cd, class GridImp, template <int, int, class> class EntityImp>
+  typename std::enable_if<cd != EntityType::codimension, size_t>::type
+  numDofs(const Entity<cd, EntityType::dimension, GridImp, EntityImp>& /*entity*/) const
+  {
+    return 0;
+  }
+
+  using BaseType::globalIndices;
+
+  void globalIndices(const EntityType& entity, DynamicVector<size_t>& ret) const
+  {
+    const size_t offset = (*offset_)[mapper_->subIndex(entity, 0, 0)];
+    const auto local_size = numDofs(entity);
+    if (ret.size() < local_size)
+      ret.resize(local_size, 0);
+    for (size_t ii = 0; ii < local_size; ++ii)
+      ret[ii] = offset + ii;
+  } // ... globalIndices(...)
+
+  size_t mapToGlobal(const EntityType& entity, const size_t local_index) const
+  {
+    if (local_index >= numDofs(entity))
+      DUNE_THROW(Exception, "numDofs(entity) = " << numDofs(entity) << "\n   local_index = " << local_index);
+    const size_t offset = (*offset_)[mapper_->subIndex(entity, 0, 0)];
+    return offset + local_index;
+  }
+
+private:
+  const std::shared_ptr<std::map<GeometryType, std::shared_ptr<FiniteElementType>>> finite_elements_;
+  const std::shared_ptr<BackendType> mapper_;
+  std::shared_ptr<std::vector<size_t>> offset_;
+  size_t max_num_dofs_;
+  size_t size_;
+}; // class FixedOrderScalarDiscontinuousMapper
+
+
 } // namespace GDT
 } // namespace Dune