Thank you for your response. Below, I’ve provided a part of my .msh file along with the modified code snippet that includes more detailed boundary condition definitions and the use of boundary markers.
Part of the .msh File:
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
9
2 92 “inlet1”
2 93 “inlet2”
2 94 “outlet1”
2 95 “outlet2”
2 96 “interface1”
2 97 “interface2”
2 98 “walls”
3 90 “channels”
3 91 “chamber”
$EndPhysicalNames
$Nodes
1340
1 5 2.5 4.19999999925
2 4.5 3 4.19999999925
3 5 -2.5 4.19999999925
4 4.5 -3 4.19999999925
5 -4.5 -3 4.19999999925
6 -5 -2.5 4.2
7 -5 2.5 4.19999999925
8 -4.5 3 4.19999999925
9 -0.3 -7.347638122934e-17 4.2
10 4.49999999925 3 0.2
…
code:
import dolfinx
import meshio
from dolfinx import fem
from dolfinx.io import gmshio
from mpi4py import MPI
from dolfinx.io.gmshio import model_to_mesh
import ufl
from dolfinx.cpp.mesh import cell_entity_type
from dolfinx.cpp.mesh import to_type
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological,
locate_dofs_geometrical)
from dolfinx.io import XDMFFile
from petsc4py import PETSc
from ufl import (FacetNormal, FiniteElement, Identity, TestFunction, TrialFunction, VectorElement,
div, dot, ds, dx, inner, lhs, grad, rhs, sym)
from dolfinx.fem.petsc import LinearProblem
from dolfinx import default_scalar_type
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
import numpy as np
proc = MPI.COMM_WORLD.rank
path = '/path/'
filename = 'mesh'
mesh, cell_markers, facet_markers = gmshio.read_from_msh(path + filename + '.msh', MPI.COMM_WORLD, gdim=3)
# parameters
mu = 0.001003e-6
rho = 1.0e-9
nu = mu/rho
flow_rate = 1.0e-8
# Define boundary conditions using facet markers
boundary_numbers = {'inlet1': 92, 'inlet2': 93, 'walls': 98, 'outlet1': 94, 'outlet2': 95}
dx = ufl.Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_markers)
inlet1_area = fem.assemble_scalar(fem.form(1 * ds(boundary_numbers['inlet1'])))
inl_velocity = flow_rate/(rho*inlet1_area)
tdim = mesh.topology.dim
fdim = tdim - 1
# Define function spaces for velocity and pressure
P2 = ufl.VectorElement('CG', mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement('CG', mesh.ufl_cell(), 1)
element = ufl.MixedElement([P2, P1])
W = fem.FunctionSpace(mesh, element)
V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()
# Define the velocity part of the boundary condition
velocity_bc_value = fem.Constant(mesh, (0.0, 0.0, -1.0))
noslip_bc_value = fem.Constant(mesh, (0.0, 0.0, 0.0))
pressure_bc_value = fem.Constant(mesh, 0.0)
# Locate degrees of freedom associated with the specified boundary
bcs = []
boundary_indices = fem.locate_dofs_topological(V, fdim, facet_markers.indices[facet_markers.values == boundary_numbers['inlet1']])
bcs.append(fem.dirichletbc(velocity_bc_value, boundary_indices, V))
boundary_indices = fem.locate_dofs_topological(V, fdim, facet_markers.indices[facet_markers.values == boundary_numbers['inlet2']])
bcs.append(fem.dirichletbc(velocity_bc_value, boundary_indices, V))
boundary_indices = fem.locate_dofs_topological(V, fdim, facet_markers.indices[facet_markers.values == boundary_numbers['walls']])
bcs.append(fem.dirichletbc(noslip_bc_value, boundary_indices, V))
boundary_indices = fem.locate_dofs_topological(Q, fdim, facet_markers.indices[facet_markers.values == boundary_numbers['outlet1']])
bcs.append(fem.dirichletbc(pressure_bc_value, boundary_indices, Q))
boundary_indices = fem.locate_dofs_topological(Q, fdim, facet_markers.indices[facet_markers.values == boundary_numbers['outlet2']])
bcs.append(fem.dirichletbc(pressure_bc_value, boundary_indices, Q))
# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Constant(mesh, PETSc.ScalarType((0, 0, 0)))
a = (
+ nu * inner(grad(u), grad(v))*dx
- 1/rho * inner(p, div(v))*dx
- q * div(u) * dx
)
L = inner(f, v) * dx
# Solve the variational problem
problem = LinearProblem(a, L, bcs)
u_p = problem.solve()
# Split the solution
u, p = u_p.split()
V_out = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1)
Q_out = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W_out = FunctionSpace(mesh, V_out * Q_out)
u_p_out = Function(W_out)
u_out, p_out = u_p_out.split()
u_out.interpolate(u)
p_out.interpolate(p)
with XDMFFile(MPI.COMM_WORLD, path + filename + "_velocity1.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(u_out)
with XDMFFile(MPI.COMM_WORLD, path + filename + "_pressure1.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(p_out)
I suspect that the issue might be related to how I’m defining or applying these boundary conditions, or possibly something nuanced with the mesh itself that I’m overlooking. I’d deeply appreciate any insights or advice on how to ensure that these boundary conditions are correctly applied to achieve accurate simulation results.
Thank you once again for your help!
