# NAN occurs when a custom grid is computed

Hi everyone, I currently create a 1D grid using numpy generated data, but when solving a heat transfer problem prompts NAN, the code is as follows:

``````vertices = np.concatenate([np.linspace(0, 0.5, num=100),np.linspace(0.5, 1, num=200)])

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)):

for j in range(len(vertices)-1):

editor.close()

V = FunctionSpace(mesh, 'P', 1)

.............

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)

.............
``````

Thank you very much!

Without a working example it’s difficult to give advice. For example we don’t know what `k(T)`, `mu` and `c` are.

Here is the complete code:

``````vertices = np.concatenate([np.linspace(0, 0.5, num=100),np.linspace(0.5, 1, num=200)])

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)):

for j in range(len(vertices)-1):

editor.close()

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)
``````

Looks like you have a duplicate point at x=0.5 in

``````vertices = np.concatenate([np.linspace(0, 0.5, num=100),np.linspace(0.5, 1, num=200)])
``````

which gives you a degenerate cell with no volume. If I account for this the system is solvable.