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