Inhomogeneous Stiffness Tensor UFL error

Hello everyone,

I try to solve the problem of linear elasticity on a simple domain with spatially varying stiffness tensor C(x) (non-isotropic). So far, I insert a scalar valued function f_C(x) into the first three main diagonal entries of C(x). This works perfectly fine if f_C(x) does not contain any if/else statements which I need to assign different values to C(x) in different parts of the domain. If f_C(x) does contain if/else statements I always get the following error message:
“ValueError: UFL conditions cannot be evaluated as bool in a Python context.”

My code example is the following:

# imports
from dolfinx import mesh, fem, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

# create domain
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))


# define boundary conditions+
def clamped_boundary(x):
    bool_x = np.isclose(x[0], 0)
    bool_y = np.isclose(x[1], 2)
    return np.logical_or(bool_x, bool_y)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

# define measure on boundary
ds = ufl.Measure("ds", domain=domain)

# create anisotropic and inhomogeneous stiffness tensor C(x)
X = ufl.SpatialCoordinate(domain)
def f_C(x):
    if x[0] <= 0.5:
        return 1
    else:
        return x[0]*10

# insert some dummy values   
C = ufl.as_tensor([[1 + f_C(X), 1.05, 0.96, 0,0,0],
                   [1.05, 2 + 0.1*f_C(X) , 0.51,0,0,0],
                   [0.96,0.51, 2.86  + 0.1*f_C(X),0,0,0],
                   [0,0,0,3.15,0,0],
                   [0,0,0,0,3.35,0],
                   [0,0,0,0,0,4.1]])

# define helper function that switch from voigt to matrix notation
def eps_to_voigt(e):
    return ufl.as_vector([e[0,0], e[1,1], e[2,2], 2*e[1,2], 2*e[0,2], 2*e[0,1]])

def voigt_to_stress(s):
    return ufl.as_tensor([[s[0], s[5], s[4]],
                          [s[5], s[1], s[3]],
                          [s[4], s[3], s[2]]])

def epsilon(u):
    return ufl.sym(ufl.grad(u))  

def sigma(u, C):
    return voigt_to_stress( ufl.dot(C, eps_to_voigt(epsilon(u))) ) 

# variational formulation
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# body force f
f = fem.Constant(domain, default_scalar_type((0.0, -0.02, 0)))
# tranction force T
T = fem.Constant(domain, default_scalar_type((0.0, 0.0, 0)))

# variational problem
a = ufl.inner(sigma(u, C), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

# solve PDE
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

Thanks for any kind of help :slight_smile:

Use ufl.conditional, as shown in

Thanks @dokken , that works perfectly :+1: .

One more question, is there a more efficient and elegant way to input an arbitrary stiffness tensor C(x) into the equations of linear elasticity, that explicitly depends on the spatial variable x? All I can find in the tutorials are examples for isotropic, homogeneous cases where C(x) only depends on the Lamé coefficients.

Thank you.

Do you mean another way than using spatial coordinate and condtitionals for the various components?

You could interpolate C(x) into an appropriate tensor function space. Then you can use if/elese statements in the Python function you pass to the interpolate-call.

Yes exactly, assume I have a piecewise constant field for my stiffness tensor given by two (or more) tensors C_1, C_2 (assuming orthotropic material so C_i is a 6x6 matrix).
All I could find to interpolate a Python function onto a suitable subspace is similar to the example below which raises: " ValueError: setting an array element with a sequence."

space_C = fem.functionspace(domain, ("Lagrange", 1, (6,6)))
C_hat = fem.Function(space_C)

def C(x, tensor_list = [c_1,c_2]):
    # x is assumed to have shape (3, num_points)
    C = np.zeros((6, 6, x[0].size))
    for idx, x_coord in enumerate(x[0]):
        if x_coord <= 1:
            C[:,:, idx] = tensor_list[0]
        else:
            C[:,:, idx] = tensor_list[1]
   
    return ufl.as_tensor(C)

C_hat.interpolate(C)

I am relatively new to FeniCS and therefore struggling to find a suitable way of defining a general Python function that takes my spatial variable x of shape (3, num_points) and returns C(x) of shape (6,6, num_points) (even not sure if that is the correct shape) to be used inside the interpolate-call.

Here is a minimal example for a 6x6 tensor that covers various conditionals:

from mpi4py import MPI
import dolfinx
import numpy as np

mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)

shape = (6, 6)
V = dolfinx.fem.functionspace(mesh, ("DG", 1, shape))

value1 = 1.3
value2 = 0.8
value3 = 0.2
tol = 1e-12


def f(x):
    values = np.zeros((shape[0], shape[1], x.shape[1]), dtype=np.float64)
    # Add value1 to the tensor at (0,0) wherever x<=1
    cond1 = x[0] <= 1 + tol
    values[0, 0] += value1 * cond1

    # Add value2 to the tensor at (2,3) wherever y < 0.5 and z<0.4
    cond2 = (x[1] < 0.5 + tol) & (x[2] < 0.4 + tol)
    values[2, 3] += value2 * cond2
    # Add value3 to the tensor at (4,3)
    values[4, 3] += value3
    return values.reshape(-1, x.shape[1])


u = dolfinx.fem.Function(V)
u.interpolate(f)


with dolfinx.io.VTXWriter(mesh.comm, "tensor.bp", [u]) as bp:
    bp.write(0)

Thanks @dokken that really helped :+1: