Merge branch 'feature/projection' of https://gitlab.delpinux.fr/code/pugs into feature/projection

src/algebra/TinyMatrix.hpp

@@ -43,6 +43,18 @@ class [[nodiscard]] TinyMatrix
     }
   }
 
+  template <typename... Args>
+  PUGS_FORCEINLINE constexpr void
+  _unpackVariadicInput(const TinyVector<M, T>& t, Args&&... args) noexcept
+  {
+    for (size_t j = 0; j < M; ++j) {
+      m_values[_index(j, N - 1 - sizeof...(args))] = t[j];
+    }
+    if constexpr (sizeof...(args) > 0) {
+      this->_unpackVariadicInput(std::forward<Args>(args)...);
+    }
+  }
+
  public:
   PUGS_INLINE
   constexpr bool
@@ -322,6 +334,13 @@ class [[nodiscard]] TinyMatrix
     this->_unpackVariadicInput(t, std::forward<Args>(args)...);
   }
 
+  template <typename... Args>
+  PUGS_INLINE explicit constexpr TinyMatrix(const TinyVector<M, T>& t, Args&&... args) noexcept
+  {
+    static_assert(sizeof...(args) == N - 1, "wrong number of parameters");
+    this->_unpackVariadicInput(t, std::forward<Args>(args)...);
+  }
+
   // One does not use the '=default' constructor to avoid (unexpected)
   // performances issues
   PUGS_INLINE

src/scheme/FluxingAdvectionSolver.cpp

@@ -192,7 +192,65 @@ template <>
 FaceValue<double>
 FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
 {
-  throw NotImplementedError("ViensViensViens");
+  const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
+  FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity());
+  auto old_xr = m_old_mesh->xr();
+  auto new_xr = m_new_mesh->xr();
+  parallel_for(
+    m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+      const auto& face_to_node = face_to_node_matrix[face_id];
+      if (face_to_node.size() == 4) {
+        const Rd& x0 = old_xr[face_to_node[0]];
+        const Rd& x1 = old_xr[face_to_node[1]];
+        const Rd& x2 = old_xr[face_to_node[2]];
+        const Rd& x3 = old_xr[face_to_node[3]];
+
+        const Rd& x4 = new_xr[face_to_node[0]];
+        const Rd& x5 = new_xr[face_to_node[1]];
+        const Rd& x6 = new_xr[face_to_node[2]];
+        const Rd& x7 = new_xr[face_to_node[3]];
+
+        const Rd& a1 = x6 - x1;
+        const Rd& a2 = x6 - x3;
+        const Rd& a3 = x6 - x4;
+
+        const Rd& b1 = x7 - x0;
+        const Rd& b2 = x5 - x0;
+        const Rd& b3 = x2 - x0;
+
+        TinyMatrix<3> M1(a1 + b1, a2, a3);
+        TinyMatrix<3> M2(b1, a2 + b2, a3);
+        TinyMatrix<3> M3(a1, b2, a3 + b3);
+
+        algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12;
+      } else if (face_to_node.size() == 3) {
+        const Rd& x0 = old_xr[face_to_node[0]];
+        const Rd& x1 = old_xr[face_to_node[1]];
+        const Rd& x2 = old_xr[face_to_node[2]];
+
+        const Rd& x3 = new_xr[face_to_node[0]];
+        const Rd& x4 = new_xr[face_to_node[1]];
+        const Rd& x5 = new_xr[face_to_node[2]];
+
+        const Rd& a1 = x5 - x1;
+        const Rd& a2 = x5 - x2;
+        const Rd& a3 = x5 - x3;
+
+        const Rd& b1 = x5 - x0;
+        const Rd& b2 = x4 - x0;
+        const Rd& b3 = x2 - x0;
+
+        TinyMatrix<3> M1(a1 + b1, a2, a3);
+        TinyMatrix<3> M2(b1, a2 + b2, a3);
+        TinyMatrix<3> M3(a1, b2, a3 + b3);
+
+        algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12;
+      } else {
+        throw NotImplementedError("Not implemented for non quad faces");
+      }
+    });
+
+  return algebraic_fluxing_volume;
 }
 
 template <size_t Dimension>

tests/test_TinyMatrix.cpp

@@ -32,10 +32,14 @@ TEST_CASE("TinyMatrix", "[algebra]")
   REQUIRE(TinyMatrix<5, 4, int>::NumberOfColumns == 4);
 
   TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
+  TinyVector<3, int> x1(1, 5, 9);
+  TinyVector<3, int> x2(2, 6, 10);
+  TinyVector<3, int> x3(3, 7, 11);
+  TinyVector<3, int> x4(4, 8, 12);
   REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(0, 3) == 4) and   //
            (A(1, 0) == 5) and (A(1, 1) == 6) and (A(1, 2) == 7) and (A(1, 3) == 8) and   //
            (A(2, 0) == 9) and (A(2, 1) == 10) and (A(2, 2) == 11) and (A(2, 3) == 12)));
-
+  REQUIRE(A == TinyMatrix<3, 4, int>(x1, x2, x3, x4));
   TinyMatrix<3, 4, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7, 1, 3, 6);
 
   SECTION("checking for opposed matrix")