# Generation of Cylindrical Mesh

Hi Everyone,

I am fairly new to Fenics and I am currently trying to model the heat equation in a 3D cylinder. I have previously developed a model using a cuboid so I could get a hang of the structure of Fenics, but now I am trying to replace this cuboid with a cylinder. In this problem, I am assuming that the metal cylinder is being illuminated on its upper surface by a guassian laser and I want to see how it heats up over time and if a stable temperature distribution exists.

For the cylindrical mesh, I have simply done

``````R = 4.0     #mm
geometry = mshr.Cylinder(Point(0,0,0),Point(0,0,-8),R,R) #8mm long target

# Create mesh
mesh = generate_mesh(geometry, 32)          # generate a mesh from the given geometry
``````

This compiles within my code and does not flag up errors. However, when my simulation begins to run, I notice that the interation stops after 0 iterations:

``````Time 2.0
*** ===================================================== ***
*** FFC: quadrature representation is deprecated! It will ***
*** likely be removed in 2018.2.0 release. Use uflacs     ***
*** ===================================================== ***
issue_deprecation_warning()
*** ===================================================== ***
*** FFC: quadrature representation is deprecated! It will ***
*** likely be removed in 2018.2.0 release. Use uflacs     ***
*** ===================================================== ***
issue_deprecation_warning()
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-03)
Newton solver finished in 0 iterations and 0 linear solver iterations.
``````

Obviously, this shouldnâ€™t happen and indeed does not happen when I use a cuboid. Does anyone have any idea why this happens and how it can be rectified?

Thanks in advance, my full code is below:

``````from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
from ufl import as_tensor
from ufl import Index
import math
import mshr

parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
set_log_level(20)

# SOLVE THE HEAT TRANSFER PROBLEM FOR A 3D TARGET

#--Dimensions in mm----
x0 = 0
y0 = 0
z0 = 0
x1 = 4
y1 = 4
z1 = 8
#----------------
t_start = 0.0
t_end = 100
nstep = 100
dtime = (t_end - t_start)/ nstep
#-----------------
#Configure for dimensions of system

pi = 3.141592653589793
T_am = 0.0 #ambient vacuum temperature
T_a4 = T_am**4
epsilon = 0.2  # material emissivity
sigma = 5.67E-14 # W/(m2.K4)
es = epsilon*sigma
Pmax = 20 #peak power in W
w = 0.5
area = pi*w*w
depth = 5e-4 #0.5um pentration depth
volume = area * depth #define a disk through which the laser is absorbed
Pden = Pmax / volume

Laser = Expression('2*A*Pden*exp(-2*pow((x[0] - 0), 2)/(w*w)-2*pow((x[1]-0), 2)/(w*w)-2*pow((x[2]-z1), 2)/(w2*w2) )',degree=3, A=epsilon, Pden=Pden,w= 0.5, w2= depth,z1=z1) #power (w2 localises the z-coordinates)
#A =! emissivity for system in thermal equillibrium
#----------THERMAL--PROPERTIES--[Element: Y]--------
kappa = 0.0172         #conductivity [W/mm K]
c = 0.3         #Specific Heat Capacity [J/gK]
rho = 4.47e-3           #Density [g/mm^3]
const = kappa /(c * rho)
tau_T = 0.0 #Temperature Gradient Lag
tau_q = 0.0 #Heat Flux Lag
#---------------------------------------
R = 4.0     #mm
geometry = mshr.Cylinder(Point(0,0,0),Point(0,0,-8),R,R) #8mm long target

# Create mesh
mesh = generate_mesh(geometry, 32)          # generate a mesh from the given geometry
#mesh = BoxMesh(Point(x0,y0,z0), Point(x1,y1,z1), 10,10,20)    #create the mesh over which we shall work
Space = FunctionSpace((mesh), 'P', 1)      #define finite element function space, defined via basis functions
VectorSpace = VectorFunctionSpace(mesh, 'P', 1)
cells = MeshFunction('size_t',mesh,mesh.topology().dim()) #codimension of 0
facets = MeshFunction('size_t',mesh,mesh.topology().dim()-1) #codimension of 1
da = Measure('ds', domain=mesh, subdomain_data = facets)  #area element
dv = Measure('dx', domain=mesh, subdomain_data = cells)   #volume element

initial_T = Expression("Tini", Tini=T_am, degree=3) # extrapolate an expression for the temperature before heating
T0 = interpolate(initial_T, Space)
T = Function(Space)         # Temperature
V = TestFunction(Space)     # Test Function used for FEA
dT = TrialFunction(Space)   # Temperature Derivative
q0 = Function(VectorSpace)
i = Index()
G = as_tensor(T.dx(i), (i))  #gradient of T
G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step

q = as_tensor(dtime/(dtime + tau_q) * (tau_q/dtime*q0[i] - kappa*(1+tau_T/dtime)*G[i] + kappa*tau_T/dtime*G0[i]),(i)) #heat
F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i) - rho*Laser*V ) * dv + es*(T**4 - T_a4)*V*da   #final form to solve
Gain = derivative(F, T, dT)    # Gain will be usedas the Jacobian required to determine the evolution of a dynamic system

file_T = File('target3D/solution.pvd')
for t in np.arange(t_start,t_end,dtime):
print( "Time", t)
solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "quadrature","quadrature_degree": 2} )
file_T << (T,t)
q0 = project(q, VectorSpace)
T0.assign(T)      #change so that T0 is equal to T for the next time step
``````

By changing the ambient vacuum temperature from 0 to a non-zero value, you can observe that it takes one iteration to solve the problem (and the solution becomes the ambient vacuum temperature). This means that the initial guess as a zero value is the solution to the problem. This indicates that you have changed something in your variational formulation when moving from a cube to a cylinder.