Can this be used for parameter identification?

[mattragoza]

I have observed displacement data and I want to recover the elasticity of the material by optimizing the elasticity parameter field through the forward model. Can this package compute sensitivity with respect to the material parameters for optimization?

[meyer-nils]

Hi Matt,

Absolutely. We already used this package to optimize fiber orientation fields by minimizing the error between simulation and observed displacement fields [1].

Is your material model simply linear isotropic elasticity and you want to find E and nu in each element?

Nils

[1] https://opus.bibliothek.uni-augsburg.de/opus4/files/114634/114634.pdf

[mattragoza]

Thanks for the reply and reference. Yes, we are assuming isotropic linear elasticity and want to recover spatially inhomogeneous E. We are using a DL model to estimate a map of E and drive the simulation, so we need to solve the FEM model many times during training, with automatic differentiation, and preferably with GPU acceleration.

[meyer-nils]

Okay, you may try something along these lines:

from tqdm import tqdm

# Initial stiffness guess
E0 = 1.0
# Let's assume that Poisson ratio is constant
nu = 0.3

# Set up intial model
C0 = Isotropic(E=E0, nu=0.3).C()
model = import_mesh("<your_mesh>", C=C0)

def compute_loss(E):
    lbd = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
    mu = E / (2.0 * (1.0 + nu))
    z = torch.zeros_like(lbd)
    model.C = torch.stack(
        [
            torch.stack([lbd + 2.0 * mu, lbd, lbd, z, z, z], dim=-1),
            torch.stack([lbd, lbd + 2.0 * mu, lbd, z, z, z], dim=-1),
            torch.stack([lbd, lbd, lbd + 2.0 * mu, z, z, z], dim=-1),
            torch.stack([z, z, z, mu, z, z], dim=-1),
            torch.stack([z, z, z, z, mu, z], dim=-1),
            torch.stack([z, z, z, z, z, mu], dim=-1),
        ],
        dim=-1,
    )
    u, _ = model.solve()
    return ((u-u_target)**2).sum()

# Optimize
E = E0 * torch.ones(len(model.elements))
E.requires_grad = True
optimizer = torch.optim.Adam([E], lr=0.1)

for _ in tqdm(range(100)):
    optimizer.zero_grad()
    objective = compute_loss(E)
    objective.backward()
    optimizer.step()

If you interpolate your E from a neural network written in pytorch, you should be able to seamlessly backpropagate the error through torchfem and the network. I did not experiment too much with GPU acceleration, but it should work in principle, too.

[mattragoza]

Thank you, this is great. Slightly offtopic question, but is there a way to evaluate the FEM displacement field at arbitrary points, for instance if I wanted to render it like an image?

[meyer-nils]

There is no built-in function for that yet. However, it should be a straightforward task:

1. Determine the element, in which the query point is located (e.g. by finding the closest element center).
2. Use that elements shape functions to evaluate the displacement at the query point.