Revert "Support dynamic quadrature in Laplacian"

Q-Minh · Q-Minh · commit 0da6085d5204 · 2024-11-07T20:14:11.000-05:00
This reverts commit 649c3fb.
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<IndexVectorX const> const& eg,
-        Eigen::Ref<MatrixX const> const& wg,
-        Eigen::Ref<MatrixX const> const& GNeg,
-        int dims);
+        Eigen::Ref<MatrixX const> const& detJe,
+        Eigen::Ref<MatrixX const> const& GNe,
+        int dims,
+        int qOrder);
 
     Laplacian(Laplacian const&)            = delete;
     Laplacian& operator=(Laplacian const&) = delete;
@@ -40,6 +40,9 @@ class Laplacian
     int mMeshDims;
     int mMeshOrder;
     int mOrder;
+    int mQuadratureOrder;
+
+    static auto constexpr kMaxQuadratureOrder = 4;
 
   private:
     void* mLaplacian;
@@ -53,52 +56,58 @@ 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("eg"),
-            pyb::arg("wg"),
+            pyb::arg("detJe"),
             pyb::arg("GNe"),
-            pyb::arg("dims") = 1,
+            pyb::arg("dims")             = 1,
+            pyb::arg("quadrature_order") = 1,
             "Construct the symmetric part of the Laplacian operator on mesh mesh, using "
-            "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.")
+            "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.")
         .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(
-            "deltag",
+            "deltaE",
             [](Laplacian const& L) { return L.ElementLaplacians(); },
-            [](Laplacian& L, Eigen::Ref<MatrixX const> const& deltag) {
-                L.ElementLaplacians() = deltag;
+            [](Laplacian& L, Eigen::Ref<MatrixX const> const& deltaE) {
+                L.ElementLaplacians() = deltaE;
             },
-            "|#element nodes|x|#element nodes * #quad.pts.| matrix of element Laplacians")
+            "|#element nodes|x|#element nodes * #elements| matrix of element Laplacians")
         .def_property_readonly("shape", &Laplacian::Shape)
         .def("to_matrix", &Laplacian::ToMatrix);
 }
 
 Laplacian::Laplacian(
     Mesh const& M,
-    Eigen::Ref<IndexVectorX const> const& eg,
-    Eigen::Ref<MatrixX const> const& wg,
+    Eigen::Ref<MatrixX const> const& detJe,
     Eigen::Ref<MatrixX const> const& GNe,
-    int dims)
+    int dims,
+    int qOrder)
     : eMeshElement(M.eElement),
       mMeshDims(M.mDims),
       mMeshOrder(M.mOrder),
       mOrder(),
+      mQuadratureOrder(),
       mLaplacian(nullptr)
 {
-    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;
-    });
+    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);
 }
 
 CSCMatrix Laplacian::ToMatrix() const
@@ -120,20 +129,20 @@ std::tuple<Index, Index> Laplacian::Shape() const
 
 MatrixX const& Laplacian::ElementLaplacians() const
 {
-    MatrixX* deltagPtr;
+    MatrixX* deltaEPtr;
     Apply([&]<class LaplacianType>(LaplacianType* laplacian) {
-        deltagPtr = std::addressof(laplacian->deltag);
+        deltaEPtr = std::addressof(laplacian->deltaE);
     });
-    return *deltagPtr;
+    return *deltaEPtr;
 }
 
 MatrixX& Laplacian::ElementLaplacians()
 {
-    MatrixX* deltagPtr;
+    MatrixX* deltaEPtr;
     Apply([&]<class LaplacianType>(LaplacianType* laplacian) {
-        deltagPtr = std::addressof(laplacian->deltag);
+        deltaEPtr = std::addressof(laplacian->deltaE);
     });
-    return *deltagPtr;
+    return *deltaEPtr;
 }
 
 int const& Laplacian::dims() const
@@ -163,11 +172,16 @@ Laplacian::~Laplacian()
 template <class Func>
 void Laplacian::Apply(Func&& f) const
 {
-    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);
-    });
+    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);
+        });
 }
 
 } // 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, egL, wgL, GNegL = pbat.fem.laplacian(mesh)
+    L, detJeL, GNeL = 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, egL, wgL, GNegL = pbat.fem.laplacian(
+    L, detJeL, GNeL = 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.| array of elements corresponding to quadrature points. Defaults to None.
-        wg (np.ndarray, optional): |#quad.pts.| array of quadrature weights. Defaults to None.
+        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.
         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, np.ndarray): 
+        (scipy.sparse.csc_matrix, 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, np.ndarray): 
+        (scipy.sparse.csc_matrix | pbat.fem.Gradient, 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, np.ndarray): 
+        (pbat.fem.HyperElasticPotential, np.ndarray, np.ndarray): 
     """
     if eg is None or wg is None or GNeg is None:
         wg = _fem.inner_product_weights(
@@ -122,37 +122,32 @@ def laplacian(
         mesh,
         dims: int = 1,
         quadrature_order: int = 1,
-        eg: np.ndarray = None,
-        wg: np.ndarray = None,
-        GNeg: np.ndarray = None,
+        detJe: np.ndarray = None,
+        GNe: 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.
-        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.
+        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.
         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, np.ndarray): 
+        (pbat.fem.Laplacian | scipy.sparse.csc_matrix, np.ndarray, np.ndarray): 
     """
-    if eg is None or wg is None or GNeg is None:
-        wg = _fem.inner_product_weights(
+    if detJe is None:
+        detJe = _fem.jacobian_determinants(
             mesh, quadrature_order=quadrature_order)
-        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(
+    if GNe is None:
+        GNe = _fem.shape_function_gradients(
             mesh, quadrature_order=quadrature_order)
-    L = _fem.Laplacian(mesh, eg, wg, GNeg, dims=dims,
+    L = _fem.Laplacian(mesh, detJe, GNe, dims=dims,
                        quadrature_order=quadrature_order)
     L = L.to_matrix() if as_matrix else L
-    return L, eg, wg, GNeg
+    return L, detJe, GNe
 
 
 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>;
-        VectorX const wg      = fem::InnerProductWeights<kQuadratureOrder>(mesh).reshaped();
-        CSCMatrix L           = LaplacianType(mesh, eg, wg, GNe).ToMatrix();
+        using LaplacianType   = fem::SymmetricLaplacianMatrix<Mesh, kQuadratureOrder>;
+        MatrixX const detJe   = fem::DeterminantOfJacobian<kQuadratureOrder>(mesh);
+        CSCMatrix L           = LaplacianType(mesh, detJe, 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,10 +47,11 @@ TEST_CASE("[fem] HyperElasticPotential")
 
         MeshType const M(V, C);
         VectorX const x       = M.X.reshaped();
-        MatrixX const wg      = fem::InnerProductWeights<kQuadratureOrder>(M).reshaped();
+        MatrixX const WG      = fem::InnerProductWeights<kQuadratureOrder>(M);
+        VectorX const wg      = WG.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.size() / M.E.cols())
+                                    .replicate(1, WG.rows())
                                     .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,16 +44,11 @@ TEST_CASE("[fem] LaplacianMatrix")
                     return (kOrder - 1) + (kOrder - 1);
             }();
 
-            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);
+            using LaplacianMatrix = fem::SymmetricLaplacianMatrix<Mesh, kQuadratureOrder>;
+            MatrixX const detJe   = fem::DeterminantOfJacobian<kQuadratureOrder>(mesh);
+            MatrixX const GNe     = fem::ShapeFunctionGradients<kQuadratureOrder>(mesh);
             CHECK(math::CLinearOperator<LaplacianMatrix>);
-            LaplacianMatrix matrixFreeLaplacian(mesh, eg, wg, GNeg, outDims);
+            LaplacianMatrix matrixFreeLaplacian(mesh, detJe, GNe, 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
