Different behaviour of a mesh coming from generate_mesh() resp. UnitSquareMesh()?

mesh
1

Hi all,

I’ve regarded the following problem:
I have a domain (unit square) and an included subdomain (a smaller rectangle inside). Now I want to refine only the elements inside the subdomain, whose faces are on the boundary of the subdomain. To this end I created the following attempt

from dolfin import *
from mshr import *
from fenics import *

# Constants related to the geometry
sw = Point ( 0.0, 0.0 )
ne = Point ( 1.0, 1.0 )

domain = Rectangle ( sw, ne ) 
#domain.set_subdomain(1, domain)

smalldomain = Rectangle(Point(0.25,0.25),Point(0.75,0.75))
domain.set_subdomain(2, smalldomain)




#  Mesh the region.
#
#
mesh = generate_mesh ( domain, 5)
#mesh= RectangleMesh(sw,ne, 2, 2,"crossed")



#mesh = UnitSquareMesh(20, 20, "crossed")
#mesh = generate_mesh(domain, 20,"cgal")
    
# Define boundaries     
# Define boundaries 
boundary_markers = MeshFunction('size_t', mesh, 2, mesh.domains())
boundaries = MeshFunction('size_t', mesh, 1, mesh.domains())

# Use the cell domains associated with each facet to set the boundary

for j in range(4):
    mesh.init(1, 2) # Initialise facet to cell connectivity
    markers = CellFunction('bool', mesh, False)
        
for f in facets(mesh):
    domains = []
    for c in cells(f):
        domains.append(boundary_markers[c])

    domains = list(set(domains))
    if len(domains) > 1:
        boundaries[f] = 2
        markers[c]=True
        

mesh = refine(mesh, markers)

My questions are:

  1. The code above works if I use “generate_mesh()” to create a mesh. If I use “RectangleMesh()” or somethin else, then if len(domains) > 1: is always false. Why do we have this different behaviour?

Unfortunately, I have to use a structured mesh instead of the unstructered mesh coming from “create_mesh”, so I want this attempt to work by using “RectangleMesh(,“crossed”)”.

  1. The same thing occurs after refinement: If I refine the mesh one time, then we have again if len(domains) > 1: is always false. What’s happening here?

I know that my solution is by far not elegant but in principle it should work.

Many thanks in advance,

Alex

2

I found out that the problems may come from the relation from the domain and the mesh:
By using some meshing operation with respect to some points (RectangleMesh(), UnitSquareMesh()…) these meshes are not related to a special domain.
But if we use generate_mesh(my_domain,n_resolution) this mesh is related to “my_domain” and domain related things like MeshFunction(‘size_t’, mesh, 2, mesh.domains()) are working.

After a refinement process a new mesh is created which is no longer related to “my_domain” and therefore the above tools which are related to some domain won’t work.

Therefore the important question is: How is possible to relate a given mesh to some domain without creating a new mesh via generate_mesh()?

3

After further research I think I have to reinitialise the class MeshDomains() after refinement. But how should I do that? One naive try is to call MeshDomains.init() but this leads to the error

    MeshDomains.init()
TypeError: unbound method MeshDomains_init() must be called with MeshDomains instance as first argument (got nothing instead)
Aborted

Any ideas?
Thank you very much!

4

Hi,

I have a similar issue.
I want to make a domain with a rectangular and a (initial) sphere inside. The length of the rectangle (L=0.1) is larger than the height (H=0.01,0.01 for now). I need to use generate_mesh since I have a subdomain. The code is the following:

import dolfin
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
import matplotlib
from matplotlib import rc # use LaTeX text

# Plotting parameters
matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
matplotlib.rcParams['text.usetex'] = True
rc('text', usetex=True)
rc('font', family = 'serif')


# Problem parameters
T = 10.0                  # final time
num_steps = int(500)      # number of time steps
##num_steps = int(1.0e3)     # number of time steps
dt = T / num_steps        # time step size
mu = 1                # dynamic viscosity [N s/m2]
rho = 1                # density [kg/m3]
L = 1                  # channel length [m]
H = 1e-1                  # channel height [m]
r = 0.25e-2               # cell radius [m]
x0,y0 = 0.5e-0,0.5e-2    # cell centre [m]
Dp = 8                # pressure difference

