Invalid domain id

Question

Hi guys,

I wrote a little simulation on python using MacOSX and it worked well. Now I want to copy the simulation to ubuntu, but it doesnt work:

Invalid domain id <Mesh of topological dimension 2 (triangles) with 12345 vertices and 123456 cells, ordered>
[...]
ufl.log.UFLException: Invalid domain id <Mesh of topological dimension 2 (triangles) with 12345 vertices and 123456 cells, ordered>.

Does somebody understand why I cant just copy the .py-file to ubuntu?
It is a 2D-problem in which I declared the domains like:

mesh = RectangleMesh(0.001, 0, a, b, na, nb)
boundary_parts = MeshFunction('size_t', mesh, 1)
boundary_parts.set_all(0)
tol = 1E-14
class Top(SubDomain):
   def inside(self, x, on_boundary):
      return on_boundary and abs(x[1] - b) < tol
top = Top()
top.mark(boundary_parts, 1)
ds = Measure('ds')[boundary_parts]

The traceback leads to a line where I wrote

a = [...]*dx(mesh) + [...]*ds(1)

Any advice?

Best,

MaxMeier

Comments

Hi, are you running same FEniCS version on both?

---

Yes, I am running the latest FEniCS version on both computers (1.5.0).
That's what made me wonder.
I guess there might be some issues with my definition of boundaries but I can't figure out what is wrong about them in my code.

---

Okay, could you post the complete code (but minimal) code?

---

Sure, thank you! Here is an example which works on my mac but not on ubuntu:

from dolfin import *
import sys, time, pylab, numpy

#-------------------------discretisation (time, space):-----------------------------#
nr=200
nz=200
del_t=0.001
t_ges=2
R=20
depth=3

#--------------------------------general parameters-------------------------------#
T0 = 0
boltz = Constant(5.67E-14)
kon = Constant(1E-5)

#--------------------------------heat parameters----------------------------------#
P = Constant(100)
source = Expression('P*(1+sin(x[0]))', P=P)

#--------------------------------create mesh-------------------------------------#
mesh = RectangleMesh(0.001,0,R,depth,nr,nz)
boundary_parts = MeshFunction('size_t', mesh, 1)
boundary_parts.set_all(0)

#--------------------------------create functionspace----------------------------#
V = FunctionSpace(mesh, 'Lagrange', 1)

#------------------------------define boundaries---------------------------------#
tol = 1E-14
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] - depth) < tol
top=Top()
top.mark(boundary_parts, 1)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1]) < tol
bot=Bottom()
bot.mark(boundary_parts, 2)

class Mid(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 0.001) < tol
mid=Mid()
mid.mark(boundary_parts, 3)

class Edge(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - R) < tol
edge=Edge()
edge.mark(boundary_parts, 4)

ds = Measure('ds')[boundary_parts]

#------------------------------define functions---------------------------------#
du = TrialFunction(V)
u=Function(V)
u_1=Function(V)
v = TestFunction(V)
r = Expression('x[0]')
theta = 0.5 #1:forward euler, 0: backward euler, 0.5: crank-nicholson

#-------------------------------start conditions--------------------------------#
u0 = Constant(T0)
u = interpolate(u0,V)
u_1 = interpolate(u0,V)

#------------------------------define nonlinear problem-------------------------#
class HeatEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)

def rad(u):
    return (u**4)

#------------------------------compiler parameters-----------------------------#
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"


#------------------------------weak formulation--------------------------------#
a = (u*v + del_t*theta*(Dx(u,0)*Dx(v,0)+Dx(u,1)*Dx(v,1) - (1/r)*v*Dx(u,0)))*dx(mesh)+ del_t*(u-u0)*kon*v*theta*ds(1) + del_t*kon*(u-u0)*v*theta*ds(2) + del_t*theta*boltz*(rad(u)-u0**4)*v*ds(1) + del_t*theta*boltz*(rad(u)-u0**4)*v*ds(2)

L = (u_1*v - del_t*(1-theta)*(Dx(u_1,0)*Dx(v,0)+Dx(u_1,1)*Dx(v,1) - (1/r)*v*Dx(u_1,0)))*dx(mesh) + del_t*v*(source - (1-theta)*(kon)*(u_1-u0) - (1-theta)*(boltz)*(rad(u_1)-u0**4))*ds(1) + del_t*(1-theta)*(-(kon)*(u_1 - u0) - (boltz)*(rad(u_1)-u0**4))*v*ds(2)

F = a - L               #=0

J = derivative(F,u,du)  #gateaux derivative

problem = HeatEquation(J,F)

#-----------------------------solver and solver parameters--------------------#
solver = NewtonSolver()
#solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6


#------------------------------solve problem---------------------------------#
t = del_t
while t <= t_ges:
    u_1.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    plot(u, mode='color', window_width=900, window_height=600)
    t+=del_t

interactive()

---

Finally I could fix it on my own; I still don't get why it works now and didn't work before.
All I had to do was to change

a=...*dx(mesh)

to

dx = Measure('dx')[boundary_parts]
a = ...*dx()

Of course, I had to change L equally.