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

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.

Alex

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()?

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!

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, 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)
# 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

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:

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

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

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.

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.

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)

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.

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?

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

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 *
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!

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

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

Thanks!

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

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

# Channel parameters
L = 100
H = 1

lc=1e-1

# Define surface

gmsh.model.geo.synchronize()

# Define boundaries
inlet_marker, outlet_marker, walls_marker = 1,2,3
gmsh.model.setPhysicalName(s1,inlet_marker,"Inlet")
gmsh.model.setPhysicalName(s1,outlet_marker,"Outlet")
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")
```````

You should probably not use `SaveAll`

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?

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
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.05)
gmsh.model.mesh.generate(gdim)
``````

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!