# Define geometries
domain = Rectangle(Point(0,0),Point(L,H))
#circle = Circle(Point(-0.5,0.5),r)
circles = [[r,x0, y0]]; circles = np.array(circles)

for circle in circles:
    circleD = Circle(Point(circle[1:]), circle[0], 60)

# Set subdomains
domain.set_subdomain(1,circleD)

# Create mesh
mesh=generate_mesh(domain,60)
#markers = MeshFunction('size_t',mesh,2,mesh.domains())
coor = mesh.coordinates()
xcoor, ycoor = coor[:,0], coor[:,1]

# Track boundaries
class InnerBoundary_circle(SubDomain):
    def __init__(self, r, x0, y0):
        self.r  = r
        self.x0 = x0
        self.y0 = y0
        super().__init__()
    def on_circle_condition(self, x, y):
        return np.power(x-self.x0, 2) + \
               np.power(y-self.y0, 2) - \
               np.power(self.r, 2)
    def point_belongs_to_circle(self, x):
        x_, y_ = x[:]
        if near(self.on_ciricle_condition(x_, y_), 0., eps= 1.e-12):
            return True
        else: return False
    def inside(self, x, on_boundary):
        if not on_boundary:
            return on_boundary
        else:
            return self.point_belongs_to_circle(x)

# Open figure
plt.figure(figsize=(8,3),dpi=400)
colors = ('k', 'r')
# Axis
plt.xlabel(r'$x$', fontsize = 12)
plt.ylabel(r'$y$', fontsize = 12, rotation='horizontal')

xc,yc,idc=[],[],[]
# Plot
plt.plot(xcoor,ycoor,'ok',ms=1,label='Flow',zorder=1)
for i,(x,y) in enumerate(coor):
    if which_circle(x,y,circles)!=0:
        plt.scatter((x,), (y,),s=1,color=colors[which_circle(x,y,circles)],zorder=2)
        plt.plot([], [], 'or',ms=1,label='Cell',zorder=3)    
# Keep cell boundary nodes and indices
        xc.append(x)
        yc.append(y)
        result=np.where(xcoor==x)
        idc.append(result[0][0])
# Legend
plt.legend(bbox_to_anchor=(0.1,1.01),edgecolor='k',markerscale=4,ncol=2,fontsize=11)
# Save plot
plt.savefig('domain.pdf',dpi=400)
#plt.show()

The issue I have is that as I decrease H to 0.1 and 0.01, the mesh/cell points are decreasing between 0 and H. Therefore, I end up with no points “inside” the rectangle, just at the boundaries. Please, print the PDFs to see.

If I use RectangleMesh then I have points everywhere, however I cannot introduce a subdomain there.

I am very confused how generate_mesh works. Any ideas?

Many thanks,
Eva

5

I would strongly advice to use something else than mshr to create your mesh, as it is no longer maintained. I would suggest using pygmsh or the Gmsh (gui or Python API).
See for instance:

Or:
https://gmsh.info/

I have made a tutorial on how to use the Gmsh Python API:

And section 3 of

on how to create an xdmf file that can be loaded into dolfin

2 Likes
6

Hi dokken,

Thank you for these references.
I felt that the geometry I need is fairly simple, hence I was questioning why I get these differences.

Is there any way I can ad subdomains using the mesh methods you mentioned? I didn’t see anything regarding that in the tutorials.

Many thanks,
Eva

7

Gmsh supports physical tags, which can be loaded into dolfin, see for instance:
the second section of:
Mesh generation and conversion with GMSH 4.6.0 and PYGMSH 6.1.1 | Jørgen S. Dokken on how to create physical markers,
and
Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #7 by dokken
for how to get the physical tags from msh file to dolfin,
or Electromagnetics example — FEniCS-X tutorial on how to create phyiscal markers using the GMSH python API.

8

Hi,

I am trying to run pygmsh for a simple rectangle and I get an error for the following line:

from pygmsh.built_in.geometry import Geometry

the error:

from pygmsh.built_in.geometry import Geometry
ModuleNotFoundError: No module named 'pygmsh.built_in'

