Thank you for the response. Here is a simplified version of the code that I am attempting to run.
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, io
from dolfinx.mesh import locate_entities, locate_entities_boundary
import ufl
import basix.ufl
from petsc4py import PETSc
import os
# Create mesh and tag subdomains
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array([0.0, 0.0]), np.array([1.0, 1.0])],
[32, 32],
cell_type=mesh.CellType.triangle)
# Subdomain tagging: 1 = solid, 2 = fluid
def solid_cells(x):
return x[0] <= 0.5
def fluid_cells(x):
return x[0] > 0.5
solid_cells_indices = locate_entities(domain, domain.topology.dim, solid_cells)
fluid_cells_indices = locate_entities(domain, domain.topology.dim, fluid_cells)
cell_tags = np.zeros(domain.topology.index_map(domain.topology.dim).size_local, dtype=np.int32)
cell_tags[solid_cells_indices] = 1
cell_tags[fluid_cells_indices] = 2
domain.topology.create_connectivity(domain.topology.dim, 0)
cell_tag = mesh.meshtags(domain, domain.topology.dim, np.arange(len(cell_tags)), cell_tags)
dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_tag)
# Mark interface facets for Nitsche terms
fdim = domain.topology.dim - 1
def interface(x):
return np.isclose(x[0], 0.5)
interface_facets = locate_entities_boundary(domain, fdim, interface)
interface_values = np.full(len(interface_facets), 3, dtype=np.int32)
facet_tag = mesh.meshtags(domain, fdim, interface_facets, interface_values)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
# Define finite elements
degree = 2
cell_type = domain.basix_cell()
lagrange_s = basix.ufl.element("Lagrange", cell_type, degree, shape=(domain.topology.dim,))
lagrange_f = basix.ufl.element("Lagrange", cell_type, degree, shape=(domain.topology.dim,))
# - For fluid pressure, a scalar element:
lagrange_p = basix.ufl.element("Lagrange", cell_type, degree - 1)
# Create a mixed element: ordering: [solid displacement, fluid velocity, fluid pressure]
mixed_element = basix.ufl.mixed_element([lagrange_s, lagrange_f, lagrange_p])
# Create function spaces
V_s = fem.functionspace(domain, lagrange_s) # Solid displacement space
V_f = fem.functionspace(domain, lagrange_f) # Fluid velocity space
Q_f = fem.functionspace(domain, lagrange_p) # Fluid pressure space
W = fem.functionspace(domain, mixed_element) # Mixed space
# Define trial and test functions, and split components
u = ufl.TrialFunction(W)
v = ufl.TestFunction(W)
u_s, u_f, p_f = ufl.split(u)
v_s, v_f, q_f = ufl.split(v)
# Build explicit jump terms so that components are formed as standard UFL vectors.
dim = domain.topology.dim
jump_u = ufl.as_vector([u_s[i] - u_f[i] for i in range(dim)])
jump_v = ufl.as_vector([v_s[i] - v_f[i] for i in range(dim)])
# Material parameters and constitutive functions
mu_s = 1.0
lmbda_s = 1.0
mu_f = 1.0
gamma = 10 * degree**2 # Nitsche penalty
def eps(u):
return ufl.sym(ufl.grad(u))
def sigma_s(u):
return lmbda_s * ufl.tr(eps(u)) * ufl.Identity(dim) + 2 * mu_s * eps(u)
def sigma_f(u, p):
return 2 * mu_f * eps(u) - p * ufl.Identity(dim)
n = ufl.FacetNormal(domain)
# Define variational forms
# Solid weak form (on subdomain 1)
a_solid = ufl.inner(sigma_s(u_s), eps(v_s)) * dx(1)
# Fluid weak form (on subdomain 2)
a_fluid = (ufl.real(ufl.inner(sigma_f(u_f, p_f), eps(v_f))) * dx(2)
- ufl.div(v_f) * p_f * dx(2)
- q_f * ufl.div(u_f) * dx(2))
# Nitsche interface terms (on ds(3))
nitsche = (
- ufl.inner(sigma_s(u_s)*n, jump_v) * ds(3)
- ufl.inner(sigma_f(u_f, p_f)*n, jump_v) * ds(3)
- ufl.inner(jump_u, sigma_s(v_s)*n) * ds(3)
- ufl.inner(jump_u, sigma_f(v_f, 0)*n) * ds(3)
+ gamma * ufl.inner(jump_u, jump_v) * ds(3)
)
a = a_solid + a_fluid + nitsche
# Right-hand side: zero source terms
zero_vector = fem.Constant(domain, (0.0, 0.0))
L = ufl.inner(zero_vector, v_s) * dx + ufl.inner(zero_vector, v_f) * dx
# Apply boundary conditions
def left(x):
return np.isclose(x[0], 0.0)
def right(x):
return np.isclose(x[0], 1.0)
def walls(x):
return np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)
left_facets = locate_entities_boundary(domain, fdim, left)
right_facets = locate_entities_boundary(domain, fdim, right)
wall_facets = locate_entities_boundary(domain, fdim, walls)
# Map BCs onto subspaces of the mixed space
left_dofs = fem.locate_dofs_topological((W.sub(0), V_s), fdim, left_facets)
wall_dofs = fem.locate_dofs_topological((W.sub(1), V_f), fdim, wall_facets)
right_p_dofs = fem.locate_dofs_topological((W.sub(2), Q_f), fdim, right_facets)
u_solid_bc = fem.Function(V_s)
u_solid_bc.x.array[:] = 0.0
u_fluid_bc = fem.Function(V_f)
u_fluid_bc.x.array[:] = 0.0
p_out_bc = fem.Function(Q_f)
p_out_bc.x.array[:] = 0.0
bc_solid_fixed = fem.dirichletbc(u_solid_bc, left_dofs, W.sub(0))
bc_noslip = fem.dirichletbc(u_fluid_bc, wall_dofs, W.sub(1))
bc_pressure_outlet = fem.dirichletbc(p_out_bc, right_p_dofs, W.sub(2))
bcs = [bc_solid_fixed, bc_noslip, bc_pressure_outlet]
# Assemble and solve
problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "minres", "pc_type": "lu"})
U = problem.solve()