From 6233f08e287865edf7af142594179620243f39de Mon Sep 17 00:00:00 2001
From: Carlos Brito <carlos.brito524@gmail.com>
Date: Wed, 31 Jan 2024 17:56:52 +0100
Subject: [PATCH] Metaprogramming facilities

---
 examples/PDEs/Darcy.cpp               |  42 ++++++++-
 src/Rodin/Tuple.h                     |  51 +++++++---
 src/Rodin/Utility/Extract.h           |  26 +++++
 src/Rodin/Utility/ParameterPack.h     |   6 ++
 src/Rodin/Utility/Product.h           |  50 ++++++++++
 src/Rodin/Variational.h               |   1 +
 src/Rodin/Variational/Problem.h       | 131 ++++++++++++++++----------
 src/Rodin/Variational/ShapeFunction.h |  48 ++++------
 8 files changed, 266 insertions(+), 89 deletions(-)
 create mode 100644 src/Rodin/Utility/Extract.h
 create mode 100644 src/Rodin/Utility/Product.h

diff --git a/examples/PDEs/Darcy.cpp b/examples/PDEs/Darcy.cpp
index 897b4cb03..2b24cf7c5 100644
--- a/examples/PDEs/Darcy.cpp
+++ b/examples/PDEs/Darcy.cpp
@@ -13,6 +13,30 @@ using namespace Rodin;
 using namespace Rodin::Geometry;
 using namespace Rodin::Variational;
 
+    template <class T>
+    struct IsTrialFunctionReferenceWrapper
+    {
+      static constexpr Boolean Value = false;
+    };
+
+    template <class T>
+    struct IsTrialFunctionReferenceWrapper<std::reference_wrapper<T>>
+    {
+      static constexpr Boolean Value = IsTrialFunction<T>::Value;
+    };
+
+    template <class T>
+    struct IsTestFunctionReferenceWrapper
+    {
+      static constexpr Boolean Value = false;
+    };
+
+    template <class T>
+    struct IsTestFunctionReferenceWrapper<std::reference_wrapper<T>>
+    {
+      static constexpr Boolean Value = IsTestFunction<T>::Value;
+    };
+
 template <class T>
 struct IsFloatingPoint
 {
@@ -25,6 +49,13 @@ auto foo(const T& v)
   return 2 * v;
 }
 
+
+template <class L, class R>
+struct Pair
+{
+};
+
+
 int main(int argc, char** argv)
 {
   // Load mesh
@@ -38,7 +69,16 @@ int main(int argc, char** argv)
   TrialFunction u(vh);
   TestFunction  v(vh);
 
-  Problem problem(u, v, u, v, v);
+  Tuple us{std::ref(u), std::ref(v), std::ref(u), std::ref(v), std::ref(v), std::ref(u), std::ref(u)};
+
+  using Miaow = Utility::Product<Tuple<int, char, bool>, Tuple<float, double>>::Type<Pair>;
+  auto res = us.filter<IsTrialFunctionReferenceWrapper>();
+
+  Miaow t;
+  auto tt = t;
+  std::cout <<  res.size() << std::endl;
+  // Problem problem(u, v, u, v, v);
+
 
   // trw.apply([](auto&& v){ std::cout << v << " "; });
   std::cout << std::endl;
diff --git a/src/Rodin/Tuple.h b/src/Rodin/Tuple.h
index 0ef050a9f..d554dd24d 100644
--- a/src/Rodin/Tuple.h
+++ b/src/Rodin/Tuple.h
@@ -1,3 +1,9 @@
+/*
+ *          Copyright Carlos BRITO PACHECO 2021 - 2022.
+ * Distributed under the Boost Software License, Version 1.0.
+ *       (See accompanying file LICENSE or copy at
+ *          https://www.boost.org/LICENSE_1_0.txt)
+ */
 #ifndef RODIN_TUPLE_H
 #define RODIN_TUPLE_H
 
@@ -48,6 +54,13 @@ namespace Rodin
       constexpr
       void apply(Function&& func)
       {}
+
+      inline
+      constexpr
+      size_t size() const
+      {
+        return 0;
+      }
   };
 
   Tuple() -> Tuple<>;