I have installed everything properly.

What can be the problem?

Thank you.

9

As the tutorial specifies, you have to use version 6.1.1.

If you want to use a more recent pygmsh version, I would suggest having a look at the demos in the pygmsh repository (GitHub - nschloe/pygmsh: Gmsh for Python)

11

There are quite a few differences which I mention here for someone who will face the same issues:

  1. Instead of pygmsh.built_in.geometry import Geometry, simply use geometry=pygmsh.geo.Geometry().
  2. No lcar is needed in the add_point so just add_point(x,x,x), resolution.
  3. add_lines changed to add_curves, same for any line command

For the following, please dokken, correct me if I am wrong or if the change I “found” will affect the code:

  1. In add_physical the label needs to be string therefore label=volume_marker instead of just volume_marker which is a number.

I have not found a solution for that though:

  1. In generate_mesh the attribute prune_z_0 does not exist. Do I need that still?
  2. I cannot export/save files within generate_mesh (with geo_filename="mesh.geo"). Is there any other way to save these files?
12

As these questions are directly related to interface changes in pygmsh, the two last questions are best to ask in the issue tracker for GitHub - nschloe/pygmsh: Gmsh for Python.
In general, you can remove the z-coordinate from any 2D mesh when you are converting it with meshio
(I.e. if the read in points with meshio are of the shape (num_points, 3), remove the last column before saving it to xdmf).
I personally stopped using pygmsh after v 6.1.1, as it now is more or less the same as the GMSH python API, which I have referenced in previous posts on this topic.

My suggested pipeline for new users is the following:

  1. Generate an msh file using either GMSH graphical interface or GMSH python interface.
  2. Convert the mesh from msh to XDMF using meshio.

Optionally, replace gmsh with pygmsh or any other mesh generator that can produce a file that is readable by meshio.

There are various mesh formats available for representing unstructured meshes. 
meshio can read and write all of the following and smoothly converts between them:

    Abaqus, ANSYS msh, AVS-UCD, CGNS, DOLFIN XML, Exodus, FLAC3D, H5M, Kratos/MDPA, Medit, MED/Salome, Nastran (bulk data), Neuroglancer precomputed format, Gmsh (format versions 2.2, 4.0, and 4.1), OBJ, OFF, PERMAS, PLY, STL, Tecplot .dat, TetGen .node/.ele, SVG (2D only, output only), SU2, UGRID, VTK, VTU, WKT (TIN), XDMF.
13

In the case of gmsh, I try to load the msh file I created in fenics but I get the following error:

from dolfinx.io import extract_gmsh_geometry, extract_gmsh_topology_and_markers, ufl_mesh_from_gmsh ModuleNotFoundError: No module named 'dolfinx'

I have downloaded gmsh_helpers.py and I am trying to run

import dolfin
from fenics import *
from mshr import *
from gmsh_helpers import read_from_msh
import numpy as np
import matplotlib.pyplot as plt

# Mesh import
mesh, cell_tags, facet_tags = read_from_msh("mesh.msh", cell_data=True, facet_data=True, gdim=2)

Regarding the import and manipulation of meshes from gmsh to fenics, do I need to use meshio etc as I see on other threads here, or simply reading it like this will implement the boundaries etc?

I also haven’t seen any tutorial about subdomains using gmsh. Please correct me if I am wrong.

Thanks!

14

Dolfin and Dolfin-x are two different versions of FEniCS, and the do not have the same functionality.
For dolfin (fenics) the only option for loadning a mesh from GMSH is by writing an msh file and using meshio to do this.
For dolfinx it is possible to directly interface with the gmsh API and convert it without loading (which is what you are trying to do).

As I mentioned above, the following tutorial (Using dolfin-x for simulation, but the GMSH python API for generating the mesh) uses subdomains
https://jorgensd.github.io/dolfinx-tutorial/chapter3/em.html

15

I am trying to implement the mesh from this tutorial . Is this not for fenics? So do I need to use meshio?

Thanks!

16

