How do I save a custom grid？

General
1

Hi everyone, my problem should be a grid problem, when I create a grid and calculate it in the following way, the program works fine, but if I create a grid first and save it in the following way, and then import it in the code, the program always does not converge:

with XDMFFile("mesh/mesh.xdmf") as outfile: 
    outfile.write(mesh)

mesh = Mesh()
with XDMFFile("mesh/mesh.xdmf") as infile:
    infile.read(mesh)

######################################
# 创建网格
vertices_start = np.concatenate([np.linspace(0, 0.5, num=100),np.linspace(0.5, 1, num=200)])

vertices = []
for i in vertices_start:
    if i not in vertices:
        vertices.append(i)
mesh = Mesh()
editor = MeshEditor()

editor.open(mesh, 'interval', 1, 1)
editor.init_vertices(len(vertices))  
editor.init_cells(len(vertices)-1)    

for i in range(len(vertices)):
    editor.add_vertex(i, np.array([vertices[i]]))
    
for j in range(len(vertices)-1):
    editor.add_cell(j, np.array([j, j+1]))

editor.close()
######################################
# 导入网格
#mesh = Mesh()    # 不收敛
#with XDMFFile("mesh/mesh.xdmf") as infile:
    #infile.read(mesh)
######################################
V = FunctionSpace(mesh, 'P', 1)

def k(T):  
    return 174.9274 - 0.1067 * T + 5.0067e-5 * T**2 - 7.8349e-9 * T**3

mu = Constant(19300)  
c = Constant(137)    
q = Expression('1e7', degree=1)   
T_ini = Expression('298.15', degree=1)    
f = Constant(0.0)  

t_total = 10   
num_steps = 500    
dt = t_total / num_steps  

class BoundaryX0(SubDomain):   
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 0.00) < DOLFIN_EPS

class BoundaryX1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 1) < DOLFIN_EPS

bx0 = BoundaryX0() 
bx1 = BoundaryX1()

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1) 

bx0.mark(boundary_markers, 0) 
bx1.mark(boundary_markers, 1)

boundary_conditions = {0: {'Neumann': q}, 1: {'Dirichlet': T_ini}}

bcs = []  
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
        bcs.append(bc)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
        
T_n = interpolate(T_ini, V)

T = Function(V)
v = TestFunction(V)

integrals_N = []  
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        q = boundary_conditions[i]['Neumann']
        integrals_N.append(q*v*ds(i))

F = mu*c*T*v*dx + dt*k(T)*dot(grad(T), grad(v))*dx - mu*c*T_n*v*dx - dt*f*v*dx - dt*sum(integrals_N)

J = derivative(F, T)
problem = NonlinearVariationalProblem(F, T, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['absolute_tolerance'] = 1e-6
prm['newton_solver']['relative_tolerance'] = 1e-7
prm['newton_solver']['maximum_iterations'] = 100
prm['newton_solver']['relaxation_parameter'] = 1.0

t = 0
plt.figure(figsize=(12,8))
for n in range(num_steps):
    
    t += dt    
  
    solver.solve()
      
    T_n.assign(T)

Thank you very much！

2

I cannot reproduce the problem (the code converges for me regardless of reading the mesh from file or not).
What version of DOLFIN are you running?
How did you install dolfin?
What is the actual output of running the code?

3

Sorry, I answered late, the previous problem has been solved, but there is a new problem, The absolute error is 0, the relative error is -nan，A similar problem occurs with heat flux, which I think has to do with ds in 1D,The code is as follows:

from fenics import *
import numpy as np

vertices_start = np.concatenate([np.linspace(0, 1e-6, num=200),np.linspace(1e-6, 50e-6, num=500)])

vertices = []
for i in vertices_start:
    if i not in vertices:
        vertices.append(i)
mesh = Mesh()
editor = MeshEditor()

editor.open(mesh, 'interval', 1, 1)
editor.init_vertices(len(vertices))  
editor.init_cells(len(vertices)-1)    

for i in range(len(vertices)):
    editor.add_vertex(i, np.array([vertices[i]]))
    
for j in range(len(vertices)-1):
    editor.add_cell(j, np.array([j, j+1]))

editor.close()

L = FunctionSpace(mesh, 'CG', 1)

def k(T):    
    return 174.9274 - 0.1067 * T + 5.0067e-5 * T**2 - 7.8349e-9 * T**3

mu = Constant(19300)  
c = Constant(137) 
T_ini = Expression('298.15', degree=1)  
T_bottom = Expression('298.15', degree=1)
T_top = Expression('500', degree=1)
f = Constant(0.0)  

t_total = 10    
num_steps = 100   
dt = Constant(t_total / num_steps)    # 时间步长(s)

class BoundaryX0(SubDomain):    
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 0.0) < DOLFIN_EPS

class BoundaryX1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 50e-6) < DOLFIN_EPS

bx0 = BoundaryX0()   
bx1 = BoundaryX1()

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)    

bx0.mark(boundary_markers, 0)    
bx1.mark(boundary_markers, 1)

boundary_conditions = {0: {'Dirichlet': T_top}, 1: {'Dirichlet': T_bottom}}

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

bcs = []   
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(L, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
        bcs.append(bc)
        
T_n = interpolate(T_ini, L)

T = Function(L)
v = TestFunction(L)

F = mu(T)*c(T)*T*v*dx + dt*k(T)*dot(grad(T), grad(v))*dx - mu(T)*c(T)*T_n*v*dx - dt*f*v*dx

xdmffile_T = XDMFFile('heat_transfer/temperature.xdmf')

xdmffile_T.parameters["flush_output"] = True

J = derivative(F, T)
problem = NonlinearVariationalProblem(F, T, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['absolute_tolerance'] = 1e-2
prm['newton_solver']['relative_tolerance'] = 1e-4
prm['newton_solver']['maximum_iterations'] = 100
prm['newton_solver']['relaxation_parameter'] = 1.0

t = 0
while t < t_total:
    
    t += float(dt)
    solver.solve()
    xdmffile_T.write(T, t)
    T_n.assign(T)
4

Also, when I reduce the size of the model (very small about 1e-6m), the relative and absolute errors will become larger, is this reasonable? Thank you!