Boundaries not getting marked

[Peter_James]

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?

Thanks in advance.

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
      
class circle_load(SubDomain):
    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_load = circle_load()
circle_bound1 = circle_bound1()
circle_bound2 = circle_bound2()
#circle_conc = circle_boundary - circle_load
centre.mark(boundaries,5)
circle_boundary.mark(boundaries,1)
circle_load.mark(boundaries,2)
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()

[marcomusy]

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


class circle_load(SubDomain):
    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_load = circle_load()
bcircle_bound1 = circle_bound1()

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


from vtkplotter.dolfin import plot
plot(boundaries)

PS: to install vtkplotter you can use:

pip install vtkplotter

[Peter_James]

@marcomusy

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.)

[Leo]

I would recommend to mark your boundaries in your GMSH file as Physical Line before converting your mesh. Please consider this:
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:

[marcomusy]

…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)