from dolfin import *
import mshr
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.tri as tri
def build_space(N_circle, N_bulk, u_in):
"""Prepare data for DGF benchmark. Return function
space, list of boundary conditions and surface measure
on the cylinder."""
# Define domain
center = Point(0.2, 0.2) # Point 是FEniCS库中的类，用于创建几何点
radius = 0.05
L = 2.2
W = 0.41
# 矩形域左下点为(0.0, 0.0)，右上点为(L, W)
geometry = mshr.Rectangle(Point(0.0, 0.0), Point(L, W)) \
- mshr.Circle(center, radius, N_circle)
# Build mesh
mesh = mshr.generate_mesh(geometry, N_bulk)
# Construct facet markers
bndry = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
for f in facets(mesh):
mp = f.midpoint()
if near(mp[0], 0.0): # inflow
bndry[f] = 1
elif near(mp[0], L): # outflow
bndry[f] = 2
elif near(mp[1], 0.0) or near(mp[1], W): # walls
bndry[f] = 3
elif mp.distance(center) <= radius: # cylinder
bndry[f] = 5
# Build function spaces (Taylor-Hood)
P2 = VectorElement("P", mesh.ufl_cell(), 2) # 2阶P型有限元素，用于速度场
P1 = FiniteElement("P", mesh.ufl_cell(), 1) # 1阶P型有限元素，用于压力场
TH = MixedElement([P2, P1]) # 混合元素，包含速度场和压力场
W = FunctionSpace(mesh, TH) # 函数空间，包含速度场和压力场
# Prepare Dirichlet boundary conditions
bc_walls = DirichletBC(W.sub(0), (0, 0), bndry, 3) # 对墙壁的速度边界条件，设为0
bc_cylinder = DirichletBC(W.sub(0), (0, 0), bndry, 5) # 对圆柱的速度边界条件，设为0
bc_in = DirichletBC(W.sub(0), u_in, bndry, 1) # 对 inflow 的速度边界条件，设为 u_in
bcs = [bc_cylinder, bc_walls, bc_in] # 合并所有边界条件
# Prepare surface measure on cylinder
ds_circle = Measure("ds", subdomain_data=bndry, subdomain_id=5) # 圆柱表面的积分测量，用于计算圆柱表面的压力项
return mesh, W, bcs, ds_circle
def solve_navier_stokes(W, nu, bcs):
"""Solve steady Navier-Stokes and return the solution"""
# Define variational forms
v, q = TestFunctions(W)
w = Function(W)
u, p = split(w)
F = nu*inner(grad(u), grad(v))*dx + dot(dot(grad(u), u), v)*dx \
- p*div(v)*dx - q*div(u)*dx
# Solve the problem
solve(F == 0, w, bcs)
return w
u_in = Expression(("4.0*U*x[1]*(0.41 - x[1])/(0.41*0.41)", "0.0"),
degree=2, U=0.3)
nu = Constant(0.001)
# Discretization parameters
N_circle = 128
N_bulk = 256
# Prepare function space, BCs and measure on circle
mesh, W, bcs, ds_circle = build_space(N_circle, N_bulk, u_in)
# Solve Stokes
w = solve_navier_stokes(W, nu, bcs)
center = Point(0.2, 0.2)
radius = 0.05
L = 2.2
W = 0.41
bndry = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
for f in facets(mesh):
mp = f.midpoint()
if near(mp[0], 0.0): # inflow
bndry[f] = 1
elif near(mp[0], L): # outflow
bndry[f] = 2
elif near(mp[1], 0.0) or near(mp[1], W): # walls
bndry[f] = 3
elif mp.distance(center) <= radius: # cylinder
bndry[f] = 5
inflow_coords = []
outflow_coords = []
wall_coords = []
cylinder_coords = []
for f in facets(mesh):
if bndry[f] == 1: # inflow
for vertex in vertices(f):
inflow_coords.append(vertex.point().array())
elif bndry[f] == 2: # outflow
for vertex in vertices(f):
outflow_coords.append(vertex.point().array())
elif bndry[f] == 3: # walls
for vertex in vertices(f):
wall_coords.append(vertex.point().array())
elif bndry[f] == 5: # cylinder
for vertex in vertices(f):
cylinder_coords.append(vertex.point().array())
inflow_coords = np.array(inflow_coords)[:,:2]
outflow_coords = np.array(outflow_coords)[:,:2]
wall_coords = np.array(wall_coords)[:,:2]
cylinder_coords = np.array(cylinder_coords)[:,:2]
Above is the solution for the Navier-Stokes cylinder flow problem using dolfin. I need to reproduce it, and I am using DOLFINx 0.10.0. In the dolfin version, there is a code snippet:
for f in facets(mesh):
if bndry[f] == 1: # inflow
for vertex in vertices(f):
inflow_coords.append(vertex.point().array())
elif bndry[f] == 2: # outflow
for vertex in vertices(f):
outflow_coords.append(vertex.point().array())
elif bndry[f] == 3: # walls
for vertex in vertices(f):
wall_coords.append(vertex.point().array())
elif bndry[f] == 5: # cylinder
for vertex in vertices(f):
cylinder_coords.append(vertex.point().array())
How can I similarly achieve this in DOLFINx, extracting coordinate points and their corresponding values according to different boundary markers? Below is the code I wrote using DOLFINx, but note that I am solving a two-dimensional scalar wave equation instead：
# define the geometry
gmsh.initialize()
L = 1
H = 1
c_x = c_y = 0.2
r = 0.05
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0 # 模型所在的进程
if mesh_comm.rank == model_rank:
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
if mesh_comm.rank == model_rank:
fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()
# define the fluid marker
fluid_marker = 1
volumes = gmsh.model.getEntities(dim=gdim)
assert len(volumes) == 1
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker) # 第一个实体的维度、标签、物理组标记
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid") # 第一个实体的维度、物理组标记、物理组名称
# define the boundary markers
inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
boundaries = gmsh.model.getBoundary(volumes, oriented=False) # False表示不考虑边界的方向
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H / 2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H / 2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(
center_of_mass, [L / 2, 0, 0]
):
walls.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
# generate the mesh
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")
# extract the mesh data
mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
mesh = mesh_data.mesh