Handling Function returns from variational equations

General
1

Hey all,

How are Function values handled when solving the variational equation?

For example if we consider the Hyperelasticty Tutorial:

v = ufl.TestFunction(V)
u = fem.Function(V)

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

Is it possible to write something similar to the following that works:

class MaterialCoordinates():
        # (0): Defining rotation
        def __init__(self, rot=ROT):
            self.rot = ROT
        # (1): Rotate coordinates into new field
        def eval(self, x):
            self.xMat = np.cos(self.rot) * x[0] + np.sin(self.rot) * x[1]
            self.yMat = -np.sin(self.rot) * x[0] + np.sin(self.rot) * x[1]
            self.zMat = x[2]
            return (self.xMat, self.yMat, self.zMat)

matCoords = MaterialCoordinates()
u_mat = matCoords.eval(x=u)
F = ufl.variable(I + ufl.grad(u_mat))

Before the earlier lines so that the deformation gradient is with respect to a rotated coordinate system? i.e. how can I interact with outputs of u on each iteration to change calculation of F?

2

Any reason you can’t use something like

x = ufl.SpatialCoordinate(mesh)
rot = np.pi / 2
mat_x = ufl.cos(rot) * x[0]
mat_y = ufl.sin(rot) * x[1]
mat_z = x[2]
mat = ufl.as_vector((mat_x, mat_y, mat_z))

?

If you can’t use it, please post a complete code (ref Read before posting: How do I get my question answered?), because no-one can run what you wrote, and thus I couldn’t even test the snippet that I sent you.

3

Hey @francesco-ballarin, I appreciate the response. What you provided is fine for a ufl.SpatialCoordinte(mesh) but my question is regarding u (u = fem.Function(V)) as the output of the Hyperelasticity Tutorial.

I’m interested in if I can conduct a similar operation to my snippet (or yours) on u to convert it into a different coordinate scheme before calculating F or other kinematic variables (i.e. C, E).

The code of interest is the tutorial code, but per completeness here is a MWE edited copy with my snippet inserted which produces an error:

from dolfinx import log, default_scalar_type, io
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import matplotlib.pyplot as plt
import pyvista
from dolfinx import nls
import numpy as np
import ufl

from mpi4py import MPI
from dolfinx import fem, mesh, plot
L = 20.0
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("Lagrange", 2))

def left(x):
    return np.isclose(x[0], 0)


def right(x):
    return np.isclose(x[0], L)


fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

###########################
class MaterialCoordinates():
        # (0): Defining rotation
        def __init__(self, rot=np.pi/4):
            self.rot = np.pi/4
        # (1): Rotate coordinates into new field
        def eval(self, x):
            self.xMat = np.cos(self.rot) * x[0] + np.sin(self.rot) * x[1]
            self.yMat = -np.sin(self.rot) * x[0] + np.sin(self.rot) * x[1]
            self.zMat = x[2]
            return (self.xMat, self.yMat, self.zMat)

matCoords = MaterialCoordinates()
u_mat = matCoords.eval(x=u)
F = ufl.variable(I + ufl.grad(u_mat))
###########################

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

problem = NonlinearProblem(F, u, bcs)

solver = NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

num_its, converged = solver.solve(u)
if converged:
    print(f"Converged in {num_its} iterations.")
else:
    print(f"Not converged after {num_its} iterations.")

# += Export
with io.VTXWriter(MPI.COMM_WORLD, "test.bp", u, engine="BP4") as vtx:
    vtx.write(0.0)
    vtx.close()

WARNING:py.warnings:/Users/murrayla/Documents/main_PhD/FEniCSx/FEniCSx/Material_Field/test.py:14: DeprecationWarning: This method is deprecated. Use FunctionSpace with an element shape argument instead
  V = fem.VectorFunctionSpace(domain, ("Lagrange", 2))

Traceback (most recent call last):
  File "/Users/murrayla/Documents/main_PhD/FEniCSx/FEniCSx/Material_Field/test.py", line 65, in <module>
    F = ufl.variable(I + ufl.grad(u_mat))
                         ^^^^^^^^^^^^^^^
  File "/Users/murrayla/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/operators.py", line 345, in grad
    f = as_ufl(f)
        ^^^^^^^^^
  File "/Users/murrayla/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/constantvalue.py", line 501, in as_ufl
    raise ValueError(
ValueError: Invalid type conversion: (Sum(Product(FloatValue(0.7071067811865476), Indexed(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, hexahedron, 1, gll_warped, unset, False), (3,)), 0), blocked element (Basix element (P, hexahedron, 2, gll_warped, unset, False), (3,))), 0), MultiIndex((FixedIndex(0),)))), Product(FloatValue(0.7071067811865475), Indexed(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, hexahedron, 1, gll_warped, unset, False), (3,)), 0), blocked element (Basix element (P, hexahedron, 2, gll_warped, unset, False), (3,))), 0), MultiIndex((FixedIndex(1),))))), Sum(Product(FloatValue(-0.7071067811865475), Indexed(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, hexahedron, 1, gll_warped, unset, False), (3,)), 0), blocked element (Basix element (P, hexahedron, 2, gll_warped, unset, False), (3,))), 0), MultiIndex((FixedIndex(0),)))), Product(FloatValue(0.7071067811865475), Indexed(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, hexahedron, 1, gll_warped, unset, False), (3,)), 0), blocked element (Basix element (P, hexahedron, 2, gll_warped, unset, False), (3,))), 0), MultiIndex((FixedIndex(1),))))), Indexed(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, hexahedron, 1, gll_warped, unset, False), (3,)), 0), blocked element (Basix element (P, hexahedron, 2, gll_warped, unset, False), (3,))), 0), MultiIndex((FixedIndex(2),)))) can not be converted to any UFL type.
4

Then just use u instead of x = ufl.SpatialCoordinate(mesh) in my snippet above.

1 Like
5

Hey @francesco-ballarin, great. much more elegant than I expected. Thank you.