Reapply "Support dynamic quadrature in Laplacian"

Q-Minh · Q-Minh · commit 85b3c4cc2ee9 · 2024-11-07T20:15:11.000-05:00
This reverts commit 0da6085.
diff --git a/bindings/pypbat/fem/Laplacian.cpp b/bindings/pypbat/fem/Laplacian.cpp
@@ -15,10 +15,10 @@ class Laplacian
   public:
     Laplacian(
         Mesh const& M,
-        Eigen::Ref<MatrixX const> const& detJe,
-        Eigen::Ref<MatrixX const> const& GNe,
-        int dims,
-        int qOrder);
+        Eigen::Ref<IndexVectorX const> const& eg,
+        Eigen::Ref<MatrixX const> const& wg,
+        Eigen::Ref<MatrixX const> const& GNeg,
+        int dims);
 
     Laplacian(Laplacian const&)            = delete;
     Laplacian& operator=(Laplacian const&) = delete;
@@ -40,9 +40,6 @@ class Laplacian
     int mMeshDims;
     int mMeshOrder;
     int mOrder;
-    int mQuadratureOrder;
-
-    static auto constexpr kMaxQuadratureOrder = 4;
 
   private:
     void* mLaplacian;
@@ -56,58 +53,52 @@ void BindLaplacian(pybind11::module& m)
         .def(
             pyb::init<
                 Mesh const&,
+                Eigen::Ref<IndexVectorX const> const&,
                 Eigen::Ref<MatrixX const> const&,
                 Eigen::Ref<MatrixX const> const&,
-                int,
                 int>(),
             pyb::arg("mesh"),
-            pyb::arg("detJe"),
+            pyb::arg("eg"),
+            pyb::arg("wg"),
             pyb::arg("GNe"),
-            pyb::arg("dims")             = 1,
-            pyb::arg("quadrature_order") = 1,
+            pyb::arg("dims") = 1,
             "Construct the symmetric part of the Laplacian operator on mesh mesh, using "
-            "precomputed jacobian determinants detJe and shape function gradients GNe evaluated at "
-            "quadrature points given by the quadrature rule of order quadrature_order. The "
-            "discretization is based on Galerkin projection. The dimensions dims can be set to "
-            "accommodate vector-valued functions.")
+            "precomputed shape function gradients GNeg evaluated at quadrature points g at "
+            "elements eg with weights wg. The discretization is based on Galerkin projection. The "
+            "dimensions dims can be set to accommodate vector-valued functions.")
         .def_property(
             "dims",
             [](Laplacian const& L) { return L.dims(); },
             [](Laplacian& L, int dims) { L.dims() = dims; })
         .def_readonly("order", &Laplacian::mOrder)
-        .def_readonly("quadrature_order", &Laplacian::mQuadratureOrder)
         .def_property(
-            "deltaE",
+            "deltag",
             [](Laplacian const& L) { return L.ElementLaplacians(); },
-            [](Laplacian& L, Eigen::Ref<MatrixX const> const& deltaE) {
-                L.ElementLaplacians() = deltaE;
+            [](Laplacian& L, Eigen::Ref<MatrixX const> const& deltag) {
+                L.ElementLaplacians() = deltag;
             },
-            "|#element nodes|x|#element nodes * #elements| matrix of element Laplacians")
+            "|#element nodes|x|#element nodes * #quad.pts.| matrix of element Laplacians")
         .def_property_readonly("shape", &Laplacian::Shape)
         .def("to_matrix", &Laplacian::ToMatrix);
 }
 
 Laplacian::Laplacian(
     Mesh const& M,
-    Eigen::Ref<MatrixX const> const& detJe,
+    Eigen::Ref<IndexVectorX const> const& eg,
+    Eigen::Ref<MatrixX const> const& wg,
     Eigen::Ref<MatrixX const> const& GNe,
-    int dims,
-    int qOrder)
+    int dims)
     : eMeshElement(M.eElement),
       mMeshDims(M.mDims),
       mMeshOrder(M.mOrder),
       mOrder(),
