Bondary conditions are not applied when applying a mshr-mesh

[Former user]

My boundary conditions are not applied when applying a mesh made with mshr, code as follows

import mshr as m

N = 250

r1 = m.Rectangle(Point(0.0, 0.0), Point(L, H))
r2 = m.Rectangle(Point(4.5*D, 14*D/2. -0.5*D), Point(5.5*D, 14*D/2. + 0.5*D))
domain = r1 - r2

mesh = m.generate_mesh(domain, N, "cgal")

When I change it to use FEniCS' RectangleMesh instead just making r1, the BCs are applied just fine.

Also I get the following error when running in parallel and using the mshr-mesh:

self.u_tilde.vector().axpy(1.5,self.u2.vector())
RuntimeError:     self.u_tilde.vector().axpy(1.5,self.u2.vector())

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***  [e-mail]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to access vector of degrees of freedom.
*** Reason:  Cannot access a non-const vector from a subfunction.
*** Where:   This error was encountered inside Function.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Again, using FEniCS' RectangleMesh fixes this problem too.

Am I doing something wrong? In the dumped VTK-files the mshr-mesh looks just fine!

[Benjamin Dam Kehlet]

I suspect the first case to vertices being moved a bit by the mesh generetor (see https://bitbucket.org/benjamik/mshr/issue/14/csgcgalmeshgenerator2d-keep-location-of ). Do you have a small complete example so I can verify that. If this is the case then it will be fixed very soon.

I have no idea about the latter case, but I will take a look. Maybe the mesh is not being correctly distributed to the processes?

[Former user]

The second problem: you are correct, the mesh is not distributed at all. I tried plotting it, one process plots the whole mesh.

As for the first problem run this small script and you get the message that "No facets matching domain...":

from dolfin import *
from mshr import Rectangle, generate_mesh

D = 0.01
L = 20.5*D
H = 14*D

N = 250

r1 = Rectangle(Point(0., 0.), Point(L, H))
r2 = Rectangle(Point(4.5*D, 14*D/2. -0.5*D), Point(5.5*D, 14*D/2. + 0.5*D))
domain = r1 - r2
mesh = Mesh(generate_mesh(domain, N))

V = VectorFunctionSpace(mesh, "CG", 2)
u1 = Function(V)

bcu_dict = {"on_boundary && x[0] > DOLFIN_EPS && std::abs(0.205-x[0]) > \
        DOLFIN_EPS": (0,0), \
        "on_boundary && x[0] < DOLFIN_EPS": (1,0)}

bcu = [DirichletBC(V, bcu_dict[key], key) for key in bcu_dict]

[bc.apply(u1.vector()) for bc in bcu]

plot(u1)
interactive()

[Benjamin Dam Kehlet]

Ok, thanks. I'll take a look. As a temporary workaround (for the second problem) you can probably save the mesh to file and then load it in parallel.

[Benjamin Dam Kehlet]

The first case should be resolved now. I disabled the snap rounding if there are no subdomains. This should be safe. I also decreased the pixel size of the snap rounding, so it should be within the tolerance of dolfin::near(). I'll take a look at the issue in parallel soon. (Note that the return value of generate_mesh() is a dolfin mesh so you have to call the constructor. mesh = generate_mesh(...) should suffice.)

[Benjamin Dam Kehlet]

Fix distribution of generated mesh in parallel. Resolves issue #15.

→ <<cset 98ebc9c52e73>>

[Claire Lestringant]

Hello,

I am facing the same kind of issue with mshr generated mesh. I use the following python script (file meshgen.py) to generate mesh :

import dolfin
from mshr import *
def set_mesh(ref,h,w):
    p1 = dolfin.Point(0.,0.)
    p2 = dolfin.Point(w,h)
    p3 = dolfin.Point(w,1/2.0*h)
    domain =   Rectangle(p1,p2)
    domain.set_subdomain(1, Rectangle(p1, p3))
    mesh2d = generate_mesh(domain, ref)
    return (mesh2d)

And then I apply following boundary conditions :

from meshgen import set_mesh
from dolfin import *

h = 1.
w = 1.2
mesh= set_mesh(10,h,w)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

class Corner(SubDomain):
  def inside(self, x, on_boundary):
    return x[0]==0.0 and x[1]==0.0

class Border(SubDomain):
  def inside(self, x, on_boundary):
    return x[1]==0.0

border = Border()
corner = Corner()
bc1 = DirichletBC(V, Constant((0.,0.)),corner,method="pointwise")
bc2 = DirichletBC(V.sub(1), 0., border)
bcs = [bc1, bc2]

And I get the following error code while trying to solve a variational problem over this mesh, including boundary conditions defined above :

*** Warning: Found no facets matching domain for boundary condition.

[yanyt12]

hello, lestringant: I have the same problem to you, and how did you solve it?