@@ -79,7 +92,7 @@ namespace Rodin
       constexpr
       void apply(Function&& func)
       {
-        applyImpl(std::forward<Function>(func), std::index_sequence_for<Ts...>{});
+        applyImpl(std::forward<Function>(func), std::index_sequence_for<T, Ts...>{});
       }
 
       template <std::size_t Index>
@@ -111,16 +124,30 @@ namespace Rodin
       constexpr
       auto map(Func&& func) const
       {
-        return mapImpl(std::index_sequence_for<Ts...>(), std::forward<Func>(func));
+        return mapImpl(std::index_sequence_for<T, Ts...>(), std::forward<Func>(func));
+      }
+
+      inline
+      constexpr
+      Tuple concatenate(const Tuple<>& other) const
+      {
+        return *this;
       }
 
       template <typename ... Gs>
       inline
       constexpr
-      Tuple<Ts..., Gs...> concatenate(const Tuple<Gs...>& other) const
+      Tuple<T, Ts..., Gs...> concatenate(const Tuple<Gs...>& other) const
       {
         return concatenateImpl(other,
-            std::index_sequence_for<Ts...>(), std::index_sequence_for<Gs...>());
+            std::index_sequence_for<T, Ts...>(), std::index_sequence_for<Gs...>());
+      }
+
+      inline
+      constexpr
+      size_t size() const
+      {
+        return sizeof...(Ts) + 1;
       }
 
     private:
@@ -129,26 +156,28 @@ namespace Rodin
       constexpr
       void applyImpl(Function&& func, std::index_sequence<Indices...>)
       {
-        (func(std::get<Indices>(*this)), ...);
+        (func(get<Indices>()), ...);
       }
 
-      template <std::size_t ... Is, typename Func>
+      template <size_t I0, std::size_t ... Is, typename Func>
       inline
       constexpr
-      auto mapImpl(std::index_sequence<Is...>, Func&& func) const
+      auto mapImpl(std::index_sequence<I0, Is...>, Func&& func) const
       {
-        return Tuple<std::invoke_result_t<Func, Ts>...>(func(get<Is>())...);
+        return Tuple<
+          std::invoke_result_t<Func, T>,
+          std::invoke_result_t<Func, Ts>...>(func(get<I0>()), func(get<Is>())...);
       }
 
       template <typename ... Gs, std::size_t... Indices1, std::size_t... Indices2>
       inline
       constexpr