-      mQuadratureOrder(),
       mLaplacian(nullptr)
 {
-    M.ApplyWithQuadrature<kMaxQuadratureOrder>(
-        [&]<pbat::fem::CMesh MeshType, auto QuadratureOrder>(MeshType* mesh) {
-            using LaplacianType = pbat::fem::SymmetricLaplacianMatrix<MeshType, QuadratureOrder>;
-            mLaplacian          = new LaplacianType(*mesh, detJe, GNe, dims);
-            mOrder              = LaplacianType::kOrder;
-            mQuadratureOrder    = LaplacianType::kQuadratureOrder;
-        },
-        qOrder);
+    M.Apply([&]<pbat::fem::CMesh MeshType>(MeshType* mesh) {
+        using LaplacianType = pbat::fem::SymmetricLaplacianMatrix<MeshType>;
+        mLaplacian          = new LaplacianType(*mesh, eg, wg, GNe, dims);
+        mOrder              = LaplacianType::kOrder;
+    });
 }
 
 CSCMatrix Laplacian::ToMatrix() const
@@ -129,20 +120,20 @@ std::tuple<Index, Index> Laplacian::Shape() const
 
 MatrixX const& Laplacian::ElementLaplacians() const
 {
-    MatrixX* deltaEPtr;
+    MatrixX* deltagPtr;
     Apply([&]<class LaplacianType>(LaplacianType* laplacian) {
-        deltaEPtr = std::addressof(laplacian->deltaE);
+        deltagPtr = std::addressof(laplacian->deltag);
     });
-    return *deltaEPtr;
+    return *deltagPtr;
 }
 
 MatrixX& Laplacian::ElementLaplacians()
 {
-    MatrixX* deltaEPtr;
+    MatrixX* deltagPtr;
     Apply([&]<class LaplacianType>(LaplacianType* laplacian) {
-        deltaEPtr = std::addressof(laplacian->deltaE);
+        deltagPtr = std::addressof(laplacian->deltag);
     });
-    return *deltaEPtr;
+    return *deltagPtr;
 }
 
 int const& Laplacian::dims() const
@@ -172,16 +163,11 @@ Laplacian::~Laplacian()
 template <class Func>
 void Laplacian::Apply(Func&& f) const
 {
-    ApplyToMeshWithQuadrature<kMaxQuadratureOrder>(
-        mMeshDims,
-        mMeshOrder,
-        eMeshElement,
-        mQuadratureOrder,
-        [&]<pbat::fem::CMesh MeshType, auto QuadratureOrder>() {
-            using LaplacianType = pbat::fem::SymmetricLaplacianMatrix<MeshType, QuadratureOrder>;
-            LaplacianType* laplacian = reinterpret_cast<LaplacianType*>(mLaplacian);
-            f.template operator()<LaplacianType>(laplacian);
-        });
+    ApplyToMesh(mMeshDims, mMeshOrder, eMeshElement, [&]<pbat::fem::CMesh MeshType>() {
+        using LaplacianType      = pbat::fem::SymmetricLaplacianMatrix<MeshType>;
+        LaplacianType* laplacian = reinterpret_cast<LaplacianType*>(mLaplacian);
+        f.template operator()<LaplacianType>(laplacian);
+    });
 }
 
 } // namespace fem
diff --git a/python/examples/diffusion_smoothing.py b/python/examples/diffusion_smoothing.py
@@ -22,7 +22,7 @@
     # Construct Galerkin laplacian, mass and gradient operators
     mesh = pbat.fem.Mesh(
         V.T, F.T, element=pbat.fem.Element.Triangle, order=1)
-    L, detJeL, GNeL = pbat.fem.laplacian(mesh)
+    L, egL, wgL, GNegL = pbat.fem.laplacian(mesh)
     M, detJeM = pbat.fem.mass_matrix(mesh, dims=1)
     # Setup heat diffusion
     dt = 0.016
diff --git a/python/examples/laplace.py b/python/examples/laplace.py
@@ -11,7 +11,7 @@ def harmonic_field(V: np.ndarray, C: np.ndarray, order: int, eps: float = 0.1):
     mesh = pbat.fem.Mesh(
         V.T, C.T, element=pbat.fem.Element.Tetrahedron, order=order)
     quadrature_order = max(int(2*(order-1)), 1)
