Modeling heterogeneous material expansion with mixed elements

[cborau]

I’m trying to model the heterogeneous dilation of a 3D solid, composed of a two-phase material that has different properties. For that, I’m creating mixed elements P and Q with different mechanical properties and expansion coefficients.
I have several questions.
1- I’m not sure I have formulated correctly the variatonal formulation. As I understand it (or need it), each P element should share the nodes with its Q counterpart.
2- I don’t know how to stablish the dirichlet boundary conditions as I’m struggling with the nomenclature. I understand that I need to define u_D as a 2x3 vector of zeros, but I’m not able to do so.

This is the code I have so far:

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner, sym, nabla_div, Identity, grad
from dolfinx import default_scalar_type

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
P_el = element("Lagrange", domain.basix_cell(), 1, shape=(3,))
Q_el = element("Lagrange", domain.basix_cell(), 1, shape=(3,))
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

# Define material properties for P elements
E_p = 10.0  # Young's modulus
nu_p = 0.3  # Poisson's ratio
alpha_p = 1e-5  # Coefficient of thermal expansion
delta_T_p = -100  # Temperature change

lambda_p = (E_p * nu_p / ((1 + nu_p) * (1 - 2 * nu_p)))
mu_p = (E_p / (2 * (1 + nu_p)))
kappa_p = delta_T_p * (alpha_p * (2 * mu_p + 3 * lambda_p))

# Define material properties for P elements
E_q = 20.0  # Young's modulus
nu_q = 0.4  # Poisson's ratio
alpha_q = 2e-5  # Coefficient of thermal expansion
delta_T_q = -100  # Temperature change

lambda_q = (E_q * nu_q / ((1 + nu_q) * (1 - 2 * nu_q)))
mu_q = (E_q / (2 * (1 + nu_q)))
kappa_q = delta_T_q * (alpha_q * (2 * mu_q + 3 * lambda_q))


(u_p, u_q) = TrialFunctions(V)
(v_p, v_q) = TestFunctions(V)

x = SpatialCoordinate(domain)
# f_p = ? # list of individual kappa values for P elements
# f_q = ? # list of individual kappa values for Q elements
f_p = fem.Constant(domain, default_scalar_type((kappa_p, kappa_p, kappa_p)))
f_q = fem.Constant(domain, default_scalar_type((kappa_q, kappa_q, kappa_q)))

dx = Measure("dx", domain)

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

def sigma(u):
    return lambda_p * nabla_div(u) * Identity(len(u)) + 2 * mu_p * epsilon(u) + lambda_q * nabla_div(u) * Identity(len(u)) + 2 * mu_q * epsilon(u) 


a = inner(sigma(u_p), epsilon(v_p)) * dx
L = inner(f_p, v_p) * dx  + inner(f_q, v_q) * dx 

# Get subspace of V
V0 = V.sub(0)
fdim = domain.topology.dim - 1
facets_x0 = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[0], 0.0))
Q, _ = V0.collapse()
dofs_x0 = fem.locate_dofs_topological((V0, Q), fdim, facets_x0)

u_D = np.array([0, 0, 0], dtype=default_scalar_type) # THIS IS WRONG
bc = fem.dirichletbc(u_D, dofs_x0, V0)

If ran like that, it throws this error:

Exception has occurred: TypeError
__init__(): incompatible function arguments. The following argument types are supported:
    1. __init__(self, g: ndarray[dtype=float64, writable=False, order='C'], dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C'], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
    2. __init__(self, g: dolfinx.cpp.fem.Constant_float64, dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C'], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
    3. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    4. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: list[ndarray[dtype=int32, writable=False, shape=(*), order='C']], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None

Invoked with types: dolfinx.cpp.fem.DirichletBC_float64, ndarray, list, dolfinx.cpp.fem.FunctionSpace_float64
TypeError: __init__(): incompatible function arguments. The following argument types are supported:
    1. __init__(self, g: ndarray[dtype=float64, writable=False, order='C'], dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C'], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
    2. __init__(self, g: dolfinx.cpp.fem.Constant_float64, dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C'], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
    3. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    4. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: list[ndarray[dtype=int32, writable=False, shape=(*), order='C']], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None

Invoked with types: dolfinx.cpp.fem.DirichletBC_float64, ndarray, list, dolfinx.fem.function.FunctionSpace

During handling of the above exception, another exception occurred:

  File "/home/cborau/fenicsx_projects/test13.py", line 79, in <module>
    bc = fem.dirichletbc(u_D, dofs_x0, V0)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: __init__(): incompatible function arguments. The following argument types are supported:
    1. __init__(self, g: ndarray[dtype=float64, writable=False, order='C'], dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C'], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
    2. __init__(self, g: dolfinx.cpp.fem.Constant_float64, dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C'], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
    3. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    4. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: list[ndarray[dtype=int32, writable=False, shape=(*), order='C']], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None

Invoked with types: dolfinx.cpp.fem.DirichletBC_float64, ndarray, list, dolfinx.cpp.fem.FunctionSpace_float64

Thanks in advance

[dokken]

For the boundar conditions, please consider

How to impose spatially varying dirichletbc boundary conditions for one direction in a linear elasticity problem?

This is because the last keyword argument V is only used when you collapse a function into a sub space. Then you send in two arguments in locate_dofs_topological, then uncollapsed and collapsed subspace. This returns you two arrays of degrees of freedom in each of these spaces. Since the collapsed function has no direct relationship to its parent object, one supplied the subspace to ensure that the BCs are imposed on the correct part of the variational form

and/or

Boundary condition to single direction in the subspace of Mixed function space

You are taking the first sub-space of the mixed element, which is V, not V.sub(0) You can check this by printing ME.element.num_sub_elements, which gives you 5, not 6. Thus, you would have to change: To ME0y, _ = ME.sub(0).sub(1).collapse() bottom_value2 = fem.Function(ME0y) bottom_value2.x.array[:]= 0 right_corner_vertex = mesh.locate_entities_boundary(msh,0,right_corner) right_dof = fem.locate_dofs_topological ((ME.sub(0).sub(1),ME0y),0,right_corner_vertex) bc_right_corner = fem.dirichl…

as well as
https://jsdokken.com/dolfinx-tutorial/chapter3/component_bc.html
and
https://docs.fenicsproject.org/dolfinx/v0.8.0/python/demos/demo_stokes.html#non-blocked-direct-solver