simkit.deformation_jacobian_p2

Deformation-Jacobian operator for quadratic (P2) simplicial elements.

Functions

deformation_jacobian_p2(→ scipy.sparse.csc_matrix)

Sparse operator J mapping P2 nodal positions to per-cubature F.

Module Contents

simkit.deformation_jacobian_p2.deformation_jacobian_p2(V2: numpy.ndarray, T2: numpy.ndarray, bary: numpy.ndarray, weights: numpy.ndarray = None) → scipy.sparse.csc_matrix

Sparse operator J mapping P2 nodal positions to per-cubature F.

Builds the constant linear map such that F = (J @ x).reshape(-1, dim, dim) gives the deformation gradient at every cubature point of every element, where x is V2 flattened in C order.

Is this matrix constant?

Yes. Like the P1 operator it is a fixed sparse matrix, built once from the rest geometry and the (fixed) reference-space quadrature points, and is independent of the deformed state. The key difference from P1: P2 shape-function gradients are linear in the reference coordinates, so F is no longer constant within an element – it varies from one cubature point to the next. That is why J maps to a stack of per-cubature-point F blocks (t * n_quad of them) rather than one F per element. Because the P2 midpoint nodes sit exactly at edge midpoints, the geometric rest map is still affine, so the reference->rest Jacobian is constant per element and is taken from the corner nodes alone.

How P2 changes the integration of an energy

With P1, each element has a single constant F, so the total elastic energy is a one-point quadrature:

E = Σ_t vol_t · ψ(F_t)

With P2, F varies inside the element, so the integral becomes a sum over the n_quad cubature points:

E = Σ_t Σ_q w_{t,q} · ψ(F_{t,q})

No energy code changes are needed – treat each cubature point as a pseudo-element. Using this operator together with the weights from gauss_legendre_quadrature():

x = V2.reshape(-1, 1)
F = (J @ x).reshape(-1, dim, dim)      # (t*n_quad, dim, dim)
psi = some_energy_element_F(F, mu)     # unchanged P1 energy code
E = float((weights.reshape(-1, 1) * psi).sum())

and likewise forces are J.T @ (weights ⊙ P) and the Hessian is J.T @ blockdiag(weights ⊙ d2psi) @ J. The quadrature order used to build bary/weights must be high enough for the chosen energy’s integrand (e.g. degree 2 for linear elasticity, whose density is quadratic in F).

param V2:

Quadratic-mesh vertex positions (dim is 2 or 3).

type V2:

np.ndarray (n2, dim)

param T2:

Quadratic connectivity (6 nodes for triangles, 10 for tets), as produced by linear_to_quadratic_elements().

type T2:

np.ndarray (t, n_nodes)

param bary:

Barycentric coordinates of the cubature points (from gauss_legendre_quadrature()).

type bary:

np.ndarray (t, n_quad, dim+1)

param weights:

Cubature weights. Accepted for interface symmetry but not used to build the operator: the operator is independent of the weights, which enter later as the energy’s per-cubature vol.

type weights:

np.ndarray (t, n_quad), optional

returns:

J – Constant deformation-Jacobian operator.

rtype:

scipy.sparse.csc_matrix (t*n_quad*dim*dim, n2*dim)