Support dynamic quadrature in Gradient operator

Author: Q-Minh

bindings/pypbat/fem/Gradient.cpp

@@ -13,7 +13,10 @@ namespace fem {
 class Gradient
 {
   public:
-    Gradient(Mesh const& M, Eigen::Ref<MatrixX const> const& GNe, int qOrder);
+    Gradient(
+        Mesh const& M,
+        Eigen::Ref<IndexVectorX const> const& eg,
+        Eigen::Ref<MatrixX const> const& GNeg);
 
     Gradient(Gradient const&)            = delete;
     Gradient& operator=(Gradient const&) = delete;
@@ -30,9 +33,6 @@ class Gradient
     int mMeshOrder;
     int mDims;
     int mOrder;
-    int mQuadratureOrder;
-
-    static auto constexpr kMaxQuadratureOrder = 2;
 
   private:
     void* mGradient;
@@ -44,37 +44,38 @@ void BindGradient(pybind11::module& m)
 
     pyb::class_<Gradient>(m, "Gradient")
         .def(
-            pyb::init<Mesh const&, Eigen::Ref<MatrixX const> const&, int>(),
+            pyb::init<
+                Mesh const&,
+                Eigen::Ref<IndexVectorX const> const&,
+                Eigen::Ref<MatrixX const> const&>(),
             pyb::arg("mesh"),
-            pyb::arg("GNe"),
-            pyb::arg("quadrature_order") = 1,
+            pyb::arg("eg"),
+            pyb::arg("GNeg"),
             "Construct Gradient operator from mesh mesh, using precomputed shape function "
-            "gradients GNe at quadrature points given by quadrature rule of order quadrature_order")
+            "gradients GNeg at quadrature points at elements eg.")
         .def_readonly("dims", &Gradient::mDims)
         .def_readonly("order", &Gradient::mOrder, "Polynomial order of the gradient")
-        .def_readonly("quadrature_order", &Gradient::mQuadratureOrder)
         .def_property_readonly("shape", &Gradient::Shape)
         .def("to_matrix", &Gradient::ToMatrix);
 }
 
-Gradient::Gradient(Mesh const& M, Eigen::Ref<MatrixX const> const& GNe, int qOrder)
+Gradient::Gradient(
+    Mesh const& M,
+    Eigen::Ref<IndexVectorX const> const& eg,
+    Eigen::Ref<MatrixX const> const& GNeg)
     : eMeshElement(M.eElement),
       mMeshDims(M.mDims),
       mMeshOrder(M.mOrder),
       mDims(),
       mOrder(),
-      mQuadratureOrder(),
       mGradient(nullptr)
 {
-    M.ApplyWithQuadrature<kMaxQuadratureOrder>(
-        [&]<pbat::fem::CMesh MeshType, auto QuadratureOrder>(MeshType* mesh) {
-            using GradientType = pbat::fem::Gradient<MeshType, QuadratureOrder>;
-            mGradient          = new GradientType(*mesh, GNe);
-            mDims              = GradientType::kDims;
-            mOrder             = GradientType::kOrder;
-            mQuadratureOrder   = GradientType::kQuadratureOrder;
-        },
-        qOrder);
+    M.Apply([&]<pbat::fem::CMesh MeshType>(MeshType* mesh) {
+        using GradientType = pbat::fem::Gradient<MeshType>;
+        mGradient          = new GradientType(*mesh, eg, GNeg);
+        mDims              = GradientType::kDims;
+        mOrder             = GradientType::kOrder;
+    });
 }
 
 CSCMatrix Gradient::ToMatrix() const
@@ -103,16 +104,11 @@ Gradient::~Gradient()
 template <class Func>
 void Gradient::Apply(Func&& f) const
 {
-    ApplyToMeshWithQuadrature<kMaxQuadratureOrder>(
-        mMeshDims,
-        mMeshOrder,
-        eMeshElement,
-        mQuadratureOrder,
-        [&]<pbat::fem::CMesh MeshType, auto QuadratureOrder>() {
-            using GradientType     = pbat::fem::Gradient<MeshType, QuadratureOrder>;
-            GradientType* gradient = reinterpret_cast<GradientType*>(mGradient);
-            f.template operator()<GradientType>(gradient);
-        });
+    ApplyToMesh(mMeshDims, mMeshOrder, eMeshElement, [&]<pbat::fem::CMesh MeshType>() {
+        using GradientType     = pbat::fem::Gradient<MeshType>;
+        GradientType* gradient = reinterpret_cast<GradientType*>(mGradient);
+        f.template operator()<GradientType>(gradient);
+    });
 }
 
 } // namespace fem

python/examples/elasticity.py

@@ -37,11 +37,6 @@
     M, detJeM = pbat.fem.mass_matrix(mesh, rho=rho)
     Minv = pbat.math.linalg.ldlt(M)
     Minv.compute(M)
-    # Could also lump the mass matrix like this
-    # lumpedm = M.sum(axis=0)
-    # M = sp.sparse.spdiags(lumpedm, np.array([0]), m=M.shape[0], n=M.shape[0])
-    # Minv = sp.sparse.spdiags(
-    #     1./lumpedm, np.array([0]), m=M.shape[0], n=M.shape[0])
 
     # Construct load vector from gravity field
     g = np.zeros(mesh.dims)

python/examples/heat.py

@@ -35,8 +35,9 @@
     # Construct Galerkin laplacian, mass and gradient operators
     n = V.shape[0]
     M, detJeM = pbat.fem.mass_matrix(mesh, dims=1)
-    G, GNeG = pbat.fem.gradient(mesh)
-    D, GNeD = pbat.fem.divergence(mesh, GNe=GNeG)
+    G, egG, GNegG = pbat.fem.gradient(mesh)
+    wgD = pbat.fem.inner_product_weights(mesh)
+    D, wgD, egD, GNegD = pbat.fem.divergence(mesh, eg=egG, wg=wgD, GNeg=GNegG)
     L = D @ G
     # Setup 1-step heat diffusion
     h = igl.avg_edge_length(V, C)

python/examples/ipc.py

@@ -367,8 +367,7 @@ def __call__(self, xk: np.ndarray):
     # Lumped mass matrix
     rho = args.rho
     M, detJeM = pbat.fem.mass_matrix(mesh, rho=rho, lump=True)
-    Minv = sp.sparse.spdiags(
-        1./M.diagonal(), np.array([0]), m=M.shape[0], n=M.shape[0])
+    Minv = sp.sparse.diags(1./M.diagonal())
 
     # Construct load vector from gravity field
     g = np.zeros(mesh.dims)

python/pbatoolkit/fem.py

@@ -23,47 +23,57 @@
 )
 
 
-def divergence(mesh, quadrature_order: int = 1, GNe: np.ndarray = None):
+def divergence(mesh, quadrature_order: int = 1, eg: np.ndarray = None, wg: np.ndarray = None, GNeg: np.ndarray = None):
     """Construct an FEM divergence operator such that composing div(grad(u)) 
     computes the Laplacian of u, i.e. D @ G = L, where L is the FEM Laplacian 
     operator.
 
     Args:
         mesh (pbat.fem.Mesh): 
         quadrature_order (int, optional): Polynomial order to use for operator construction. Defaults to 1.
-        GNe (np.ndarray, optional): Shape function gradients at quadrature points. 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): 
     """
-    G, GNe = gradient(mesh, quadrature_order=quadrature_order,
-                      GNe=GNe, as_matrix=True)
-    qgL = _fem.inner_product_weights(
-        mesh, quadrature_order=quadrature_order).flatten(order="F")
-    QL = sp.sparse.diags_array([qgL], offsets=[0])
+    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(
+        mesh, eg=eg, GNeg=GNeg, as_matrix=True)
+    wg = _fem.inner_product_weights(
+        mesh, quadrature_order=quadrature_order).flatten(order="F") if should_compute_quadrature else wg
+    QL = sp.sparse.diags_array(np.asarray(wg).squeeze())
     QL = sp.sparse.kron(sp.sparse.eye_array(mesh.dims), QL)
     D = -G.T @ QL
-    return D, GNe
+    return D, eg, wg, GNeg
 
 
-def gradient(mesh, quadrature_order: int = 1, GNe: np.ndarray = None, as_matrix: bool = True):
+def gradient(mesh, quadrature_order: int = 1, eg: np.ndarray = None, GNeg: np.ndarray = None, as_matrix: bool = True):
     """Construct an FEM gradient operator
 
     Args:
         mesh (pbat.fem.Mesh): 
         quadrature_order (int, optional): Polynomial order to use for operator construction. Defaults to 1.
-        GNe (np.ndarray, optional): Shape function gradients at quadrature points. Defaults to None.
+        eg (np.ndarray, optional): Elements corresponding to quadrature points. 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:
         (scipy.sparse.csc_matrix | pbat.fem.Gradient, np.ndarray): 
     """
-    if GNe is None:
-        GNe = _fem.shape_function_gradients(
+
+    if GNeg is None:
+        GNeg = _fem.shape_function_gradients(
             mesh, quadrature_order=quadrature_order)
-    G = _fem.Gradient(mesh, GNe, quadrature_order=quadrature_order)
+        n_quadpts_per_element = GNeg.shape[1] / (mesh.dims * mesh.E.shape[1])
+        eg = np.linspace(0, mesh.E.shape[1]-1,
+                         num=mesh.E.shape[1], dtype=np.int64)
+        eg = np.repeat(eg, n_quadpts_per_element)
+
+    G = _fem.Gradient(mesh, eg, GNeg)
     G = G.to_matrix() if as_matrix else G
-    return G, GNe
+    return G, eg, GNeg
 
 
 def hyper_elastic_potential(
@@ -170,9 +180,8 @@ def mass_matrix(
     if as_matrix:
         M = M.to_matrix()
         if lump:
-            lumpedm = M.sum(axis=0)
-            M = sp.sparse.spdiags(lumpedm, np.array(
-                [0]), m=M.shape[0], n=M.shape[0])
+            lumpedm = np.asarray(M.sum(axis=0)).squeeze()
+            M = sp.sparse.diags_array(lumpedm)
     return M, detJe
 
 
@@ -199,7 +208,7 @@ def load_vector(
     nelems = mesh.E.shape[1]
     wg = np.tile(mesh.quadrature_weights(quadrature_order), nelems)
     qgf = detJe.flatten(order="F") * wg
-    Qf = sp.sparse.diags_array([qgf], offsets=[0])
+    Qf = sp.sparse.diags_array(qgf)
     Nf = _fem.shape_function_matrix(mesh, quadrature_order=quadrature_order)
     if len(fe.shape) == 1:
         fe = np.tile(fe[:, np.newaxis], Qf.shape[0])

source/pbat/fem/Gradient.cpp

@@ -36,16 +36,20 @@ TEST_CASE("[fem] Gradient")
         using Mesh                      = fem::Mesh<Element, kDims>;
         Mesh mesh(V, C);
 
-        MatrixX const GNe  = fem::ShapeFunctionGradients<kQuadratureOrder>(mesh);
-        using GradientType = fem::Gradient<Mesh, kQuadratureOrder>;
+        MatrixX const GNe                = fem::ShapeFunctionGradients<kQuadratureOrder>(mesh);
+        auto const nQuadPointsPerElement = GNe.cols() / (kDims * mesh.E.cols());
+        IndexVectorX const eg = IndexVectorX::LinSpaced(mesh.E.cols(), Index(0), mesh.E.cols() - 1)
+                                    .replicate(1, nQuadPointsPerElement)
+                                    .transpose()
+                                    .reshaped();
+        using GradientType = fem::Gradient<Mesh>;
         CHECK(math::CLinearOperator<GradientType>);
-        GradientType G{mesh, GNe};
+        GradientType G{mesh, eg, GNe};
 
         auto const n                = G.InputDimensions();
         auto const m                = G.OutputDimensions();
         auto const numberOfElements = mesh.E.cols();
-        auto constexpr kQuadPts     = GradientType::QuadratureRuleType::kPoints;
-        CHECK_EQ(m, kDims * kQuadPts * numberOfElements);
+        CHECK_EQ(m, kDims * nQuadPointsPerElement * numberOfElements);
 
         VectorX const ones = VectorX::Ones(n);
         VectorX gradOnes   = VectorX::Zero(m);