Thank you again for the great hint dokken.

This is a part of my code:

```
W = VectorFunctionSpace(mesh, "Lagrange", 1)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V) # For linear case
test_u = TrialFunction(W)
#u = Function(V) # For nonlinear case
v = TestFunction(V)
uh = Function(V)
g = Constant(0.0)
c = Constant(1.0e0)
x = SpatialCoordinate(mesh)
f = sin(2.0*pi*x[0])*cos(2.0*pi*x[1])
n = FacetNormal(mesh)
# Tangential gradient operator:
def grad_t(u):
return grad(u) - dot(grad(u),n)*n
# Randbedingungen-----------------------------------------------------------
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (-1.0, 1.0)) and between(x[0], (4.0, -4.0)))
class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ((between(x[1], (2.0, 2.0)) and between(x[0], (1.0, 3.0))) and (between(x[1], (0.1, 0.0)) and between(x[0], (2.0, 2.0))) and (between(x[1], (2.0, 0.0)) and between(x[0], (0.1, 0.0))) and (between(x[1], (5.0, -3.0)) and between(x[0], (2.0, 0.0))))
##Facetfunction
#facetfunc = MeshFunction("size_t", mesh, mesh.geometry().dim()-1)
#bnd = Right()
#bnd.mark(facetfunc, 3)
#File("facetfunc.pvd").write(facetfunc)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (5.0, -3.0)) and between(x[0], (4.0, -4.0)))
class Top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (1.0, 3.0)) and between(x[0], (-1.0, 1.0)))
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(4)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
# Define Dirichlet boundary condition at top boundary
bc2 = DirichletBC(V, 0.0, boundaries, 2)
#Take a look at the solution:
ds=Measure('ds', domain=mesh, subdomain_data=boundaries)
area_gamma1 = assemble(1*ds(4))
leitfaehigkeit=8.6e6
a = inner(grad(u), grad(v))*dx
L = g*v*ds(1)+g*v*ds(3)+leitfaehigkeit*(1/area_gamma1)*(-zz[1, 3])*v*ds(4)
solve(a==L,uh, bc2)
p = plot(uh)
plot(grad(uh))
cbar = plt.colorbar(p)
plt.show()
plt.savefig('solution' + str(1) + '.png')
```

zz is a numerical solution from an ODE-System.