Part one of the tutorial is a pure GMSH tutorial.
Part two of the tutorial, with title 2. How to load this mesh into dolfin-X without IO is on how to load the mesh into dolfin-x without meshio. As you are using dolfin, you need to use meshio`, as shown Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #7 by dokken

17

Okay, I assumed these tutorials were for fenics.

I created a msh mesh file using gmsh, however when I try to read it with mesh=meshio.read("mesh.msh"), I get the following:

File "cell_gmsh.py", line 57, in <module> mesh=meshio.read("mesh.msh") File "/home/user/.local/lib/python3.8/site-packages/meshio/_helpers.py", line 69, in read return reader_map[file_format](filename) File "/home/user/.local/lib/python3.8/site-packages/meshio/gmsh/main.py", line 19, in read mesh = read_buffer(f) File "/home/user/.local/lib/python3.8/site-packages/meshio/gmsh/main.py", line 48, in read_buffer return reader.read_buffer(f, is_ascii, data_size) File "/home/user/.local/lib/python3.8/site-packages/meshio/gmsh/_gmsh41.py", line 109, in read_buffer return Mesh( File "/home/user/.local/lib/python3.8/site-packages/meshio/_mesh.py", line 52, in __init__ assert len(data) == len(cells), ( AssertionError: Incompatible cell data. 9 cell blocks, but 'gmsh:physical' has 4 blocks.

Did I create a wrong msh file?

import gmsh
import meshio

gmsh.initialize()
gmsh.model.add("channel")

# Channel parameters
L = 100
H = 1

# Add points
lc=1e-1
p1=gmsh.model.geo.addPoint(0,0,0,lc)
p2=gmsh.model.geo.addPoint(0,H,0,lc)
p3=gmsh.model.geo.addPoint(L,H,0,lc)
p4=gmsh.model.geo.addPoint(L,0,0,lc)

# Add lines
l1=gmsh.model.geo.addLine(p1,p2)
l2=gmsh.model.geo.addLine(p2,p3)
l3=gmsh.model.geo.addLine(p3,p4)
l4=gmsh.model.geo.addLine(p4,p1)

# Define surface
cl1=gmsh.model.geo.addCurveLoop([l1,l2,l3,l4])
s1=gmsh.model.geo.addPlaneSurface([1])

gmsh.model.geo.synchronize()

# Define boundaries
inlet_marker, outlet_marker, walls_marker = 1,2,3
gmsh.model.addPhysicalGroup(s1,[l1],inlet_marker)
gmsh.model.setPhysicalName(s1,inlet_marker,"Inlet")
gmsh.model.addPhysicalGroup(s1,[l3],outlet_marker)
gmsh.model.setPhysicalName(s1,outlet_marker,"Outlet")
gmsh.model.addPhysicalGroup(s1,[l2,l4],walls_marker)
gmsh.model.setPhysicalName(s1,walls_marker,"Walls")

# Generate and write mesh
gmsh.model.mesh.generate(2)
gmsh.option.setNumber("Mesh.SaveAll", 1)
gmsh.write("mesh.msh")
`
18

You should probably not use SaveAll

19

But that will save only the points in the boundaries, according to this tutorial. Don’t I need all the elements inside the rectangle too, in order to solve the problem?

20

Therefore, create a physical marker for each entity you would like to save, i.e. Add a Physical group for each surface you would like to mesh.
For instance, to mesh a simple disk:

import gmsh
gmsh.initialize()
membrane = gmsh.model.occ.addDisk(0, 0, 0, 1, 1)
gmsh.model.occ.synchronize()
gdim = 2
status = gmsh.model.addPhysicalGroup(gdim, [membrane], 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.05)
gmsh.model.mesh.generate(gdim)
21

Great, thank you for that!

I am going through the thread and I have a few questions:

  1. How can I tell if my mesh has triangles or tetrahedrons? I assume it is something on the gmsh creation?
  2. I use the create_mesh function presented here. Is there any way I can get the mesh node coordinates?
  3. Why do we need two xdmf filesL the triangle and the facet (cell type=lines)?
  4. It is not clear to me how I can read the mesh and use the markers I made in gmsh in order to add boundary conditions. Later, I will need to add a sudbomain in my mesh and problem, so I think that will answer my future question too.

Many thanks!