[spaces] add RestrictedSpace

dune/gdt/exceptions.hh

@@ -29,6 +29,14 @@ class projection_error : public operator_error
 {
 };
 
+class space_error : public Dune::Exception
+{
+};
+
+class restricted_space_error : public space_error
+{
+};
+
 
 } // namespace GDT
 } // namespace Dune

dune/gdt/playground/spaces/mapper/restricted.hh

+// This file is part of the dune-gdt project:
+//   https://github.com/dune-community/dune-gdt
+// Copyright 2010-2017 dune-gdt developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+//      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+//          with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+//   Felix Schindler (2017)
+
+#ifndef DUNE_GDT_SPACES_MAPPER_RESTRICTED_HH
+#define DUNE_GDT_SPACES_MAPPER_RESTRICTED_HH
+
+#include <dune/grid/common/rangegenerators.hh>
+
+#include <dune/xt/common/exceptions.hh>
+#include <dune/xt/grid/type_traits.hh>
+
+#include <dune/gdt/spaces/mapper/interfaces.hh>
+#include <dune/gdt/spaces/interface.hh>
+
+namespace Dune {
+namespace GDT {
+
+
+// forward
+template <class UnrestrictedSpace, class RestrictionGridLayer>
+class RestrictedMapper;
+
+
+namespace internal {
+
+
+template <class UnrestrictedSpace, class RestrictionGridLayer>
+class RestrictedMapperTraits
+{
+  static_assert(XT::Grid::is_layer<RestrictionGridLayer>::value, "");
+  static_assert(is_space<UnrestrictedSpace>::value, "");
+
+public:
+  typedef RestrictedMapper<UnrestrictedSpace, RestrictionGridLayer> derived_type;
+  typedef typename UnrestrictedSpace::MapperType BackendType;
+  typedef XT::Grid::extract_entity_t<RestrictionGridLayer> EntityType;
+};
+
+
+} // namespace internal
+
+
+template <class UnrestrictedSpace, class RestrictionGridLayer>
+class RestrictedMapper
+    : public MapperInterface<internal::RestrictedMapperTraits<UnrestrictedSpace, RestrictionGridLayer>>
+{
+  typedef MapperInterface<internal::RestrictedMapperTraits<UnrestrictedSpace, RestrictionGridLayer>> BaseType;
+  typedef RestrictedMapper<UnrestrictedSpace, RestrictionGridLayer> ThisType;
+
+public:
+  using typename BaseType::BackendType;
+  using typename BaseType::EntityType;
+
+  RestrictedMapper(const UnrestrictedSpace& unrestricted_space, RestrictionGridLayer restriction_grid_layer)
+    : unrestricted_space_(unrestricted_space)
+    , grid_layer_(restriction_grid_layer)
+    , max_num_dofs_(0)
+    , map_to_restricted_()
+    , map_to_unrestricted_()
+  {
+    std::set<size_t> restricted_indices;
+    DynamicVector<size_t> unrestricted_indices(backend().maxNumDofs(), 0);
+    for (auto&& entity : elements(grid_layer_)) {
+      max_num_dofs_ = std::max(max_num_dofs_, backend().numDofs(entity));
+      backend().globalIndices(entity, unrestricted_indices);
+      for (size_t ii = 0; ii < backend().numDofs(entity); ++ii)
+        restricted_indices.insert(unrestricted_indices[ii]);
+    }
+    map_to_unrestricted_ = std::vector<size_t>(restricted_indices.size());
+    size_t counter = 0;
+    for (const auto& index : restricted_indices) {
+      map_to_restricted_[index] = counter;
+      map_to_unrestricted_[counter] = index;
+      ++counter;
+    }
+    if (map_to_restricted_.size() != map_to_unrestricted_.size())
+      DUNE_THROW(InvalidStateException, "This should not happen!");
+  }
+
+  RestrictedMapper(const ThisType& other) = default;
+  RestrictedMapper(ThisType&& source) = default;
+
+  ThisType& operator=(const ThisType& other) = delete;
+  ThisType& operator=(ThisType&& source) = delete;
+
+  template <class U, class V>
+  void restrict(const XT::LA::VectorInterface<U>& unrestricted_vector,
+                XT::LA::VectorInterface<V>& restricted_vector) const
+  {
+    if (unrestricted_vector.size() != backend().size())
+      DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match,
+                 "unrestricted_vector.size() = " << unrestricted_vector.size()
+                                                 << "\n   unrestricted_space.mapper().size() = "
+                                                 << backend().size());
+    if (restricted_vector.size() != size())
+      DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match,
+                 "restricted_vector.size() = " << restricted_vector.size() << "\n   size() = " << size());
+    // the actual work
+    for (size_t restricted_DoF = 0; restricted_DoF < map_to_unrestricted_.size(); ++restricted_DoF)
+      restricted_vector[restricted_DoF] = unrestricted_vector[map_to_unrestricted_[restricted_DoF]];
+  } // ... restrict(...)
+
+  template <class V>
+  typename V::derived_type restrict(const XT::LA::VectorInterface<V>& unrestricted_vector) const
+  {
+    typename V::derived_type restricted_vector(size(), 0.);
+    restrict(unrestricted_vector, restricted_vector);
+    return restricted_vector;
+  }
+
+  template <class U, class V>
+  void extend(const XT::LA::VectorInterface<U>& restricted_vector,
+              XT::LA::VectorInterface<V>& unrestricted_vector) const
+  {
+    if (restricted_vector.size() != size())
+      DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match,
+                 "restricted_vector.size() = " << restricted_vector.size() << "\n   size() = " << size());
+    if (unrestricted_vector.size() != backend().size())
+      DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match,
+                 "unrestricted_vector.size() = " << unrestricted_vector.size()
+                                                 << "\n   unrestricted_space.mapper().size() = "
+                                                 << backend().size());
+    // the actual work
+    unrestricted_vector *= 0.;
+    for (size_t restricted_DoF = 0; restricted_DoF < map_to_unrestricted_.size(); ++restricted_DoF)
+      unrestricted_vector[map_to_unrestricted_[restricted_DoF]] = restricted_vector[restricted_DoF];
+  } // ... extend(...)
+
+  template <class V>
+  typename V::derived_type extend(const XT::LA::VectorInterface<V>& restricted_vector) const
+  {
+    typename V::derived_type unrestricted_vector(backend().size(), 0.);
+    extend(restricted_vector, unrestricted_vector);
+    return unrestricted_vector;
+  }
+
+  const BackendType& backend() const
+  {
+    return unrestricted_space_.mapper();
+  }
+
+  size_t size() const
+  {
+    return map_to_restricted_.size();
+  }
+
+  size_t maxNumDofs() const
+  {
+    return max_num_dofs_;
+  }
+
+  size_t numDofs(const EntityType& entity) const
+  {
+    const auto error_message = check_entity(entity);
+    if (error_message.size() > 0)
+      DUNE_THROW(restricted_space_error, error_message);
+    return backend().numDofs(entity);
+  }
+
+  void globalIndices(const EntityType& entity, DynamicVector<size_t>& ret) const
+  {
+    const auto error_message = check_entity(entity);
+    if (error_message.size() > 0)
+      DUNE_THROW(restricted_space_error, error_message);
+    backend().globalIndices(entity, ret);
+    for (size_t ii = 0; ii < numDofs(entity); ++ii)
+      ret[ii] = map_to_restricted_.at(ret[ii]);
+  }
+
+  size_t mapToGlobal(const EntityType& entity, const size_t& localIndex) const
+  {
+    const auto error_message = check_entity(entity);
+    if (error_message.size() > 0)
+      DUNE_THROW(restricted_space_error, error_message);
+    return map_to_restricted_.at(backend().mapToGlobal(entity, localIndex));
+  }
+
+private:
+  std::string check_entity(const EntityType& entity) const
+  {
+    std::stringstream error_message;
+    if (!grid_layer_.indexSet().contains(entity)) {
+      if (unrestricted_space_.grid_layer().indexSet().contains(entity))
+        error_message << "Entity not contained in restriction grid layer, but contained in the unrestricted grid layer "
+                         "with index "
+                      << unrestricted_space_.grid_layer().indexSet().index(entity) << "!";
+      else
+        error_message << "Entity neither contained in restriction grid layer nor in the unrestricted grid layer!";
+    }
+    return error_message.str();
+  } // ... check_entity(...)
+
+  const UnrestrictedSpace unrestricted_space_;
+  const RestrictionGridLayer grid_layer_;
+  size_t max_num_dofs_;
+  std::map<size_t, size_t> map_to_restricted_;
+  std::vector<size_t> map_to_unrestricted_;
+}; // class RestrictedMapper
+
+
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_SPACES_MAPPER_RESTRICTED_HH

