Inhomogeneous Stiffness Tensor UFL error

Errors
,
1

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:

2

Use ufl.conditional, as shown in

3

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.

4

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.

5

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.

6

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)
7

Thanks @dokken that really helped :+1: