Add utility to convert a preBasis to a flat preBasis

CMakeLists.txt

@@ -23,8 +23,6 @@ add_subdirectory("libs")
 add_subdirectory("test")
 
 dune_target_link_libraries(amdis fmt)
-target_compile_options(amdis PUBLIC "-Wall" "-Wpedantic")
-
 
 if (MTL_FOUND)
   dune_target_link_libraries(amdis MTL::MTL)
@@ -48,7 +46,8 @@ endif (PETSc_FOUND)
 
 option(ENABLE_ALL_WARNINGS "enable all meaningful warnings" OFF)
 if (ENABLE_ALL_WARNINGS)
-  target_compile_options(amdis PUBLIC "-Wextra" "-Wnon-virtual-dtor"
+  target_compile_options(amdis PUBLIC "-Wall" "-Wpedantic"
+                                      "-Wextra" "-Wnon-virtual-dtor"
                                       "-Wold-style-cast" "-Wcast-align"
                                       "-Woverloaded-virtual" "-Wconversion")
 endif (ENABLE_ALL_WARNINGS)

amdis/AdaptiveGrid.hpp

@@ -312,6 +312,11 @@ namespace AMDiS
     unsigned long changeIndex_ = 0;
   };
 
+  // deduction guide
+  template <class HostGrid>
+  AdaptiveGrid(HostGrid const&)
+    -> AdaptiveGrid<HostGrid>;
+
 
   template <class HostGrid>
   class AdaptiveGridFamily

amdis/DOFVector.hpp

@@ -222,7 +222,7 @@ namespace AMDiS
   // deduction guides
   template <class GB>
   DOFVector(GB&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
-    -> DOFVector<ParallelGlobalBasis<typename Underlying_t<GB>::PreBasis>>;
+    -> DOFVector<Underlying_t<GB>>;
 
   template <class GV, class PBF>
   DOFVector(GV const& gridView, PBF const& pbf, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
@@ -240,7 +240,7 @@ namespace AMDiS
    * \ref DataTransferOperation for more options.
    **/
   template <class ValueType = double, class GB>
-  DOFVector<ParallelGlobalBasis<typename Underlying_t<GB>::PreBasis>, ValueType>
+  DOFVector<Underlying_t<GB>, ValueType>
   makeDOFVector(GB&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
   {
     return {FWD(basis), op};

amdis/Observer.hpp

@@ -3,6 +3,7 @@
 #include <set>
 #include <utility>
 
+#include <amdis/Output.hpp>
 #include <amdis/common/ConceptsBase.hpp>
 #include <amdis/common/Index.hpp>
 #include <amdis/common/TypeTraits.hpp>
@@ -100,13 +101,22 @@ namespace AMDiS
       : public ObserverInterface<Event>
   {
   public:
-    template <class Notifier>
-    Observer(Notifier const& notifier)
-      : notifier_(const_cast<Notifier*>(&notifier))
+    template <class N,
+      REQUIRES(std::is_base_of_v<Notifier<Event>,N>)>
+    Observer(N const& notifier)
+      : notifier_(const_cast<N*>(&notifier))
     {
       notifier_->attach(this);
     }
 
+    template <class N,
+      REQUIRES(not std::is_base_of_v<Notifier<Event>,N>)>
+    Observer(N const& notifier)
+    {
+      warning("Ignoring Notifier. Use AdaptiveGrid wrapper.");
+    }
+
+
     /// Destructor, detaches from the notifier
     virtual ~Observer()
     {

amdis/functions/CMakeLists.txt

 install(FILES
+    FlatPreBasis.hpp
     FunctionFromCallable.hpp
     GlobalIdSet.hpp
     HierarchicNodeToRangeMap.hpp

amdis/functions/FlatPreBasis.hpp

+#pragma once
+
+#include <dune/common/version.hh>
+#include <dune/functions/functionspacebases/basistags.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/flatmultiindex.hh>
+
+#include <amdis/Output.hpp>
+#include <amdis/common/TypeTraits.hpp>
+#include <amdis/typetree/FiniteElementType.hpp>
+
+namespace AMDiS
+{
+  // Convert the index-merging strategy to FlatLexicographic or FlatInterleaved for
+  // composite or power bases, respectively.
+  template <class PreBasis,
+            class MultiIndex = Dune::Functions::FlatMultiIndex<std::size_t>>
+  struct FlatPreBasis
+  {
+    using type = PreBasis;
+
+    template <class PB>
+    static type create(PB const& preBasis)
+    {
+      return {preBasis.gridView()};
+    }
+
+    static PreBasis const& create(PreBasis const& preBasis)
+    {
+      return preBasis;
+    }
+  };
+
+  template <class PreBasis,
+            class MultiIndex = Dune::Functions::FlatMultiIndex<std::size_t>>
+  using FlatPreBasis_t = typename FlatPreBasis<PreBasis, MultiIndex>::type;
+
+  template <class PreBasis>
+  auto flatPreBasis(PreBasis const& preBasis)
+  {
+    return FlatPreBasis<PreBasis>::create(preBasis);
+  }
+
+  // specialization for Lagrange basis that needs an additional `order` parameter.
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+  template <class GV, class MI, class MultiIndex>
+  struct FlatPreBasis<Dune::Functions::LagrangePreBasis<GV,-1,MI>, MultiIndex>
+  {
+    using type = Dune::Functions::LagrangePreBasis<GV,-1,MultiIndex>;
+#else
+  template <class GV, class MI, class R, class MultiIndex>
+  struct FlatPreBasis<Dune::Functions::LagrangePreBasis<GV,-1,MI,R>, MultiIndex>
+  {
+    using type = Dune::Functions::LagrangePreBasis<GV,-1,MultiIndex,R>;
+#endif
+
+    template <class PB>
+    static type create(PB const& preBasis)
+    {
+      auto node = preBasis.makeNode();
+      node.bind(*preBasis.gridView().template begin<0>());
+      return {preBasis.gridView(), (unsigned int)(polynomialDegree(node))};
+    }
+
+    static type const& create(type const& preBasis)
+    {
+      return preBasis;
+    }
+  };
+
+  namespace Impl
+  {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+    // NOTE: dirty hack to get access to protected member variable, due to missing
+    // function in dune-functions-2.6
+    template <class PreBasis>
+    class DerivedPreBasis {
+    public:
+      // specialization for PowerPreBasis
+      template <class PB, class SPB = typename PB::SubPreBasis>
+      static SPB const& subPreBasis(PB const& pb)
+      {
+        return access(pb).subPreBasis_;
+      }
+
+      // specialization for CompositePreBasis
+      template <std::size_t I, class PB, class SPB = typename PB::template SubPreBasis<I>>
+      static SPB const& subPreBasis(PB const& pb, Dune::index_constant<I>)
+      {
+        return std::get<I>(access(pb).subPreBases_);
+      }
+
+      template <class T>
+      class Accessor : public T { friend class DerivedPreBasis<PreBasis>; };
+
+      template <class T>
+      static Accessor<T> const& access(T const& obj) { return static_cast<Accessor<T> const&>(obj); }
+    };
+
+    template <class PreBasis, class... Index>
+    static auto const& subPreBasis(PreBasis const& preBasis, Index... ii)
+    {
+      return DerivedPreBasis<PreBasis>::subPreBasis(preBasis,ii...);
+    }
+#else
+    template <class PreBasis, class... Index>
+    static auto const& subPreBasis(PreBasis const& preBasis, Index... ii)
+    {
+      return preBasis.subPreBasis(ii...);
+    }
+#endif
+  } // end namespace Impl
+
+
+  // specialization for composite bases
+  template <class MI, class IMS, class... SPB, class MultiIndex>
+  struct FlatPreBasis<Dune::Functions::CompositePreBasis<MI, IMS, SPB...>, MultiIndex>
+  {
+    using FIMS = Dune::Functions::BasisFactory::FlatLexicographic;
+    using type = Dune::Functions::CompositePreBasis<MultiIndex, FIMS, FlatPreBasis_t<SPB, MultiIndex>...>;
+
+    template <class PreBasis>
+    static type create(PreBasis const& preBasis)
+    {
+      return create(preBasis, std::index_sequence_for<SPB...>{});
+    }
+
+    template <class PreBasis, std::size_t... I>
+    static type create(PreBasis const& preBasis, std::index_sequence<I...>)
+    {
+      test_warning(std::is_same_v<IMS,FIMS>, "Basis converted into flat index-merging strategy.");
+      return {FlatPreBasis<SPB, MultiIndex>::create(Impl::subPreBasis(preBasis,Dune::index_constant<I>{}))...};
+    }
+  };
+
+  // specialization for flat power bases
+  template <class MI, class IMS, class SPB, std::size_t C, class MultiIndex>
+  struct FlatPreBasis<Dune::Functions::PowerPreBasis<MI, IMS, SPB, C>, MultiIndex>
+  {
+    using type = Dune::Functions::PowerPreBasis<MultiIndex, IMS, SPB, C>;
+
+    template <class PreBasis>
+    static type create(PreBasis const& preBasis)
+    {
+      return {FlatPreBasis<SPB, MultiIndex>::create(Impl::subPreBasis(preBasis))};
+    }
+  };
+
+  // specialization for blocked power bases
+  template <class MI, class SPB, std::size_t C, class MultiIndex>
+  struct FlatPreBasis<Dune::Functions::PowerPreBasis
+    <MI, Dune::Functions::BasisFactory::BlockedInterleaved, SPB, C>, MultiIndex>
+  {
+    using FIMS = Dune::Functions::BasisFactory::FlatInterleaved;
+    using type = Dune::Functions::PowerPreBasis<MultiIndex, FIMS, FlatPreBasis_t<SPB, MultiIndex>, C>;
+
+    template <class PreBasis>
+    static type create(PreBasis const& preBasis)
+    {
+      warning("Basis converted into flat index-merging strategy.");
+      return {FlatPreBasis<SPB, MultiIndex>::create(Impl::subPreBasis(preBasis))};
+    }
+  };
+
+  // specialization for blocked power bases
+  template <class MI, class SPB, std::size_t C, class MultiIndex>
+  struct FlatPreBasis<Dune::Functions::PowerPreBasis
+    <MI, Dune::Functions::BasisFactory::BlockedLexicographic, SPB, C>, MultiIndex>
+  {
+    using FIMS = Dune::Functions::BasisFactory::FlatLexicographic;
+    using type = Dune::Functions::PowerPreBasis<MultiIndex, FIMS, FlatPreBasis_t<SPB, MultiIndex>, C>;
+
+    template <class PreBasis>
+    static type create(PreBasis const& preBasis)
+    {
+      warning("Basis converted into flat index-merging strategy.");
+      return {FlatPreBasis<SPB, MultiIndex>::create(Impl::subPreBasis(preBasis))};
+    }
+  };
+
+} // end namespace AMDiS

amdis/functions/ParallelGlobalBasis.hpp

@@ -28,6 +28,7 @@
 #include <amdis/Output.hpp>
 #include <amdis/common/Concepts.hpp>
 #include <amdis/common/TypeTraits.hpp>
+#include <amdis/functions/FlatPreBasis.hpp>
 #include <amdis/linearalgebra/Traits.hpp>
 #include <amdis/typetree/MultiIndex.hpp>
 
@@ -103,7 +104,7 @@ namespace AMDiS
     template <class PBF>
     ParallelGlobalBasis(std::string const& name, GridView const& gridView, PBF const& preBasisFactory)
       : ParallelGlobalBasis(name, gridView.grid(),
-          preBasisFactory.template makePreBasis<MultiIndex_t<PBF>>(gridView))
+          flatPreBasis(preBasisFactory.template makePreBasis<MultiIndex_t<PBF>>(gridView)))
     {}
 
     /// Construct this global basis with empty name
@@ -179,13 +180,13 @@ namespace AMDiS
 
   // Deduction guides
   template <class GV, class PBF>
-  ParallelGlobalBasis(std::string const& name, GV const& gridView, PBF&& preBasisFactory)
-    -> ParallelGlobalBasis<decltype(
-         preBasisFactory.template makePreBasis<MultiIndex_t<PBF>>(gridView))>;
+  ParallelGlobalBasis(std::string const& name, GV const& gridView, PBF const& preBasisFactory)
+    -> ParallelGlobalBasis<decltype(flatPreBasis(
+          preBasisFactory.template makePreBasis<MultiIndex_t<PBF>>(gridView)))>;
 
   template <class GV, class PBF>
-  ParallelGlobalBasis(GV const& gridView, PBF&& preBasisFactory)
-    -> ParallelGlobalBasis<decltype(
-         preBasisFactory.template makePreBasis<MultiIndex_t<PBF>>(gridView))>;
+  ParallelGlobalBasis(GV const& gridView, PBF const& preBasisFactory)
+    -> ParallelGlobalBasis<decltype(flatPreBasis(
+          preBasisFactory.template makePreBasis<MultiIndex_t<PBF>>(gridView)))>;
 
 } // end namespace AMDiS

examples/ellipt.cc

@@ -8,6 +8,7 @@
 using Grid = Dune::YaspGrid<GRIDDIM>;
 
 using namespace AMDiS;
+using namespace Dune::Functions::BasisFactory;
 
 int main(int argc, char** argv)
 {
@@ -28,8 +29,9 @@ int main(int argc, char** argv)
     return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
   };
 
-  using Param = LagrangeBasis<Grid, 2>;
-  ProblemStat<Param> prob("ellipt");
+  auto grid = MeshCreator<Grid>{"elliptMesh"}.create();
+
+  ProblemStat prob("ellipt", *grid, lagrange<2>());
   prob.initialize(INIT_ALL);
 
   auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);

examples/vecellipt.cc

@@ -7,15 +7,15 @@
 #include <amdis/ProblemStat.hpp>
 
 using namespace AMDiS;
-
-using ElliptParam   = YaspGridBasis<GRIDDIM, 2,2>;
-using ElliptProblem = ProblemStat<ElliptParam>;
+using namespace Dune::Functions::BasisFactory;
 
 int main(int argc, char** argv)
 {
   Environment env(argc, argv);
 
-  ElliptProblem prob("ellipt");
+  using Grid = Dune::YaspGrid<GRIDDIM>;
+  auto grid = MeshCreator<Grid>{"elliptMesh"}.create();
+  ProblemStat prob("ellipt", *grid, power<2>(lagrange<2>()));
   prob.initialize(INIT_ALL);
 
   AdaptInfo adaptInfo("adapt");

test/CMakeLists.txt

@@ -60,6 +60,9 @@ dune_add_test(SOURCES FilesystemTest.cpp
 dune_add_test(SOURCES FlatMatVecTest.cpp
   LINK_LIBRARIES amdis)
 
+dune_add_test(SOURCES FlatPreBasisTest.cpp
+  LINK_LIBRARIES amdis)
+
 dune_add_test(SOURCES GlobalIdSetTest.cpp
   LINK_LIBRARIES amdis
   MPI_RANKS 2

test/FlatPreBasisTest.cpp

+#include <amdis/AMDiS.hpp>
+#include <amdis/DOFVector.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/functions/ParallelGlobalBasis.hpp>
+
+#include <dune/grid/yaspgrid.hh>
+
+int test()
+{
+  using namespace AMDiS;
+  using namespace Dune::Functions::BasisFactory;
+
+  // construct a grid
+  Dune::YaspGrid<2> grid({1.0, 1.0}, {8,8});
+  auto gridView = grid.leafGridView();
+
+  // construct a basis
+  ParallelGlobalBasis basis1(gridView, power<2>(lagrange<1>(), blockedInterleaved()));
+  ParallelGlobalBasis basis2(gridView, power<2>(lagrange<1>(), blockedLexicographic()));
+  ParallelGlobalBasis basis3(gridView, power<2>(lagrange<1>(), flatInterleaved()));
+  ParallelGlobalBasis basis4(gridView, power<2>(lagrange<1>(), flatLexicographic()));
+  ParallelGlobalBasis basis5(gridView, power<2>(lagrange<1>()));
+  ParallelGlobalBasis basis6(gridView, composite(lagrange<1>()));
+
+  // construct a ProblemStat
+  ProblemStat prob1("prob1", grid, power<2>(lagrange<1>(), blockedInterleaved()));
+  ProblemStat prob2("prob2", grid, power<2>(lagrange<1>(), blockedLexicographic()));
+  ProblemStat prob3("prob3", grid, power<2>(lagrange<1>(), flatInterleaved()));
+  ProblemStat prob4("prob4", grid, power<2>(lagrange<1>(), flatLexicographic()));
+  ProblemStat prob5("prob5", grid, power<2>(lagrange<1>()));
+  ProblemStat prob6("prob6", grid, composite(lagrange<1>()));
+
+  // construct a DOFVector
+  DOFVector vec1(gridView, power<2>(lagrange<1>(), blockedInterleaved()));
+  DOFVector vec2(gridView, power<2>(lagrange<1>(), blockedLexicographic()));
+  DOFVector vec3(gridView, power<2>(lagrange<1>(), flatInterleaved()));
+  DOFVector vec4(gridView, power<2>(lagrange<1>(), flatLexicographic()));
+  DOFVector vec5(gridView, power<2>(lagrange<1>()));
+  DOFVector vec6(gridView, composite(lagrange<1>()));
+
+  return 0;
+}
+
+int main(int argc, char** argv)
+{
+  AMDiS::Environment env(argc, argv);
+
+  return test();
+}