1. My geo file is:
/ Gmsh project created on Thu Mar 14 01:43:28 2024
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0, 16e-3, 0, 1.0};
//+
Point(3) = {0, 17e-3, 0, 1.0};
//+
Point(4) = {0, 22e-3, 0, 1.0};
//+
Point(5) = {8e-3, 22e-3, 0, 1.0};
//+
Point(6) = {8e-3, 17e-3, 0, 1.0};
//+
Point(7) = {8e-3, 16e-3, 0, 1.0};
//+
Point(8) = {8e-3, 3e-3, 0, 1.0};
//+
Point(9) = {8e-3, 13e-3, 0, 1.0};
//+
Point(10) = {8e-3, 8e-3, 0, 1.0};
//+
Point(11) = {8e-3, 0, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 5};
//+
Line(5) = {5, 6};
//+
Line(6) = {6, 7};
//+
Line(7) = {1, 11};
//+
Line(8) = {11, 8};
//+
Line(9) = {7, 9};
//+
Circle(10) = {9, 10, 8};
//+
Line(11) = {3, 6};
//+
Line(12) = {2, 7};
//+
Curve Loop(1) = {3, 4, 5, -11};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {11, 6, -12, 2};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {1, 12, 9, 10, -8, -7};
//+
Plane Surface(3) = {3};
//+
Transfinite Curve {4, 11, 12} = 160 Using Progression 1;
//+
Transfinite Curve {3, 5} = 100 Using Progression 1;
//+
Transfinite Curve {1, 10} = 160 Using Progression 1;
//+
Transfinite Curve {9, 8} = 30 Using Progression 1;
//+
Transfinite Curve {7} = 20 Using Progression 1;
//+
Transfinite Curve {2, 6} = 20 Using Progression 1;
//+
Transfinite Surface {1};
//+
Transfinite Surface {2};
//+
Physical Curve("ce1", 13) = {3};
//+
Physical Curve("ce2", 14) = {5};
//+
Physical Curve("top", 15) = {4};
//+
Physical Curve("heat transfer", 16) = {10};
2.I converted the msh file to an xdmf file using meshio：
import meshio
mesh_from_file = meshio.read("mesh_for_gmsh/ping.msh")
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
"name_to_read": [cell_data]})
return out_mesh
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("ping_facet.xdmf", line_mesh)
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("ping.xdmf", triangle_mesh)
- Import fenics for calculation：
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
import time
t = sym.symbols('t')
mesh = Mesh()
with XDMFFile("ping.xdmf") as infile:
infile.read(mesh)
V = FunctionSpace(mesh, 'P', 1)
t_total = 1
num_steps = 10
dt = Constant(t_total / num_steps)
tol = 1E-14
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return 16e-3 + tol >= x[1] >= 0 - tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return 17e-3 + tol >= x[1] >= 16e-3 - tol
class Omega_2(SubDomain):
def inside(self, x, on_boundary):
return 22e-3 + tol >= x[1] >= 17e-3 - tol
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
dx = Measure('dx', domain = mesh, subdomain_data = materials)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_2 = Omega_2()
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)
subdomain_2.mark(materials, 3)
class blood(UserExpression):
def __init__(self, materials, **kwargs):
super().__init__(**kwargs)
self.materials = materials
def eval_cell(self, values, x, cell):
material_index = self.materials[cell.index]
if material_index == 1:
values[0] = 1
elif material_index == 2:
values[0] = 2
else:
values[0] = 3
def value_shape(self):
return ()
k = blood(materials, degree=1)
class BoundaryTop(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1] - 22e-3) < DOLFIN_EPS
class BoundaryTransfer(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] > 3e-3 - DOLFIN_EPS and x[1] > 3e-3 - DOLFIN_EPS and x[1] < 13e-3 + DOLFIN_EPS
b_Top = BoundaryTop()
b_Transfer = BoundaryTransfer()
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
b_Top.mark(boundary_markers, 2)
b_Transfer.mark(boundary_markers, 1)
bc1 = DirichletBC(V, 5, boundary_markers, 1)
bc2 = DirichletBC(V, 11, boundary_markers, 2)
bcs = [bc1, bc2]
ds = Measure('ds', domain = mesh, subdomain_data = boundary_markers)
T_ini = Expression('1', degree=1)
T_n = interpolate(T_ini, V)
T = TrialFunction(V)
v = TestFunction(V)
F = T*v*dx + dt*k*dot(grad(T), grad(v))*dx - T_n*v*dx
a, L = lhs(F), rhs(F)
xdmffile_T = XDMFFile('blood.xdmf')
xdmffile_T.parameters["flush_output"] = True
T = Function(V)
t = 0
progress = Progress('Time-stepping')
while t < t_total:
t += float(dt)
solve(a == L, T, bcs)
xdmffile_T.write(T, t)
T_n.assign(T)
set_log_level(LogLevel.PROGRESS)
4.The obtained blood.xdmf file into the paraview visualization can see that there is a 0 value:
5.Then I used to find the 0 point value in the center of the circle，I know it’s not normal, but what went wrong?：
