# Basic example does nothing

**URL:** <https://fenicsproject.discourse.group/t/basic-example-does-nothing/4971>\
**Category:** variational formulation\
**Created:** [January 26, 2021, 8:04am UTC](https://fenicsproject.discourse.group/t/basic-example-does-nothing/4971 "2021-01-26T08:04:04Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![mhechthz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mhechthz/32/2100_2.png) [@mhechthz](https://fenicsproject.discourse.group/u/mhechthz)\
**Post date:** [January 26, 2021, 8:04am UTC](https://fenicsproject.discourse.group/t/basic-example-does-nothing/4971/1 "2021-01-26T08:04:04Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mhechthz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mhechthz/32/2100_2.png) [@mhechthz](https://fenicsproject.discourse.group/u/mhechthz)\
**Post date:** [January 26, 2021, 9:28am UTC](https://fenicsproject.discourse.group/t/basic-example-does-nothing/4971/2 "2021-01-26T09:28:42Z")

</div>

Ok, I searched again for some hours and found, that the book example missed the following line:

`boundary_markers.set_all(9999)`

Therefore the inner of the domain was also assigned to boundary 0. After correction all works well.