dune/gdt/playground/spaces/restricted.hh

+// This file is part of the dune-gdt project:
+//   https://github.com/dune-community/dune-gdt
+// Copyright 2010-2017 dune-gdt developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+//      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+//          with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+//   Felix Schindler (2017)
+
+#ifndef DUNE_GDT_SPACES_RESTRICTED_HH
+#define DUNE_GDT_SPACES_RESTRICTED_HH
+
+#include <sstream>
+
+#include <dune/xt/la/type_traits.hh>
+#include <dune/xt/grid/layers.hh>
+
+#include <dune/gdt/exceptions.hh>
+#include <dune/gdt/spaces/interface.hh>
+
+#include "mapper/restricted.hh"
+
+namespace Dune {
+namespace GDT {
+
+
+// forward
+template <class UnrestrictedSpace, class RestrictionGridLayer>
+class RestrictedSpace;
+
+
+namespace internal {
+
+
+template <class UnrestrictedSpace, class RestrictionGridLayer>
+class RestrictedSpaceTraits
+{
+  static_assert(XT::Grid::is_layer<RestrictionGridLayer>::value, "");
+  static_assert(is_space<UnrestrictedSpace>::value, "");
+
+  template <class G = RestrictionGridLayer, bool v = XT::Grid::is_view<G>::value, bool p = XT::Grid::is_part<G>::value>
+  struct layer_backend_helper
+  {
+    static_assert(AlwaysFalse<G>::value, "");
+  };
+
+  template <class G>
+  struct layer_backend_helper<G, true, false>
+  {
+    static const XT::Grid::Backends value = XT::Grid::Backends::view;
+  };
+
+  template <class G>
+  struct layer_backend_helper<G, false, true>
+  {
+    static const XT::Grid::Backends value = XT::Grid::Backends::part;
+  };
+
+public:
+  typedef RestrictedSpace<UnrestrictedSpace, RestrictionGridLayer> derived_type;
+  static const int polOrder = UnrestrictedSpace::polOrder;
+  static const bool continuous = UnrestrictedSpace::continuous;
+  typedef UnrestrictedSpace BackendType;
+  typedef RestrictedMapper<UnrestrictedSpace, RestrictionGridLayer> MapperType;
+  typedef typename UnrestrictedSpace::BaseFunctionSetType BaseFunctionSetType;
+  typedef typename UnrestrictedSpace::CommunicatorType CommunicatorType;
+  typedef RestrictionGridLayer GridLayerType;
+  typedef typename UnrestrictedSpace::RangeFieldType RangeFieldType;
+  static const XT::Grid::Backends layer_backend = layer_backend_helper<>::value;
+}; // class RestrictedSpaceTraits
+
+
+} // namespace internal
+
+
+template <class UnrestrictedSpace, class RestrictionGridLayer>
+class RestrictedSpace : public SpaceInterface<internal::RestrictedSpaceTraits<UnrestrictedSpace, RestrictionGridLayer>,
+                                              UnrestrictedSpace::dimDomain,
+                                              UnrestrictedSpace::dimRange,
+                                              UnrestrictedSpace::dimRangeCols>
+{
+  typedef SpaceInterface<internal::RestrictedSpaceTraits<UnrestrictedSpace, RestrictionGridLayer>,
+                         UnrestrictedSpace::dimDomain,
+                         UnrestrictedSpace::dimRange,
+                         UnrestrictedSpace::dimRangeCols>
+      BaseType;
+  typedef RestrictedSpace<UnrestrictedSpace, RestrictionGridLayer> ThisType;
+
+public:
+  typedef internal::RestrictedSpaceTraits<UnrestrictedSpace, RestrictionGridLayer> Traits;
+  using typename BaseType::MapperType;
+  using typename BaseType::GridLayerType;
+  using typename BaseType::BackendType;
+  using typename BaseType::EntityType;
+  using typename BaseType::BaseFunctionSetType;
+  using typename BaseType::CommunicatorType;
+  using typename BaseType::PatternType;
+
+  RestrictedSpace(const UnrestrictedSpace& unrestricted_space, RestrictionGridLayer restriction_grid_layer)
+    : unrestricted_space_(unrestricted_space)
+    , grid_layer_(restriction_grid_layer)
+    , mapper_(unrestricted_space_, grid_layer_)
+  {
+  }
+
+  RestrictedSpace(const ThisType& other) = default;
+  RestrictedSpace(ThisType&& source) = default;
+
+  ThisType& operator=(const ThisType& other) = delete;
+  ThisType& operator=(ThisType&& source) = delete;
+
+  const GridLayerType& grid_layer() const
+  {
+    return grid_layer_;
+  }
+
+  const BackendType& backend() const
+  {
+    return unrestricted_space_;
+  }
+
+  const MapperType& mapper() const
+  {
+    return mapper_;
+  }
+
+  BaseFunctionSetType base_function_set(const EntityType& entity) const
+  {
+    const auto error_message = check_entity(entity);
+    if (error_message.size() > 0)
+      DUNE_THROW(restricted_space_error, error_message);
+    return unrestricted_space_.base_function_set(entity);
+  }
+
+  CommunicatorType& communicator() const
+  {
+    return unrestricted_space_.communicator();
+  }
+
+  template <class S, size_t d, size_t r, size_t rC, class C>
+  void local_constraints(const SpaceInterface<S, d, r, rC>& ansatz_space,
+                         const EntityType& entity,
+                         ConstraintsInterface<C>& ret) const
+  {
+    const auto error_message = check_entity(entity);
+    if (error_message.size() > 0)
+      DUNE_THROW(restricted_space_error, error_message);
+    return unrestricted_space_.local_constraints(ansatz_space, entity, ret);
+  }
+
+  template <class GL, class S, size_t d, size_t r, size_t rC>
+  typename std::enable_if<XT::Grid::is_layer<GL>::value, PatternType>::type
+  compute_pattern(const GL& /*grd_layr*/, const SpaceInterface<S, d, r, rC>& /*ansatz_space*/) const
+  {
+    DUNE_THROW(NotImplemented, "Yet");
+  }
+
+private:
+  std::string check_entity(const EntityType& entity) const
+  {
+    std::stringstream error_message;
+    if (!grid_layer_.indexSet().contains(entity)) {
+      if (unrestricted_space_.grid_layer().indexSet().contains(entity))
+        error_message << "Entity not contained in restriction grid layer, but contained in the unrestricted grid layer "
+                         "with index "
+                      << unrestricted_space_.grid_layer().indexSet().index(entity) << "!";
+      else
+        error_message << "Entity neither contained in restriction grid layer nor in the unrestricted grid layer!";
+    }
+    return error_message.str();
+  } // ... check_entity(...)
+
+  const UnrestrictedSpace unrestricted_space_;
+  const GridLayerType grid_layer_;
+  const MapperType mapper_;
+}; // class RestrictedSpace
+
+
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_SPACES_RESTRICTED_HH