-    L, detJeL, GNeL = pbat.fem.laplacian(
+    L, egL, wgL, GNegL = pbat.fem.laplacian(
         mesh, quadrature_order=quadrature_order)
     # Set Dirichlet boundary conditions at bottom and top of the model
     Xmin = mesh.X.min(axis=1)
diff --git a/python/pbatoolkit/fem.py b/python/pbatoolkit/fem.py
@@ -31,12 +31,12 @@ def divergence(mesh, quadrature_order: int = 1, eg: np.ndarray = None, wg: np.nd
     Args:
         mesh (pbat.fem.Mesh): 
         quadrature_order (int, optional): Polynomial order to use for operator construction. Defaults to 1.
-        eg (np.ndarray, optional): |#quad.pts.|x1 array of elements corresponding to quadrature points. Defaults to None.
-        wg (np.ndarray, optional): |#quad.pts.|x1 array of quadrature weights. Defaults to None.
+        eg (np.ndarray, optional): |#quad.pts.| array of elements corresponding to quadrature points. Defaults to None.
+        wg (np.ndarray, optional): |#quad.pts.| array of quadrature weights. Defaults to None.
         GNeg (np.ndarray, optional): |#elem. nodes|x|#dims * #quad.pts.| array of shape function gradients at quadrature points. Defaults to None.
 
     Returns:
-        (scipy.sparse.csc_matrix, np.ndarray): 
+        (scipy.sparse.csc_matrix, np.ndarray, np.ndarray): 
     """
     should_compute_quadrature = eg is None or wg is None or GNeg is None
     G, eg, GNeg = gradient(mesh, quadrature_order=quadrature_order, as_matrix=True) if should_compute_quadrature else gradient(
@@ -60,7 +60,7 @@ def gradient(mesh, quadrature_order: int = 1, eg: np.ndarray = None, GNeg: np.nd
         as_matrix (bool, optional): Return the operator as a sparse matrix. Defaults to True.
 
     Returns:
-        (scipy.sparse.csc_matrix | pbat.fem.Gradient, np.ndarray): 
+        (scipy.sparse.csc_matrix | pbat.fem.Gradient, np.ndarray, np.ndarray): 
     """
 
     if GNeg is None:
@@ -100,7 +100,7 @@ def hyper_elastic_potential(
         GNeg (np.ndarray, optional): Shape function gradients at quadrature points. Defaults to None.
 
     Returns:
-        (pbat.fem.HyperElasticPotential, np.ndarray, np.ndarray): 
+        (pbat.fem.HyperElasticPotential, np.ndarray, np.ndarray, np.ndarray): 
     """
     if eg is None or wg is None or GNeg is None:
         wg = _fem.inner_product_weights(
@@ -122,32 +122,37 @@ def laplacian(
         mesh,
         dims: int = 1,
         quadrature_order: int = 1,
-        detJe: np.ndarray = None,
-        GNe: np.ndarray = None,
+        eg: np.ndarray = None,
+        wg: np.ndarray = None,
+        GNeg: np.ndarray = None,
         as_matrix: bool = True):
     """Construct an FEM Laplacian operator
 
     Args:
         mesh (pbat.fem.Mesh): 
         dims (int, optional): Solution space dimensions. Defaults to 1.
         quadrature_order (int, optional): Polynomial order to use for operator construction. Defaults to 1.
-        detJe (np.ndarray, optional): Jacobian determinants at quadrature points. Defaults to None.
-        GNe (np.ndarray, optional): Shape function gradients at quadrature points. Defaults to None.
+        eg (np.ndarray, optional): |#quad.pts.| array of elements associated with quadrature points. Defaults to None.
+        wg (np.ndarray, optional): |#quad.pts.| array of quadrature weights. Defaults to None.
+        GNeg (np.ndarray, optional): Shape function gradients at quadrature points. Defaults to None.
         as_matrix (bool, optional): Return the operator as a sparse matrix. Defaults to True.
 
     Returns:
-        (pbat.fem.Laplacian | scipy.sparse.csc_matrix, np.ndarray, np.ndarray): 
+        (pbat.fem.Laplacian | scipy.sparse.csc_matrix, np.ndarray, np.ndarray, np.ndarray): 
     """
-    if detJe is None:
-        detJe = _fem.jacobian_determinants(
+    if eg is None or wg is None or GNeg is None:
+        wg = _fem.inner_product_weights(
             mesh, quadrature_order=quadrature_order)
-    if GNe is None:
-        GNe = _fem.shape_function_gradients(
+        eg = np.linspace(0, mesh.E.shape[1]-1,
+                         num=mesh.E.shape[1], dtype=np.int64)
+        eg = np.repeat(eg, wg.shape[0])
+        wg = wg.flatten(order="F")
+        GNeg = _fem.shape_function_gradients(
             mesh, quadrature_order=quadrature_order)
-    L = _fem.Laplacian(mesh, detJe, GNe, dims=dims,
+    L = _fem.Laplacian(mesh, eg, wg, GNeg, dims=dims,
                        quadrature_order=quadrature_order)
     L = L.to_matrix() if as_matrix else L
-    return L, detJe, GNe
+    return L, eg, wg, GNeg
 
 
 def mass_matrix(
diff --git a/source/pbat/fem/Gradient.cpp b/source/pbat/fem/Gradient.cpp
@@ -86,9 +86,9 @@ TEST_CASE("[fem] Gradient")
         CSCMatrix Lhat           = -GMT * Ihat.asDiagonal() * GM;
         VectorX const LhatValues = Lhat.coeffs();
         Lhat.coeffs().setOnes();
-        using LaplacianType   = fem::SymmetricLaplacianMatrix<Mesh, kQuadratureOrder>;
-        MatrixX const detJe   = fem::DeterminantOfJacobian<kQuadratureOrder>(mesh);
-        CSCMatrix L           = LaplacianType(mesh, detJe, GNe).ToMatrix();
+        using LaplacianType   = fem::SymmetricLaplacianMatrix<Mesh>;
+        VectorX const wg      = fem::InnerProductWeights<kQuadratureOrder>(mesh).reshaped();
+        CSCMatrix L           = LaplacianType(mesh, eg, wg, GNe).ToMatrix();
         VectorX const Lvalues = L.coeffs();
         L.coeffs().setOnes();
         Scalar const LsparsityError = (L - Lhat).squaredNorm();
diff --git a/source/pbat/fem/HyperElasticPotential.cpp b/source/pbat/fem/HyperElasticPotential.cpp
@@ -47,11 +47,10 @@ TEST_CASE("[fem] HyperElasticPotential")
 
         MeshType const M(V, C);
         VectorX const x       = M.X.reshaped();
-        MatrixX const WG      = fem::InnerProductWeights<kQuadratureOrder>(M);
-        VectorX const wg      = WG.reshaped();
+        MatrixX const wg      = fem::InnerProductWeights<kQuadratureOrder>(M).reshaped();
         MatrixX const GNeg    = fem::ShapeFunctionGradients<kQuadratureOrder>(M);
         IndexVectorX const eg = IndexVectorX::LinSpaced(M.E.cols(), Index(0), M.E.cols() - 1)
-                                    .replicate(1, WG.rows())
+                                    .replicate(1, wg.size() / M.E.cols())
                                     .transpose()
                                     .reshaped();
         ElasticPotentialType U(M, eg, wg, GNeg, x, Y, nu);
diff --git a/source/pbat/fem/LaplacianMatrix.cpp b/source/pbat/fem/LaplacianMatrix.cpp
@@ -44,11 +44,16 @@ TEST_CASE("[fem] LaplacianMatrix")
                     return (kOrder - 1) + (kOrder - 1);
             }();
 
-            using LaplacianMatrix = fem::SymmetricLaplacianMatrix<Mesh, kQuadratureOrder>;
-            MatrixX const detJe   = fem::DeterminantOfJacobian<kQuadratureOrder>(mesh);
-            MatrixX const GNe     = fem::ShapeFunctionGradients<kQuadratureOrder>(mesh);
+            using LaplacianMatrix = fem::SymmetricLaplacianMatrix<Mesh>;
+            VectorX const wg      = fem::InnerProductWeights<kQuadratureOrder>(mesh).reshaped();
+            IndexVectorX const eg =
+                IndexVectorX::LinSpaced(mesh.E.cols(), Index(0), mesh.E.cols() - 1)
+                    .replicate(1, wg.size() / mesh.E.cols())
+                    .transpose()
+                    .reshaped();
+            MatrixX const GNeg = fem::ShapeFunctionGradients<kQuadratureOrder>(mesh);
             CHECK(math::CLinearOperator<LaplacianMatrix>);
-            LaplacianMatrix matrixFreeLaplacian(mesh, detJe, GNe, outDims);
+            LaplacianMatrix matrixFreeLaplacian(mesh, eg, wg, GNeg, outDims);
 
             CSCMatrix const L = matrixFreeLaplacian.ToMatrix();
             CHECK_EQ(L.rows(), n);
diff --git a/source/pbat/fem/LaplacianMatrix.h b/source/pbat/fem/LaplacianMatrix.h
