Problems using CSG 2D meshes with subdomains

Question

Hi!

I'm new to fenics and playing with demos and tutorials I encountered problems using CSG constructed meshes with subdomains.

The problem is: I can construct meshes using CSG as long as I don't define any subdomain. But as soon as I use the set_subdomain method I got the message:

"Unable to set given rows ot identity matrix."
"some diagonal element not preallocated (try assembler option keep_diagonal)"
"This error was encountered inside PETScMatrix.cpp"

I'm using Dolfin version 1.3.0

At the end of this message you will find the first Poisson sample in the Tutorial with just a different mesh.

As you can check, the model without subdomains works fine but you just have to uncomment out the line with "geom.set_subdomain(1, Circle(0., 0., 3.5))" to see the error message appear again.

I alredy tried to define new measures for the different domains using different CellFunction() definitions but without luck.

I would be most grateful if someone can provide a clue on how to proceed!

Regards,

Víctor

from dolfin import *

# Create mesh and define function space
geom = Rectangle(0., 0., 5., 5.) - Circle(0., 0., 2.5)
## geom.set_subdomain(1, Circle(0., 0., 3.5))
mesh = Mesh(geom, 35)
V = FunctionSpace(mesh, 'Lagrange', 1)

plot(mesh)

# Define boundary conditions
u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')

def u0_boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u0, u0_boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution and mesh
plot(u)

# Dump solution to file in VTK format
file = File('poisson.pvd')
file << u

# Hold plot
interactive()

Answer 1

Hi, I'd just like to add a few observations. If you plot your mesh, there's a cluster of points near bottom left arch. If you zoom into the cluster, you'll see hanging nodes. Further with

print V.dofmap().tabulate_all_coordinates(mesh).reshape((-1, 2))

you'll see that several dofs have same physical coordinates and some of them lie outside your intended domain. The problem is also visible if you plot the mesh and press v to get labels of the vertices. Finally if you try other domain

geom = Rectangle(0., 0., 5., 5.) - Rectangle(0., 0., 2, 2)
geom.set_subdomain(1, Rectangle(0, 0, 3, 3))

there's no more clustering or hanging nodes but vertex [0,0] is still part of the mesh. Dof associated with this vertex is a source of the PETSc error. I am not sure if this behaviour is unexpected. The CSG demos use A.set_subdomain(1, B) with B a true subdomain of A in a sense that $B-A$ is an empty set and that's not what you are doing. So to make the example with rectangle work, you could try

geom = Rectangle(0., 0., 5., 5.) - Rectangle(0., 0., 2, 2)
diff = Rectangle(0, 0, 3, 3) - Rectangle(0, 0, 2, 2)
geom.set_subdomain(1, diff)

With the circle, it seems you have to leave some margin to make CGAL work

geom = Rectangle(0., 0., 5., 5.) - Circle(0., 0., 2.5)
diff = CSGIntersection(Circle(0., 0., 3.5) - Circle(0., 0., 2.51),
        Rectangle(0.01, 0.01, 3.5, 3.5))
geom.set_subdomain(1, diff)

Comments

Thanks a lot MiroK!

I was so convinced I was doing something wrong with the subdomain definition or the measure that I didn't think about cheking the mesh!

Now the code is working fine (nevertheless, for real work I will probably check external meshers that offer more flexibility on the mesh definition)