[QUESTION] Understanding warp.fem.deformation_gradient behaviour in deformed configurations

[akrivx]

Hi,

I'm implementing corotational FEM structural dynamics using warp.fem, and I'd like to clarify how fem.deformation_gradient works in the context of a deformed configuration.

My base geometry is a 3D grid, and at each simulation step, I create a deformed geometry using:

deformed_geo = u_field.make_deformed_geometry() # u_field is the displacement field (vec3d)
deformed_domain = fem.Cells(deformed_geo)

Since u_field.dof_values stores relative displacements from the base configuration, I call make_deformed_geometry using the default relative=True parameter.

I then compute per-element rotation matrices using fem.interpolate:

# in each step:
fem.interpolate(
  element_rotation_matrix_integrand,
  domain=deformed_domain,
  quadrature=fem.RegularQuadrature(deformed_domain, 0),
  values={'Re': Re} # output array of wp.mat33 with one slot per element
)

Inside element_rotation_matrix_integrand, I call:

F = fem.deformation_gradient(domain, s)

I expected F to represent the deformation gradient relative to the reference configuration, but in my case, F is always a diagonal matrix with the diagonal entries corresponding to element dimensions.

I want to make sure that I correctly compute the deformation gradient relative to the base configuration while working with deformed geometries in warp.fem.

Would you be able to clarify how this should be handled? Am I using deformation_gradient wrong? Thanks in advance!

Regards,
Chris

[gdaviet]

Hi Chris,

fem.deformation_gradient computes the deformation gradient of the domain with respect to the reference element configuration. For a DeformedGeometry, this corresponds to deformation_gradient(deformed_geo(u), s) = (I + grad(u), s) @ deformation_gradient(base_geo), s); see the script below. So in your case, I expect the diagonal scaling comes from the deformation_gradient(base_geo), s) part.

If you're simply looking to compute the deformation gradient of the displacement field with respect to the base geometry, then you can simply evaluate (I + grad(u), s) over the base geometry, no need to explicitly create a deformed geometry. See for instances examples/fem/example_mixed_elasticity.py.

import numpy as np

import warp as wp
import warp.fem as fem


@fem.integrand
def def_grad_field(s: fem.Sample, domain: fem.Domain, u: fem.Field):
    F = wp.identity(n=3, dtype=float) + fem.grad(u, s)
    return F @ fem.deformation_gradient(domain, s)


@fem.integrand
def def_grad(s: fem.Sample, domain: fem.Domain):
    return fem.deformation_gradient(domain, s)


grid = fem.Grid3D(res=wp.vec3i(2))

# Create a shear field
u_field = fem.make_polynomial_space(grid, dtype=wp.vec3).make_field()
u = u_field.dof_values.numpy()
u[:, 0] = u_field.space.node_positions().numpy()[:, 1]
u_field.dof_values = u

# compute deformation from grid over non-deformed grid
quadrature = fem.RegularQuadrature(fem.Cells(grid), order=1)
result_undef = wp.empty(quadrature.total_point_count(), dtype=wp.mat33)
fem.interpolate(def_grad_field, quadrature=quadrature, dest=result_undef, fields={"u": u_field})

print(result_undef.numpy())

# compute deformation from deformed grid 
deformed_grid = u_field.make_deformed_geometry(relative=True)
quadrature = fem.RegularQuadrature(fem.Cells(deformed_grid), order=1)
result_def = wp.empty(quadrature.total_point_count(), dtype=wp.mat33)
fem.interpolate(def_grad, quadrature=quadrature, dest=result_def)

print(result_def.numpy())

assert np.max(np.abs((result_def.numpy() - result_undef.numpy()))) < 1.e-8

[akrivx]

Hi @gdaviet ,

Excellent, thank you.

Regards,
Chris