Material property as a field

[abhiawasthi1993]

I would like to model an inhomogeneous material where the material property value would be a spatially varying field instead of a constant value throughout (current implementation). Any leads on this would be much appreciated.

[adtzlr]

Hi, do your material properties also change with the deformation?

[abhiawasthi1993]

No, they remain constant with deformation. They are just the function of space at reference configuration.

[adtzlr]

Okay. Do you have linear-elasticity in mind? I'll come up with an example later on.

[abhiawasthi1993]

Right now, I am willing to start with linear elasticity but would like to extend it to non-linear elasticity as well.

[adtzlr]

Here's the example of the README, modified to linear-elasticity and with a spatial-dependent young's modulus.

import felupe as fem

# create a hexahedron-region on a cube
mesh = fem.Cube(n=6)
region = fem.RegionHexahedron(mesh)

# add a field container (with a vector-valued displacement field)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, move=1, clamped=True)

# example: interpolate (initial) X-coordinates from mesh-points to quadrature points
X, Y, Z = fem.Field(region, dim=3, values=mesh.points).interpolate()

# the X-coordinate is now available for every quadrature point and for every cell
print(X.shape) # (8, 125)

# field-dependent young's modulus
E = 1 + X

# define the constitutive material behaviour
# and create a linear-elastic solid body
umat = fem.LinearElastic(E=E, nu=0.3)
solid = fem.SolidBody(umat, field)

# prepare a step with substeps
step = fem.Step(items=[solid], boundaries=boundaries)

# add the step to a job, evaluate
fem.Job(steps=[step]).evaluate()

# visualize the results
view = fem.View(field)
view.plot("Displacement", component=None).show()

Hope that helps!

[adtzlr]

And here is an example for nonlinear elasticity (isotropic hyperelasticity):

import felupe as fem

# create a hexahedron-region on a cube
mesh = fem.Cube(n=6)
region = fem.RegionHexahedron(mesh)

# add a field container (with a vector-valued displacement field)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)

# interpolate (initial) coordinates from mesh-points to quadrature points
X, Y, Z = fem.Field(region, dim=3, values=mesh.points).interpolate()

# define the constitutive material behaviour (with a spatial-dependent shear modulus)
# and create an isotropic-hyperelastic solid body
umat = fem.NeoHooke(mu=1 + X)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)

# prepare a step with substeps
step = fem.Step(
    items=[solid], 
    boundaries=boundaries,
    ramp={boundaries["move"]: fem.math.linsteps([0, 1], num=5)},
)

# add the step to a job, evaluate
job = fem.Job(steps=[step])
job.evaluate()

# visualize the results
view = fem.View(field)
view.plot("Displacement", component=None).show()

[abhiawasthi1993]

Well my plan is to develop an iterative solution to an inverse problem using nonlinear optimization. CasADi server that purpose. I will be using one of the gradient-based optimization modules for the solution. But I suppose the same thing could be done using Tensortrax if I define my parameters accordingly. What are your thoughts on it?

[adtzlr]

I think you could also use scipy.optimize. The thing that is very time-consuming is the jacobian of the unknowns (your material parameters). The other complicated task is to enforce stable material parameters and what to do with numeric errors. If well defined, it's a powerful technique!

[abhiawasthi1993]

Precisely. Can't the evaluation of Jacobian be done using AD only, provided it supports the evaluation of implicit derivatives?

[adtzlr]

Precisely. Can't the evaluation of Jacobian be done using AD only, provided it supports the evaluation of implicit derivatives?

Unfortunately with the current implementation, no.

AD is only used inside the constitutive material formulation in FElupe (sometimes called user material, or short, "umat"). That means for a given hyperelastic strain energy function the gradient and the hessian w.r.t. the deformation gradient are carried out by AD. The whole assembly of the system force vector and the system stiffness matrix are plain sparse arrays built on top of the result of the AD procedure. All AD-information about gradients and jacobians are lost on assembly. That means no AD-information about the jacobian is available on the "top-level" system force-vector and the stiffness matrix.

[adtzlr]

If you're fine with a more math-centric approach, scikit-fem could have that implemented in skfem.experimental.autodiff.NonLinearForm. However, you would have to ask your questions in their repo because I don't remember how to use it correctly.

[adtzlr]

Hi @abhiawasthi1993, sorry for the super-late reply. I don't know if you're still interested in this topic.

Here's a quick example where a material parameter is defined by a field and the strain energy density function of a hyperelastic solid is defined by matadi. However, I'm still not sure how to create the Jacobian (using automatic differentiation) w.r.t. the material parameter(s). I think this depends on the quantity you would like to fit, I could imagine the maximum principal log. strains or the displacement vectors. Anyway, an additional field would be required for this quantity such that the Jacobian can be extracted from the list of matrices = solid.assemble.matrix(block=False) (this feature is currently unreleased, requires latest git main branch). I could extend this example accordingly if you're still interested. Please let me know!

import felupe as fem
import numpy as np
import matadi as mat
import matadi.math as mm

mesh = fem.Rectangle(n=3)
region = fem.RegionQuad(mesh)

# a mixed-field container with displacements and material parameter field
field = fem.FieldContainer(
    fields=[fem.FieldPlaneStrain(region, dim=2), fem.Field(region, dim=1)]
)


def fun(x, K=5):
    # unpack input variables
    F, C10, ζ = x[:3]

    C = F.T @ F
    J = mm.det(F)
    I1 = mm.trace(C) * J ** (-2 / 3)

    # strain energy density function
    ψ = C10 * (I1 - 3) + K * (J - 1) ** 2 / 2

    # update of state variables (not used)
    ζ_new = ζ

    return [mm.gradient(ψ, F), mm.gradient(ψ, C10), ζ_new]


# symbolic variables, including (required - but not used) state variables
F = mat.Variable("F", 3, 3)
C10 = mat.Variable("C10")
ζ = mat.Variable("ζ")

umat = mat.MaterialTensor([F, C10, ζ], fun, statevars=1)
solid = fem.SolidBody(umat=umat, field=field)

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)

# prescribed values for the material parameters by a boundary condition
entire_model = np.ones(field[-1].region.mesh.npoints, dtype=bool)
boundaries["C10"] = fem.Boundary(field[-1], mask=entire_model)
boundaries["C10"].update(value=0.5 + np.random.rand(field[-1].region.mesh.npoints) / 10)

move = fem.math.linsteps([0, 1], num=5)
ramp = {boundaries["move"]: move}
step = fem.Step(items=[solid], ramp=ramp, boundaries=boundaries)

job = fem.Job(steps=[step]).evaluate()
ax = solid.imshow("Principal Values of Logarithmic Strain")