Hello,
I apologize for not providing a MWE, but I hope my question is clear from the title.
Is there a way to instruct the solver to ignore one dimension in a scenario where the imported mesh is only one cell thick? Specifically, so I don’t need to set up boundary conditions for the front and back sides of the mesh, as I’m not concerned with the z-dimension.
I’m looking for something similar to the “empty” patch type in OpenFOAM.
Thank you!
Hi.
You could use the create_submesh function over a 3D mesh to solve a 2D problem on a mesh immersed in a 3D space. See for example https://github.com/FEniCS/dolfinx/blob/06e90b8506d967f1836ffac7f289f49cb0fd17ec/python/test/unit/fem/test_assemble_submesh.py#L102. For example, the following solves a Poisson equation on the xz-plane extracted from a facet of the unit cube domain.
# +
import numpy as np
import ufl
from dolfinx import fem, mesh, io
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import ds, dx, grad, inner
msh3D = mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
edim = msh3D.topology.dim - 1
entities = mesh.locate_entities_boundary(msh3D, edim, lambda x: np.isclose(x[1], 0.0))
smsh_tuple = mesh.create_submesh(msh3D, edim, entities)
# smsh_tuple consist of Mesh(submsh, submsh_domain), entity_map, vertex_map, geom_map
msh = smsh_tuple[0] # mesh on the x-z plane
V = fem.functionspace(msh, ("Lagrange", 1))
# -
facets = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
np.isclose(x[2], 0.0)))
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
# +
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
g = ufl.sin(5 * x[0])
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx + inner(g, v) * ds
# -
problem = LinearProblem(a, L, bcs=[bc], petsc_options={
"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
with io.XDMFFile(msh.comm, "out_poisson/poisson2D_on_3Dmesh.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(uh)
Hope this helps.
Cheers.
You can also slightly modify a varitional form (if containing vector components) such that it extends in 3D, i.e.
# Solve Stokes eq in XY-plane of a 3D mesh
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT
from mpi4py import MPI
import dolfinx.fem.petsc
import dolfinx
import basix.ufl
import ufl
import numpy as np
nu_val = 1
# Create mesh and mixed function space
mesh = dolfinx.mesh.create_unit_cube(
MPI.COMM_WORLD, 12, 10, 1)
el_u = basix.ufl.element(
"Lagrange", mesh.basix_cell(), 2, shape=(2,)
)
el_p = basix.ufl.element(
"Lagrange",
mesh.basix_cell(),
1,
)
el_w = basix.ufl.mixed_element([el_u, el_p])
W = dolfinx.fem.functionspace(mesh, el_w)
v, q = ufl.TestFunctions(W)
# Define solution vector and/or trial function
w = dolfinx.fem.Function(W)
u, p = ufl.TrialFunctions(W)
# Define integration measures
dx = ufl.Measure("dx", domain=mesh)
ds = ufl.Measure("ds", domain=mesh)
uhat = ufl.as_vector([u[0], u[1], 0])
vhat = ufl.as_vector([v[0], v[1], 0])
# Add artificial forcing in y direction to illustrate the effect of solving in 2D
dst = dolfinx.default_scalar_type
x = ufl.SpatialCoordinate(mesh)
nu = dolfinx.fem.Constant(mesh, dst(nu_val))
f = ufl.as_vector((0,10*x[1],0))
g = dolfinx.fem.Constant(mesh, dst(0))
# Convenience functions for the variational form
def a(u, v):
return nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
def b(p, v):
return -p * ufl.div(v) * dx
def L(v, q):
return ufl.inner(f, v) * dx - ufl.inner(g, q) * dx
# Define stabilized variational form
F = a(uhat, vhat) - b(p, vhat) - b(q, uhat) - L(vhat, q)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[0], 0))
V, _ = W.sub(0).collapse()
left_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0),V), mesh.topology.dim-1, left_facets)
u_bc = dolfinx.fem.Function(V)
u_bc.interpolate(lambda x: (x[1]*(1-x[1]), np.zeros_like(x[0])))
bc_u = dolfinx.fem.dirichletbc(u_bc, left_dofs, W.sub(0))
walls = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[1], 0)|np.isclose(x[1], 1))
wall_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0),V), mesh.topology.dim-1, walls)
u_wall = dolfinx.fem.Function(V)
u_wall.interpolate(lambda x: (np.zeros_like(x[0]), np.zeros_like(x[0])))
bc_wall = dolfinx.fem.dirichletbc(u_wall, wall_dofs, W.sub(0))
bcs = [bc_u, bc_wall]
problem = dolfinx.fem.petsc.LinearProblem(ufl.lhs(F), ufl.rhs(F), bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
wh = problem.solve()
uh = wh.sub(0).collapse()
converged = problem.solver.getConvergedReason()
print(converged)
assert converged > 0
with dolfinx.io.VTXWriter(mesh.comm, "stokes.bp", [uh], engine="BP5") as writer:
writer.write(0.0)
where I have defined
uhat = ufl.as_vector([u[0], u[1], 0])
vhat = ufl.as_vector([v[0], v[1], 0])
which on a standard 2D grid would be
uhat = u
vhat = v
yielding the following with some artificial forcing in the y direction:
You could of course follow @Jesus_Vellojin proposal if you do not need the 3D geometry at a later stage with the extended field.