-      Tuple<Ts..., Gs...> concatenateImpl(
+      Tuple<T, Ts..., Gs...> concatenateImpl(
           const Tuple<Gs...>& other,
           std::index_sequence<Indices1...>,
           std::index_sequence<Indices2...>) const
       {
-        return Tuple<Ts..., Gs...>(get<Indices1>()..., other.template get<Indices2>()...);
+        return Tuple<T, Ts..., Gs...>{get<Indices1>()..., other.template get<Indices2>()...};
       }
 
       template <std::size_t Index, template <class> class Predicate>
@@ -156,7 +185,7 @@ namespace Rodin
       constexpr
       auto filterImpl() const
       {
-        if constexpr(Index == sizeof...(Ts))
+        if constexpr(Index == sizeof...(Ts) + 1)
         {
           return Tuple<>{};
         }
diff --git a/src/Rodin/Utility/Extract.h b/src/Rodin/Utility/Extract.h
new file mode 100644
index 000000000..f9c2d9939
--- /dev/null
+++ b/src/Rodin/Utility/Extract.h
@@ -0,0 +1,26 @@
+/*
+ *          Copyright Carlos BRITO PACHECO 2021 - 2022.
+ * Distributed under the Boost Software License, Version 1.0.
+ *       (See accompanying file LICENSE or copy at
+ *          https://www.boost.org/LICENSE_1_0.txt)
+ */
+#ifndef RODIN_UTILITY_EXTRACT_H
+#define RODIN_UTILITY_EXTRACT_H
+
+#include "Rodin/Tuple.h"
+
+namespace Rodin::Utility
+{
+  template <class ...>
+  struct Extract;
+
+  template <class T, class ... Ts>
+  struct Extract<Tuple<T, Ts...>>
+  {
+    template <template <class> class Extractor>
+    using Type = Tuple<typename Extractor<T>::Type, typename Extractor<Ts>::Type...>;
+  };
+}
+
+#endif
+
diff --git a/src/Rodin/Utility/ParameterPack.h b/src/Rodin/Utility/ParameterPack.h
index 72e83e78c..03022d6d8 100644
--- a/src/Rodin/Utility/ParameterPack.h
+++ b/src/Rodin/Utility/ParameterPack.h
@@ -1,3 +1,9 @@
+/*
+ *          Copyright Carlos BRITO PACHECO 2021 - 2022.
+ * Distributed under the Boost Software License, Version 1.0.
+ *       (See accompanying file LICENSE or copy at
+ *          https://www.boost.org/LICENSE_1_0.txt)
+ */
 #ifndef RODIN_UTILITY_PARAMETERPACK_H
 #define RODIN_UTILITY_PARAMETERPACK_H
 
diff --git a/src/Rodin/Utility/Product.h b/src/Rodin/Utility/Product.h
new file mode 100644
index 000000000..9d7b69185
--- /dev/null
+++ b/src/Rodin/Utility/Product.h
@@ -0,0 +1,50 @@
+/*
+ *          Copyright Carlos BRITO PACHECO 2021 - 2022.
+ * Distributed under the Boost Software License, Version 1.0.
+ *       (See accompanying file LICENSE or copy at
+ *          https://www.boost.org/LICENSE_1_0.txt)
+ */
+#ifndef RODIN_UTILITY_PRODUCT_H
+#define RODIN_UTILITY_PRODUCT_H
+
+#include "Rodin/Tuple.h"
+
+namespace Rodin::Utility
+{
+  template <class ...>
+  struct Product;
+
+  template <class H1, class H2, class ... Hs, class ... Gs>
+  struct Product<Tuple<H1, H2, Hs...>, Tuple<Gs...>>
+  {
+    template <template <class, class> class Pair>
+    using Type =
+      decltype(
+        std::declval<Tuple<Pair<H1, Gs>...>>().concatenate(
+          std::declval<Tuple<Pair<H2, Gs>...>>()).concatenate(
+            std::declval<typename Product<Tuple<Hs...>, Tuple<Gs...>>::template Type<Pair>>()));
+  };
+
+  template <class H, class ... Gs>
+  struct Product<Tuple<H>, Tuple<Gs...>>
+  {
+    template <template <class, class> class Pair>
+    using Type = Tuple<Pair<H, Gs>...>;
+  };
+
+  template <class ... Gs>
+  struct Product<Tuple<>, Tuple<Gs...>>
+  {
+    using Type = Tuple<>;
+  };
+
+  template <class ... Gs>
+  struct Product<Tuple<Gs...>, Tuple<>>
+  {
+    using Type = Tuple<>;
+  };
+}
+
+#endif
+
+
diff --git a/src/Rodin/Variational.h b/src/Rodin/Variational.h
index 6e996e245..854550d01 100644
--- a/src/Rodin/Variational.h
+++ b/src/Rodin/Variational.h
@@ -17,6 +17,7 @@
 #include "Variational/GridFunction.h"
 #include "Variational/FiniteElementSpace.h"
 
+#include "Variational/ShapeFunction.h"
 #include "Variational/TrialFunction.h"
 #include "Variational/TestFunction.h"
 
diff --git a/src/Rodin/Variational/Problem.h b/src/Rodin/Variational/Problem.h
index 18e0e7f9a..91caf712a 100644
--- a/src/Rodin/Variational/Problem.h
+++ b/src/Rodin/Variational/Problem.h
@@ -10,6 +10,7 @@
 #include <set>
 #include <variant>
 #include <functional>
+#include <boost/mp11.hpp>
 
 #include "Rodin/Alert.h"
 #include "Rodin/Geometry.h"
@@ -19,6 +20,8 @@
 #include "Rodin/Math/BlockSparseMatrix.h"
 #include "Rodin/FormLanguage/Base.h"
 #include "Rodin/Tuple.h"
+#include "Rodin/Utility/Extract.h"
+#include "Rodin/Utility/Product.h"
 
 #include "ForwardDecls.h"
 
@@ -90,6 +93,9 @@ namespace Rodin::Variational
       virtual ProblemBase* copy() const noexcept override = 0;
   };
 
+  template <>
+  class Problem<> {};
+
   /**
    * @ingroup ProblemSpecializations
    * @brief General class to assemble linear systems with `Math::SparseMatrix`
@@ -236,61 +242,86 @@ namespace Rodin::Variational
   Problem(TrialFunction<TrialFES>&, TestFunction<TestFES>&)
     -> Problem<TrialFES, TestFES, typename TrialFES::Context, Math::SparseMatrix, Math::Vector>;
 
-  template <class ... FES, class ... Ds, ShapeFunctionSpaceType ... Space>
+  template <class U1, class U2, class ... Us>
   class Problem<
-    Tuple<ShapeFunction<Ds, FES, Space>...>,
-    Context::Sequential, Math::SparseMatrix, Math::Vector>
+      Tuple<U1, U2, Us...>,
+      Context::Sequential, Math::SparseMatrix, Math::Vector>
     : public ProblemBase<Math::SparseMatrix, Math::Vector>
   {
-    template <class T>
-    struct IsTrialShapeFunctionReferenceWrapper
-    {
-      static constexpr Boolean Value = false;
-    };
-
-    template <class T>
-    struct IsTrialShapeFunctionReferenceWrapper<std::reference_wrapper<T>>
-    {
-      static constexpr Boolean Value = IsTrialShapeFunction<T>::Value;
-    };
-
-    template <class T>
-    struct IsTestShapeFunctionReferenceWrapper
-    {
-      static constexpr Boolean Value = false;
-    };
-
-    template <class T>
-    struct IsTestShapeFunctionReferenceWrapper<std::reference_wrapper<T>>
-    {
-      static constexpr Boolean Value = IsTestShapeFunction<T>::Value;
-    };
-
-    using TrialShapeFunctionTuple =
-      decltype(std::declval<
-        Tuple<
-          std::reference_wrapper<ShapeFunction<Ds, FES, Space>>...
-            >>().template filter<IsTrialShapeFunctionReferenceWrapper>());
-
-    using TestShapeFunctionTuple =
-      decltype(std::declval<
-        Tuple<
-          std::reference_wrapper<ShapeFunction<Ds, FES, Space>>...
-            >>().template filter<IsTestShapeFunctionReferenceWrapper>());
-
     public:
       using Context = Context::Sequential;
       using OperatorType = Math::SparseMatrix;
       using VectorType = Math::Vector;
       using Parent = ProblemBase<Math::SparseMatrix, Math::Vector>;
 
-      Problem(ShapeFunction<Ds, FES, Space>&... us)
+    private:
+      template <class T>
+      struct GetFES;
+
+      template <class T>
+      struct GetFES<std::reference_wrapper<T>>
+      {
+        using Type = typename FormLanguage::Traits<T>::FES;
+      };
+
+      template <class T>
+      struct IsTrialFunctionReferenceWrapper
+      {
+        static constexpr Boolean Value = false;
+      };
+
+      template <class T>
+      struct IsTrialFunctionReferenceWrapper<std::reference_wrapper<T>>
+      {
+        static constexpr Boolean Value = IsTrialFunction<T>::Value;
+      };
+
+      template <class T>
+      struct IsTestFunctionReferenceWrapper
+      {
+        static constexpr Boolean Value = false;
+      };
+
+      template <class T>
+      struct IsTestFunctionReferenceWrapper<std::reference_wrapper<T>>
+      {
+        static constexpr Boolean Value = IsTestFunction<T>::Value;
+      };
+
+      using TrialFunctionTuple =
+        decltype(std::declval<
+          Tuple<
+            std::reference_wrapper<U1>,
+            std::reference_wrapper<U2>,
+            std::reference_wrapper<Us>...>>()
+            .template filter<IsTrialFunctionReferenceWrapper>());
+
+      using TestFunctionTuple =
+        decltype(std::declval<
+          Tuple<
+            std::reference_wrapper<U1>,
+            std::reference_wrapper<U2>,
+            std::reference_wrapper<Us>...>>()
+            .template filter<IsTestFunctionReferenceWrapper>());
+
+      using TrialFESTuple = typename Utility::Extract<TrialFunctionTuple>::template Type<GetFES>;
+
+      using TestFESTuple = typename Utility::Extract<TestFunctionTuple>::template Type<GetFES>;
+
+      template <class TrialFES, class TestFES>
+      using BilinearFormType = BilinearForm<TrialFES, TestFES, OperatorType>;
+
+      using BilinearFormTuple =
+        typename Utility::Product<TrialFESTuple, TestFESTuple>::template Type<BilinearFormType>;
+
+    public:
+      Problem(U1& u1, U2& u2, Us&... us)
         : m_us(
-            Tuple{std::reference_wrapper<ShapeFunction<Ds, FES, Space>>(us)...}
-            .template filter<IsTrialShapeFunctionReferenceWrapper>()),
+            Tuple{std::ref(u1), std::ref(u2), std::ref(us)...}
+            .template filter<IsTrialFunctionReferenceWrapper>()),
           m_vs(
-            Tuple{std::reference_wrapper<ShapeFunction<Ds, FES, Space>>(us)...}
-            .template filter<IsTestShapeFunctionReferenceWrapper>())
+            Tuple{std::ref(u1), std::ref(u2), std::ref(us)...}
+            .template filter<IsTestFunctionReferenceWrapper>())
       {}
 
       Problem& assemble() override
@@ -333,8 +364,10 @@ namespace Rodin::Variational
       }
 
     private:
-      TrialShapeFunctionTuple m_us;
-      TestShapeFunctionTuple  m_vs;
+      TrialFunctionTuple m_us;
+      TestFunctionTuple  m_vs;
+
+      FormLanguage::List<BilinearFormBase<OperatorType>> m_bfs;
 
       bool            m_assembled;
       VectorType      m_mass;
@@ -342,10 +375,10 @@ namespace Rodin::Variational
       OperatorType    m_stiffness;
   };
 
-  template <class ... FES, class ... Ds, ShapeFunctionSpaceType ... Space>
-  Problem(ShapeFunction<Ds, FES, Space>&... us)
+  template <class U1, class U2, class ... Us>
+  Problem(U1& u1, U2& u2, Us&... us)
     -> Problem<
-        Tuple<ShapeFunction<Ds, FES, Space>...>,
+        Tuple<U1, U2, Us...>,
         Context::Sequential, Math::SparseMatrix, Math::Vector>;
 }
 
diff --git a/src/Rodin/Variational/ShapeFunction.h b/src/Rodin/Variational/ShapeFunction.h
index 0005d289c..5e49d7029 100644
--- a/src/Rodin/Variational/ShapeFunction.h
+++ b/src/Rodin/Variational/ShapeFunction.h
@@ -1,3 +1,9 @@
+/*
+ *          Copyright Carlos BRITO PACHECO 2021 - 2024.
+ * Distributed under the Boost Software License, Version 1.0.
+ *       (See accompanying file LICENSE or copy at
+ *          https://www.boost.org/LICENSE_1_0.txt)
+ */
 #ifndef RODIN_VARIATIONAL_SHAPEFUNCTION_H
 #define RODIN_VARIATIONAL_SHAPEFUNCTION_H
 
@@ -61,41 +67,27 @@ namespace Rodin::Variational
   };
 
   template <class T>
