Hello, everyone！Currently, I want to conduct numerical simulation for the curved domain. Since I’m not very familiar with it, I started with the simplest Poisson equation. Below is the code I used to calculate it with fenicsx.
import gmsh
from dolfinx.io import gmshio
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import numpy as np
from dolfinx import fem
import ufl
from dolfinx import default_scalar_type
import numpy as np
def solve_poisson_on_mesh(h, element_order, FE_order):
gmsh.initialize()
membrane = gmsh.model.occ.addDisk(0, 0, 0, 1, 1)
gmsh.model.occ.synchronize()
gdim = 2
gmsh.model.addPhysicalGroup(gdim, [membrane], 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)
gmsh.option.setNumber("Mesh.ElementOrder", element_order)
gmsh.option.setNumber("Mesh.SecondOrderLinear", 0)
gmsh.option.setNumber("Mesh.SecondOrderIncomplete", 0)
gmsh.model.mesh.generate(gdim)
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
print(f"网格几何度数: {domain.geometry.cmap.degree}")
print(f"几何映射类型: {domain.geometry.cmap.__class__.__name__}")
print(f"拓扑维度: {domain.topology.dim}")
print(f"顶点数量: {domain.geometry.x.shape[0]}")
print(f"单元数量: {domain.topology.index_map(domain.topology.dim).size_global}")
boundary_points = domain.geometry.x
radii = np.sqrt(boundary_points[:, 0]**2 + boundary_points[:, 1]**2)
on_boundary = np.isclose(radii, 1.0, atol=0.02)
print(f"总点数: {len(boundary_points)}")
print(f"边界点数量: {np.sum(on_boundary)}")
if np.sum(on_boundary) > 0:
print(f"边界点半径范围: [{radii[on_boundary].min():.6f}, {radii[on_boundary].max():.6f}]")
def exact_solution(x):
"""精确解:w = (1 - x² - y²) * sin(πx) * cos(πy)"""
r2 = x[0]**2 + x[1]**2
return (1 - r2) * np.sin(np.pi*x[0]) * np.cos(np.pi*x[1])
def source_term(x):
"""手动计算的源项:f = -∇²w"""
x_val = x[0]
y_val = x[1]
r2 = x_val**2 + y_val**2
term1 = (4 + 2*np.pi**2 * (1 - r2)) * np.sin(np.pi*x_val) * np.cos(np.pi*y_val)
term2 = 4*np.pi*x_val * np.cos(np.pi*x_val) * np.cos(np.pi*y_val)
term3 = -4*np.pi*y_val * np.sin(np.pi*x_val) * np.sin(np.pi*y_val)
return term1 + term2 + term3
def on_boundary(x):
return np.isclose(np.sqrt(x[0]**2 + x[1]**2), 1)
x = ufl.SpatialCoordinate(domain)
V = fem.functionspace(domain, ("Lagrange", FE_order))
# f = fem.Constant(domain, default_scalar_type(4))
f = fem.Function(V)
f.interpolate(source_term)
boundary_dofs = fem.locate_dofs_geometrical(V, on_boundary)
bc = fem.dirichletbc(default_scalar_type(0), boundary_dofs, V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
V_exact = fem.functionspace(domain, ("Lagrange", FE_order +2))
u_exact_func = fem.Function(V_exact)
u_exact_func.interpolate(exact_solution)
error_L2 = fem.form(ufl.inner(uh - u_exact_func, uh - u_exact_func) * ufl.dx)
L2_error = np.sqrt(fem.assemble_scalar(error_L2))
error_H1 = fem.form(ufl.inner(uh - u_exact_func, uh - u_exact_func) * ufl.dx +
ufl.inner(ufl.grad(uh - u_exact_func), ufl.grad(uh - u_exact_func)) * ufl.dx)
H1_error = np.sqrt(fem.assemble_scalar(error_H1))
gmsh.finalize()
return L2_error, H1_error,h
def main():
mesh_sizes = [0.2, 0.1, 0.05, 0.025, 0.0125, 0.00625]
L2_errors = []
H1_errors = []
h_list = []
print("开始收敛性分析...")
print("=" * 60)
for i, h in enumerate(mesh_sizes):
print(f"\n计算网格尺寸 h = {h}")
L2_error, H1_error, actual_h = solve_poisson_on_mesh(h, 2, 2)
L2_errors.append(L2_error)
H1_errors.append(H1_error)
h_list.append(actual_h)
print(f" 网格尺寸: {h:.3f}")
print(f" L2误差: {L2_error:.6e}")
print(f" H1误差: {H1_error:.6e}")
print("\n" + "=" * 60)
print("收敛阶分析:")
print("=" * 60)
print("网格尺寸 L2误差 H1误差 L2收敛阶 H1收敛阶")
print("-" * 75)
for i in range(len(mesh_sizes)):
if i == 0:
print(f"{h_list[i]:.3f} {L2_errors[i]:.2e} {H1_errors[i]:.2e} - -")
else:
# 计算收敛阶
rate_L2 = np.log(L2_errors[i-1] / L2_errors[i]) / np.log(h_list[i-1] / h_list[i])
rate_H1 = np.log(H1_errors[i-1] / H1_errors[i]) / np.log(h_list[i-1] / h_list[i])
print(f"{h_list[i]:.3f} {L2_errors[i]:.2e} {H1_errors[i]:.2e} {rate_L2:.3f} {rate_H1:.3f}")
if __name__ == "__main__":
print("程序开始执行...")
main()
print("程序执行完成")
The result is：
============================================================
Convergence order analysis:
============================================================
h L2 error H1error L2 convergence order H1 convergence order
---------------------------------------------------------------------------
0.200 2.23e-03 8.57e-02 - -
0.100 2.98e-04 2.32e-02 2.908 1.886
0.050 3.74e-05 5.89e-03 2.992 1.976
0.025 4.70e-06 1.48e-03 2.994 1.994
0.013 5.91e-07 3.72e-04 2.990 1.992
0.006 7.41e-08 9.31e-05 2.997 1.998
What I don’t quite understand is whether the calculation in the curved region differs from that in the straight region, and this difference lies in the import process of GMSH? Are there no differences in the others? Is this result of my operation reasonable? Is the degree of approximation controlled by the line gmsh.option.setNumber("Mesh.ElementOrder", element_order)? Looking forward to your reply!