在Python FEniCS中基于函数图像标记网格边界的方法

完全可行。核心思路是在Python中直接利用现有网格的节点坐标，通过数学条件筛选出落在目标曲线上的边界单元，手动完成面标记，无需回到Gmsh重新处理网格。

实现步骤 #

1. 导入网格并提取边界信息 #

以DOLFINx为例，先导入已有的XDMF网格文件，获取所有边界单元及其节点坐标：

from dolfinx.io import XDMFFile
from mpi4py import MPI
import dolfinx.mesh

# 读取网格文件
mesh = XDMFFile(MPI.COMM_WORLD, "你的网格文件.xdmf", "r").read_mesh()
# 获取所有边界单元（dim=1对应2D网格的边）
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, dim=1, marker=lambda x: True)

2. 定义目标曲线的判断条件 #

函数$f:(x,y)↦(2x,3y)$作用于单位圆$S^{1$（满足$a}2 + b^2 = 1$），得到的曲线是椭圆，其方程为：
$$\left(\frac{x}{2}\right)^2 + \left(\frac{y}{3}\right)^2 = 1$$
考虑数值浮点误差，我们用一个极小阈值来判断节点是否落在这条椭圆上：

import numpy as np

def is_on_target_curve(x, y, tol=1e-6):
    return np.abs((x/2)**2 + (y/3)**2 - 1) < tol

3. 筛选并标记目标边界单元 #

遍历所有边界单元，检查其所有节点是否符合椭圆条件，标记符合要求的单元：

# 初始化标记数组，默认值为0（未标记）
marker = np.zeros(len(boundary_facets), dtype=np.int32)
target_tag = 1  # 给目标椭圆分配的标记ID

# 向量化筛选（比循环更高效）
# 获取所有边界单元的节点坐标
facet_coords = mesh.geometry.x[mesh.geometry.dofmap[boundary_facets]]
x = facet_coords[:, :, 0]
y = facet_coords[:, :, 1]
# 检查每个单元的所有节点是否都在椭圆上
on_curve = np.all(is_on_target_curve(x, y), axis=1)
# 标记符合条件的单元
marker[on_curve] = target_tag

4. 创建面标记并用于PDE求解 #

将标记结果转为网格的面标签，后续即可在边界条件定义中直接使用：

from dolfinx.mesh import meshtags
from dolfinx.fem import DirichletBC, FunctionSpace, Function

# 创建面标记对象
facet_tags = meshtags(mesh, dim=1, entities=boundary_facets, values=marker)

# 示例：用标记定义Dirichlet边界条件
V = FunctionSpace(mesh, ("CG", 1))
u_D = Function(V)
u_D.interpolate(lambda x: 0.0)
bc = DirichletBC(u_D, facet_tags.find(target_tag), V)

关键注意事项 #

• 数值精度：根据网格的离散精度调整tol阈值，避免因浮点误差漏标或误标单元；
• 高阶网格：如果是二次及以上的高阶网格，需确保单元的所有节点（包括内部节点）都满足条件，避免标记不完整；
• 性能优化：大网格场景下优先使用向量化操作，替代循环遍历，提升处理速度。

内容的提问来源于stack exchange，提问作者user5998510