-  struct IsTrialShapeFunction;
-
-  template <class Derived, class FES, ShapeFunctionSpaceType Space>
-  struct IsTrialShapeFunction<ShapeFunction<Derived, FES, Space>>
+  struct IsTrialFunction
   {
-    static constexpr Boolean Value = (Space == ShapeFunctionSpaceType::Trial);
-
-    constexpr
-    IsTrialShapeFunction(const ShapeFunction<Derived, FES, Space>&) {}
+    static constexpr Boolean Value = false;
+  };
 
-    inline
-    constexpr
-    Boolean operator()()
-    {
-      return Value;
-    }
+  template <class FES>
+  struct IsTrialFunction<TrialFunction<FES>>
+  {
+    static constexpr Boolean Value = true;
   };
 
   template <class T>
-  struct IsTestShapeFunction;
-
-  template <class Derived, class FES, ShapeFunctionSpaceType Space>
-  struct IsTestShapeFunction<ShapeFunction<Derived, FES, Space>>
+  struct IsTestFunction
   {
-    static constexpr Boolean Value = (Space == ShapeFunctionSpaceType::Test);
-
-    constexpr
-    IsTestShapeFunction(const ShapeFunction<Derived, FES, Space>&) {}
+    static constexpr Boolean Value = false;
+  };
 
-    inline
-    constexpr
-    Boolean operator()()
-    {
-      return Value;
-    }
+  template <class FES>
+  struct IsTestFunction<TestFunction<FES>>
+  {
+    static constexpr Boolean Value = true;
   };
 
   template <