Boundaries not getting marked

Hi all,

I am loading a mesh from gmsh and marking the boundaries (mesh can be found here. I wish to mark the first quadrant as ‘2’, the fourth quadrant as ‘10’ and the second and third quadrant as ‘20’. Those boundaries are not getting marked properly (as seen in the paraview from loading pvd file). Please can anyone tell me the reason for not getting marked properly and how to fix it?

``````from dolfin import *
from fenics import *
import numpy as np
import math
import mshr
import dolfin
from mshr import *
import os
parameters["allow_extrapolation"]=True

dir_path = os.path.dirname(os.path.realpath(__file__))

R = 10.00 # Radius of the circular part

mesh = Mesh("mesh.xml")

class centre(SubDomain):
def inside(self,x, on_boundary):
tol = 3e-1
return (near(x[0],0.00,tol) and near(x[1],0.00,tol))

class circle_boundary(SubDomain):
def inside(self,x, on_boundary):
tol = 1.0
r=(x[0]**2+x[1]**2)**0.5
return ( R-r < tol) and on_boundary

def inside(self,x, on_boundary):
tol = 5e-1
return ((R+tol > x[0] > 0-tol and R+tol > x[1] > 0-tol) and on_boundary)

class circle_bound1(SubDomain):
def inside(self,x, on_boundary):
tol = 5e-1
return (( x[0] > 0 and  x[1] < 0) and on_boundary)

class circle_bound2(SubDomain):
def inside(self,x, on_boundary):
tol = 5e-1
return ((0> x[0] > -R and R > x[1] > -R) and on_boundary)

boundaries = MeshFunction("size_t",mesh,mesh.topology().dim()-1)
boundaries.set_all(0)
centre= centre()
circle_boundary = circle_boundary()
circle_bound1 = circle_bound1()
circle_bound2 = circle_bound2()
centre.mark(boundaries,5)
circle_boundary.mark(boundaries,1)
circle_bound1.mark(boundaries,10)
circle_bound2.mark(boundaries,20)

meshfile=File(dir_path+'/mesh.pvd')
meshfile<<mesh
#plot(boundaries, interactive=True)

boundary_file = File(dir_path+'/boundaries1.pvd')
boundary_file <<boundaries

raise SystemExit()
``````

After some cleanup with names, this is what I get:

``````from dolfin import *

R = 10.00  # Radius of the circular part

mesh = Mesh("mesh.xml")

class centre(SubDomain):
def inside(self, x, on_boundary):
tol = 3e-1
return near(x[0], 0.00, tol) and near(x[1], 0.00, tol)

class circle_boundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1.0
r = (x[0] ** 2 + x[1] ** 2) ** 0.5
return (R - r < tol) and on_boundary

def inside(self, x, on_boundary):
tol = 5e-1
return (R + tol > x[0] > 0 - tol and R + tol > x[1] > 0 - tol) and on_boundary

class circle_bound1(SubDomain):
def inside(self, x, on_boundary):
tol = 5e-1
return (x[0] > 0 and x[1] < 0) and on_boundary

boundaries = MeshFunction("size_t", mesh, 0)
boundaries.set_all(0)

bcentre = centre()
bcircle_boundary = circle_boundary()
bcircle_bound1 = circle_bound1()

bcentre.mark(boundaries, 5)
bcircle_boundary.mark(boundaries, 1)
bcircle_bound1.mark(boundaries, 10)

from vtkplotter.dolfin import plot
plot(boundaries)
``````

PS: to install vtkplotter you can use:

``````pip install vtkplotter
``````

This does not answer my question. Yes! That is what I am getting. But I wish to mark the entire quadrant. Please let me know if you have any workaround to mark the entire quadrant rather than some portion of it. (from the condition provided in the snippet, I expect to get marked the entire quadrant.)

Lets say we have a circle with radius = 1 and the quadrants are:

• First quadrant: Top Right ## Should be marked as 2
• Second quadrant: Top Left ## Should be marked as 20
• Third quadrant: Bottom Left ## Should be marked as 20
• Fourth quadrant: Bottom Right ## Should be marked as 10

Here is the GMSH file and how you can tag the quadrants:

``````Point(1) = {0, 0, 0, 1.0};
Point(2) = {1, 0, 0, 1.0};
Point(3) = {0, 1, 0, 1.0};
Point(4) = {-1, 0, 0, 1.0};
Point(5) = {0, -1, 0, 1.0};

Circle(1) = {2, 1, 3};
Circle(2) = {4, 1, 3};
Circle(3) = {4, 1, 5};
Circle(4) = {2, 1, 5};

Line Loop(5) = {2, -1, 4, -3};
Plane Surface(6) = {5};

Physical Surface(1) = {6};
Physical Line(2) = {1};
Physical Line(20) = {2,3};
Physical Line(10) = {4};
``````

The mesh is something like this:

Then save the mesh in msh format (mesh.msh). Now you can convert it to xml format by running this in terminal:

`````` dolfin-convert mesh.msh newmesh.xml
``````

Three files are created:

1. newmesh.xml
2. newmesh_facet_region.xml
3. newmesh_physical_region.xml

Then run this:

``````from dolfin import *
mesh = Mesh("newmesh.xml")
facets = MeshFunction("size_t", mesh, "newmesh_facet_region.xml")
f = File("facets.pvd")
f<< facets
``````

When you open the facets.pvd file, this is what you can see:

1 Like

…indeed your mesh looks a bit weird:
`vtkplotter mesh.xml -c black -bg white -a 0.1`

it’s not on a plane and its cell orientations are not uniform:

In any case you can play with the SubDomain to select the region you wish, e.g.

``````from dolfin import *

R = 10.00  # Radius of the circular part

mesh = Mesh("mesh.xml")

class sector1(SubDomain):
def inside(self, x, on_boundary):
return x[1]>0 and x[0]>0

class rim(SubDomain):
def inside(self, x, on_boundary):
r = sqrt(x[0]**2+x[1]**2)
return 9.99<r

boundaries = MeshFunction("size_t", mesh, 0)
boundaries.set_all(0)

bs1 = sector1()
bs1.mark(boundaries, 5)

brm = rim()
brm.mark(boundaries, 10)

from vtkplotter.dolfin import plot
plot(boundaries)
``````