Hi,

I’m totally new to FeniCS and FEM (not to PDE) and try to build up a more complex application which I reduced meanwhile to the attached code. Even if so basic I don’t get the example to work, i.e. I expect to get a linear temperature decrease from bound G to bound F.

```
"""
Time dependent Thermal Convection-Diffusion Equation
rho*cp*grad(T,t)+rho*cp*w*grad(T,x)=div(lambda*grad(T,x),x)
T = Tc on Gamma G
T = Tc/2. on Gamma F
Temperatures in degree Celsius
Boundaries:
x: 0 --> Thick
y: 0 --> LengthT (y-Axis inverted)
x=0.0 boundary G x = Thick
0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0
b 4 1
o 4 1 boundary M
u 4 1
n ... ...
d 4 1
a 4 1 y = LengthM
r 4 2
y 4 2 boundary O
... ...
S 4 2
4 2
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 y = LengthT
boundary F
"""
from fenics import *
import matplotlib.pyplot as plt
import sys
#=============================================================
# some problem parameters
rho = 7000 # density kg/m**3
lmbda = 50 # conductivity W/m K
cp = 500 # haéat capacity J/kg K
thd = lmbda/(rho*cp) # thermal diffusivity m**2/s
TG = 1550.0 # Dirichlet Boundary G degree Celsius
TF = TG/2.0 # Dirichlet Boundary F degree Celsius
Tini = 1000.0 # initial temperature
Thick = 0.2 # thickness
CenterX = Thick/2.0
LengthM = 0.6 # length from 0.0 to of first boundary
npThick = 20
LengthT = 3.0 # total length
npLength = 300
#=============================================================
# mesh, functionspace and boundaries
p0 = Point(0.0,0.0)
p1 = Point(Thick,LengthT)
mesh = RectangleMesh(p0,p1,npThick,npLength,'left')
V = FunctionSpace(mesh, 'P', 1)
# ALL OTHER BOUNDARIES ARE IN THIS EXAMPLE REMOVED
class BoundaryG(SubDomain): # Boundary 0 == G
tol = 1E-8
nr = 0
def inside(self, x, on_boundary):
if on_boundary and near(x[1],0.0,self.tol): return True
else: return False
class BoundaryF(SubDomain): # Boundary 1 == F (original Boundary 3 - all others removed)
tol = 1E-8
nr = 1
def inside(self, x, on_boundary):
if on_boundary and \
near(x[1],LengthT,self.tol): return True
else: return False
boundary_markers = MeshFunction('size_t', mesh,dim=mesh.topology().dim()-1)
boundG = BoundaryG()
boundG.mark(boundary_markers,boundG.nr)
boundF = BoundaryF()
boundF.mark(boundary_markers,boundF.nr)
# for Neumann and Robin boundary conditions (later)
ds = Measure('ds',domain=mesh,subdomain_data=boundary_markers)
boundary_file = File("t_modell/boundaries.pvd")
boundary_file << boundary_markers
# check boundary assignments
# for x in mesh.coordinates():
# if bG.inside(x,True): print("%s is on G" % x)
# elif bF.inside(x,True): print("%s is on F" % x)
#=============================================================
# boundary conditions
u = TrialFunction(V)
v = TestFunction(V)
# default bc is Neumann symmetry
bcond = {
boundG.nr:{'Dirichlet':Constant(TG)},
boundF.nr:{'Dirichlet':Constant(TF)}
}
bcs = []
for b in bcond:
bc = DirichletBC(V,bcond[b]['Dirichlet'],boundary_markers,b)
bcs.append(bc)
#=============================================================
# Define variational problem
fDC = Constant(thd) # thermal diffusivity
dt = 1.0 # time step size - here constant
fDt = Constant(dt)
# Define initial value
fIni = Constant(Tini) # Initial Values
u_n = interpolate(fIni,V) # old / initial values
#===========================================================
# PDE
#
# rho*cp*grad(T,t)+rho*cp*w*grad(T,x)=div(lambda*grad(T,x),x)
#
#===========================================================
# weak problem formulation without convection
F = u*v*dx + fDt*fDC*dot(grad(u),grad(v))*dx - u_n*v*dx
a, L = lhs(F), rhs(F)
print("a:",a)
print("L:",L)
#=============================================================
# solve time dependent problem
tFinal = 60.0 # final simulation time
# Create VTK file for saving solution
vtkfile = File('t_modell/solution.pvd')
# Time-stepping
u = Function(V)
t = 0
while t < tFinal:
# Update current time
t += dt
print(t,"s /",tFinal,'s')
# Compute solution
solve(a == L, u, bcs)
# Save to file and plot solution
vtkfile << (u, t)
plot(u)
# Update previous solution
u_n.assign(u)
# Hold plot
import matplotlib.pyplot as plt
plt.show()
```

Could anyone with some experience give me a hint? I think it has something to do with the definition of boundary conditions.

Thank’s in advance.