Define element-wise thermal expansion coefficient

variational formulation
1

I’m trying to solve a thermoelastic problem with constant delta_T but with a different alpha for each element.
When all elements are the same, the code works as I define
f = fem.Constant(domain, default_scalar_type((kappa, kappa, kappa))), where kappa = delta_T * (alpha*(2*mu + 3*lambda_))

As I understand, all I have to do is to define f as a vector function with 3 values per element, but I only managed to create kappa_local with a single value per element. I’m struggling hard with the documentation and all the examples/questions that I find seem to be outdated.

The code:

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

# Define the box domain dimensions
Lx, Ly, Lz = 1.0, 1.0, 1.0  # Length in x, y, z directions
nx, ny, nz = 10, 10, 10  # Number of elements in x, y, z directions

# Create the box mesh
domain = mesh.create_box(MPI.COMM_WORLD,
                                 [np.array([0.0, 0.0, 0.0]), np.array([Lx, Ly, Lz])],
                                 [nx, ny, nz],
                                 mesh.CellType.tetrahedron)

# Define the function space
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))

# Define boundary condition
def clamped_boundary(x):
    return np.isclose(x[0], 0)

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

# Define material properties
E = 10.0  # Young's modulus
nu = 0.3  # Poisson's ratio
alpha = 1e-5  # Coefficient of thermal expansion
delta_T = -100  # Temperature change

lambda_ = (E * nu / ((1 + nu) * (1 - 2 * nu)))
mu = (E / (2 * (1 + nu)))
kappa = delta_T * (alpha*(2*mu + 3*lambda_))

print(f"lambda = {lambda_}, mu = {mu}, kappa = {kappa}")

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

def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u) 

V0 = fem.functionspace(domain, ("DG", 0))
kappa_local = fem.Function(V0)

# Assign different kappa values to the elements
num_dofs = V.dofmap.index_map.size_local
# Replace the following with your actual method for generating kappa values
kappa_values = (kappa * np.random.rand(num_dofs)  - 0.5) * 2
#kappa_values = np.tile(kappa_values[:, np.newaxis], (1, 3))

for i in range(num_dofs):
    kappa_local.vector.setValue(i, kappa_values[i])

# Define the variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
#f = fem.Constant(domain, default_scalar_type((kappa, kappa, kappa)))
f = kappa_local
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx + ufl.dot(T, v) * ufl.Measure("ds", domain=domain) 

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

f is equal to kappa local, which is a scalar valued function, which you cannot take the inner product with v with in

3

Hi @dokken, yes, I’m sorry my question was not clear. I understand why it fails, but I do not know how to define kappa_local as a vector function (with 3 components) so that the inner product works.

4

I think I made it work. It turns out kappa_local was (after creating an appropiate function space in 3D) a vector of size (ndof * 3, ) instead of (ndof , 3).

V0 = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
kappa_local = fem.Function(V0)
ndof = V0.dofmap.index_map.size_local
kappa_values = (kappa * np.random.rand(ndof * 3)  - 0.5) * 5
for i in range(ndof * 3):
    kappa_local.vector.setValue(i, kappa_values[i])