Allow dynamic quadrature in HyperElasticPotential

Author: Q-Minh

bindings/pypbat/fem/HyperElasticPotential.cpp

@@ -52,8 +52,8 @@
 
     # Create hyper elastic potential
     Y, nu, energy = args.Y, args.nu, pbat.fem.HyperElasticEnergy.StableNeoHookean
-    hep, detJeU, GNeU = pbat.fem.hyper_elastic_potential(
-        mesh, Y=Y, nu=nu, energy=energy, detJe=detJeF)
+    hep, egU, wgU, GNeU = pbat.fem.hyper_elastic_potential(
+        mesh, Y=Y, nu=nu, energy=energy)
 
     # Set Dirichlet boundary conditions
     Xmin = mesh.X.min(axis=1)

python/examples/elasticity.py

@@ -53,7 +53,7 @@
 
     # Create hyper elastic potential
     Y, nu, psi = args.Y, args.nu, pbat.fem.HyperElasticEnergy.StableNeoHookean
-    hep, detJeU, GNeU = pbat.fem.hyper_elastic_potential(
+    hep, egU, wgU, GNeU = pbat.fem.hyper_elastic_potential(
         mesh, Y=Y, nu=nu, energy=psi, quadrature_order=4)
 
     # Set Dirichlet boundary conditions

python/examples/elasticity_higher_order.py

@@ -378,8 +378,8 @@ def __call__(self, xk: np.ndarray):
 
     # Create hyper elastic potential
     Y, nu, psi = args.Y, args.nu, pbat.fem.HyperElasticEnergy.StableNeoHookean
-    hep, detJeU, GNeU = pbat.fem.hyper_elastic_potential(
-        mesh, Y=Y, nu=nu, energy=psi, detJe=detJeF)
+    hep, egU, wgU, GNeU = pbat.fem.hyper_elastic_potential(
+        mesh, Y=Y, nu=nu, energy=psi)
 
     # Setup IPC contact handling
     F = igl.boundary_facets(C)

python/examples/ipc.py

@@ -37,7 +37,7 @@ def signal(w: float, v: np.ndarray, t: float, c: float, k: float):
     x = mesh.X.reshape(math.prod(mesh.X.shape), order='f')
     M, detJeM = pbat.fem.mass_matrix(mesh, rho=args.rho)
     Y, nu, energy = args.Y, args.nu, pbat.fem.HyperElasticEnergy.StableNeoHookean
-    hep, detJeU, GNeU = pbat.fem.hyper_elastic_potential(
+    hep, egU, wgU, GNeU = pbat.fem.hyper_elastic_potential(
         mesh, Y, nu, energy=energy)
     hep.compute_element_elasticity(x)
     U, gradU, HU = hep.eval(), hep.gradient(), hep.hessian()

python/examples/modes.py

@@ -148,7 +148,7 @@ def __init__(self, MD: pbat.fem.Mesh, MS: pbat.fem.Mesh, MT: pbat.fem.Mesh, H, l
         self.P = (self.NT.T @ (rho * self.Ig) @ self.NS).asformat('csc')
         self.rho = rho
         energy = pbat.fem.HyperElasticEnergy.StableNeoHookean
-        self.hep, self.detJeU, self.GNeU = pbat.fem.hyper_elastic_potential(
+        self.hep, self.egU, self.wgU, self.GNeU = pbat.fem.hyper_elastic_potential(
             MT, Y, nu, energy=energy)
         self.hxreg = hxreg
         self.X = MT.X
@@ -219,7 +219,7 @@ def __matmul__(self, B):
 def rest_pose_hessian(mesh, Y, nu):
     x = mesh.X.reshape(math.prod(mesh.X.shape), order='f')
     energy = pbat.fem.HyperElasticEnergy.StableNeoHookean
-    hep, detJeU, GNeU = pbat.fem.hyper_elastic_potential(
+    hep, egU, wgU, GNeU = pbat.fem.hyper_elastic_potential(
         mesh, Y, nu, energy=energy)
     hep.compute_element_elasticity(x)
     HU = hep.hessian()

python/examples/multigrid/linear_restriction_prolongation.py

@@ -16,7 +16,9 @@
 with contextlib.redirect_stdout(_strio):
     help(_fem)
 _strio.seek(0)
-setattr(__module, "__doc__", f"{getattr(__module, "__doc__")}\n\n{_strio.read()}")
+setattr(__module, "__doc__", f"{
+        getattr(__module, "__doc__")}\n\n{_strio.read()}")
+
 
 def divergence(mesh, quadrature_order: int = 1, GNe: np.ndarray = None):
     """Construct an FEM divergence operator such that composing div(grad(u)) 
@@ -68,8 +70,9 @@ def hyper_elastic_potential(
         energy=_fem.HyperElasticEnergy.StableNeoHookean,
         precompute_sparsity: bool = True,
         quadrature_order: int = 1,
-        detJe: np.ndarray = None,
-        GNe: np.ndarray = None):
+        eg: np.ndarray = None,
+        wg: np.ndarray = None,
+        GNeg: np.ndarray = None):
     """Construct an FEM hyper elastic potential
 
     Args:
@@ -79,23 +82,27 @@ def hyper_elastic_potential(
         energy (pbat.fem.HyperElasticEnergy, optional): Constitutive model. Defaults to pbat.fem.HyperElasticEnergy.StableNeoHookean.
         precompute_sparsity (bool, optional): Precompute an acceleration data structure for fast hessian construction. Defaults to True.
         quadrature_order (int, optional): Polynomial order to use for potential (and its derivatives) evaluation. 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): Elements corresponding to each quadrature weight in wg.
+        wg (np.ndarray, optional): Quadrature weights at quadrature points. Defaults to None.
+        GNeg (np.ndarray, optional): Shape function gradients at quadrature points. Defaults to None.
 
     Returns:
         (pbat.fem.HyperElasticPotential, 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)
     hep = _fem.HyperElasticPotential(
-        mesh, detJe, GNe, Y, nu, energy=energy, quadrature_order=quadrature_order)
+        mesh, eg, wg, GNeg, Y, nu, energy=energy)
     if precompute_sparsity:
         hep.precompute_hessian_sparsity()
-    return hep, detJe, GNe
+    return hep, eg, wg, GNeg
 
 
 def laplacian(

python/pbatoolkit/fem.py

@@ -38,20 +38,23 @@ TEST_CASE("[fem] HyperElasticPotential")
             else
                 return kOrder - 1;
         }();
-        using ElasticEnergyType = physics::StableNeoHookeanEnergy<3>;
-        using ElementType       = fem::Tetrahedron<kOrder>;
-        using MeshType          = fem::Mesh<ElementType, kDims>;
-        using ElasticPotentialType =
-            fem::HyperElasticPotential<MeshType, ElasticEnergyType, kQuadratureOrder>;
-        using QuadratureType = ElasticPotentialType::QuadratureRuleType;
+        using ElasticEnergyType    = physics::StableNeoHookeanEnergy<3>;
+        using ElementType          = fem::Tetrahedron<kOrder>;
+        using MeshType             = fem::Mesh<ElementType, kDims>;
+        using ElasticPotentialType = fem::HyperElasticPotential<MeshType, ElasticEnergyType>;
 
         CHECK(math::CLinearOperator<ElasticPotentialType>);
 
         MeshType const M(V, C);
-        VectorX const x     = M.X.reshaped();
-        MatrixX const detJe = fem::DeterminantOfJacobian<kQuadratureOrder>(M);
-        MatrixX const GNe   = fem::ShapeFunctionGradients<kQuadratureOrder>(M);
-        ElasticPotentialType U(M, detJe, GNe, x, Y, nu);
+        VectorX const x       = M.X.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.rows())
+                                    .transpose()
+                                    .reshaped();
+        ElasticPotentialType U(M, eg, wg, GNeg, x, Y, nu);
         Scalar const UMaterial      = U.Eval();
         VectorX const gradUMaterial = U.ToVector();
         CSCMatrix const HMaterial   = U.ToMatrix();