Rectangle() and Circle() functions deprecated

1

Hi!

What would be the equivalent of

channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)

In the 2019.1.0 version?

2

Hi @luc ,
I am assuming you are trying to use mshr. Since it is outdated and possibly no longer maintained, I would recommend switching to an external mesh generator such as gmsh.

2 Likes
3

Well, I can do that. But now I’m stuck into converting the .msh file into .xml or .xdmf (apparently the last one is preferred).

First I converted this .geo file into .msh using Gmsh. But using

from fenics import *
dolfin-convert geometry.msh geometry.xml

gives me an “invalid syntax” error on the first argument. This answer prompts me to install and use meshio instead, which I did, but gives the same syntax error when using

import meshio
meshio-convert geometry.msh geometry.xdmf

My goal is to define a circular obstacle in a rectangular domain in order to study the scattering of an incident wave. Is it the right way of doing so?

4

dolfin-convert is deprecated too and xml is legacy format. You may use xdmf instead. I would use meshio-convert instead of dolfin-convert. If you have installed Gmsh and have correctly generated the mesh, then using meshio you can convert it via

import numpy as np
import meshio
msh = meshio.read("<filename>.msh")
clls = np.vstack(( 
    cell.data for cell in msh.cells if cell.type == "triangle"
    ))  # only write 2D cells unless you need lower dimensional elements in your code
meshio.xdmf.write("channel.xdmf", meshio.Mesh(msh.points, cells = {"triangle": clls}))

I normally use pygmsh as it is intuitive for simple problems. A complete code to generate the mesh you want would then be


import numpy as np
import pygmsh as pm
import meshio

geom = pm.opencascade.Geometry()
rect = geom.add_rectangle([0., 0., 0], 2.2, 0.41)

circ = geom.add_disk([0.2, 0.2, 0], 0.05)
final = geom.boolean_difference([rect], [circ])

msh = pm.generate_mesh(geom, dim=2, geo_filename="channel.geo")
mesh_to_write = meshio.Mesh(msh.points, cells = {
    "triangle" : np.vstack((
        cell.data for cell in msh.cells if cell.type == "triangle"
    ))
})

meshio.xdmf.write("channel.xdmf", mesh_to_write)

# load the mesh in dolfin
from dolfin import *
mesh = Mesh()
XDMFFile("channel.xdmf").read(mesh)

The geo file should look like

// This code was created by pygmsh<version>.

SetFactory("OpenCASCADE");

s0 = news;

Rectangle(s0) = {0.0, 0.0, 0, 2.2, 0.41};

s1 = news;

Disk(s1) = {0.2, 0.2, 0, 0.05};

bo1[] = BooleanDifference{ Surface{s0}; Delete; } { Surface{s1}; Delete;};

You can check how to install pygmsh and some examples here. Note that this is not a drop in replacement for the Gmsh API and it is simply hiding the mesh generation that you would manually do under the hood.

4 Likes
5

Thank you very much @bhaveshshrimali, it works.
I just noticed that I can no longer visualize my plots on my jupyter notebook though. It says

NotImplementedError: It is not currently possible to manually set the aspect on 3D axes

and displays an empty 3D plot. Does it have something to do with meshio or pygmsh ?

Last question if I may: how can I refine the mesh?

Thanks in advance!

6

No this is a matplotlib issue. I believe this was fixed in a later release but I haven’t checked. In any case, I would recommend using @marcomusy’s excellent library vedo. It’s dolfin.plot function gives very nice inline plots. You can install it via pip3 install vedo or via conda

from vedo.dolfin import plot as vplot
vplot(mesh) # or vplot(u) where `u` is your function to be plotted

For refining the mesh, either load the current mesh and use dolfin’s refine function, namely

mesh = refine(mesh)

or use equivalent function in Gmsh/pygmsh

Edit: I checked with dolfin-2019.2.0.dev0 and matplotlib-3.3.0 by plotting the result of Hyperelasticity demo. It can produce quiver plots but not the desired plot of displaced mesh. So I would go ahead with vedo in 3D unless someone knows of a better solution

3 Likes
7

Thank you again for the information! Vedo works great and I congratulate @marcomusy for the time and effort. Anyway, I’ll leave this link for those who are using Vedo in a jupyter notebook and may be facing the issue of the Vtk window which won’t close.

1 Like