Creation of mesh with generate_mesh : most boundary points are not exactly on the boundary

[moritz braun]

Dear Friends

I have noted the following serious problem with the mesh generation using mshr in FEniCS. When creating a circle of radius r and then building a mesh usiing generate_mesh only 32 of the boundary points are exacty or within 1e-8 on the circle all of the remaining ones are off by +- 1% To reproduce the problem run python testCircle.py 10

This makes is impossible to do molecular muffin tine type calculations with this meshing software. I really don't feel like paying 100's of USDollars for cubit!

Any assistance will be much appreciated!

regards

Moritz

[Benjamin Dam Kehlet]

Hi Moritz

You are right, this is unfortunate, and a solution is work in progress, but requires a bit implementation work in mshr. Fortunately there is an easy workaround.

CGAL's 2D mesher works on polygons, so before invoking the mesher, mshr approximates the domain by a polygon. This polygon is what is actually being meshed be CGAL. When implemented properly, mshr should take the requested cell size into account when creating this polygon to make sure it is not too rough.

The workaround is to set the resolution of the polygon explicitly when creating the csg primitives with an additional optional argument to the constructor

c = Circle(Point(0.,0.), 1, 100)

Try experimenting with this argument. If set to high, the it will introduce short edge on the boundary so the mesh will be finer here.

Hope this saves you the money :-)

[moritz braun]

Dear Benjamin

Thanks !

I am just confused about which of the latter two arguments is the resolution and which controls the size of the polygonal side / surface for 3D

regards

Moritz

[moritz braun]

Thanks a lot

this is a good solution. It might be an idea to mention this additional argument somewhere!

regards

Moritz

[Benjamin Dam Kehlet]

The last argument to the constructor of the primitives (100 in the example above) controls the size of the polygonal side (setting it to 100 gives you a 100-sided polygon). (The second argument is the radius of the circle).

[moritz braun]

Dear Benjamin

Thanks a lot

I will now be able to make use of mshr

regards

Moritz

[Benjamin Dam Kehlet]

Good to hear!

(I'll add a note to the documentation about this)

[moritz braun]

Dear Benjamin

I am trying to understand the exact meaning of the resolution my guess is that there will be of the order of r^D elements if D is the dimension.

regards

Moritz

[moritz braun]

Dear Benjamin

After looking at 2D I was now trying the same thing for 3D but there seems to be a problem, because irrespective how I choose the third parameter, there are only 8 points that end up exactly on the boundary while the other ones are all a bit off but not as much as would happen on the cube corresponding to the 8 points. This currently rules out using meshr for 3D work

regards

Moritz

[Benjamin Dam Kehlet]

Hello Moritz

Is the domain a sphere? In that case you can use the class UnitSphereMesh in mshr

mesh = UnitSphereMesh(10)

The constructor takes one argument which behaves similarly to the resolution argument of generate_mesh() (ie. higher number gives more cells). This class basically does this: 1) Generate the boundary mesh of good quality with points exactly (ignoring roundoff errors) on the boundary. 2) Calls tetgen to mesh the volume with a flag that prevent inserting new points on the boundary.

Note that I just pushed some fixes to this class, so do pull from the repository before testing. The code may need some tuning. Let me know if something does not behave as expected.

The explanation that points are not exactly on the boundary when using the CGAL backend is that it insists on remeshing the surface. Hence it is not possible to control where the points on the boundary will be located. Since the domain that is being meshes is a polyhedral approximation (same as in 2D), the points will be slightly off the exact surface of the domain. We should add a note about this in the documentation.

[moritz braun]

Dear Benjamin

Thanks a lot!

I will try that.

However I will actually require a domain with the inner sphere being as regular and smooth as possible!

The other query I have , has do to with reading in surfaces in the object file format. Firstly it files for many files that geomview can read, and secondly I am not sure whether it does it properly.

Finally I am looking for a flag, to prevent mshr from refining a surface that is already fine enough to start with.

regards

Moritz

regards

Moritz

[moritz braun]

Dear Benjamin

I pulled from bitbucket but I do get the following error

StandardError Traceback (most recent call last) <ipython-input-3-bdcce947ecd2> in <module>()

[Benjamin Dam Kehlet]

• Not sure I understood what you mean by "regular and smooth". Could you elaborate a bit?
• Preserving the surface: This question comes back from time to time. This is possible with mshr, but you have to use the Tetgen backend. I added a short description here: https://bitbucket.org/fenics-project/mshr/wiki/FAQ
• Where did you get this error? Did mshr install successfully?

[Benjamin Dam Kehlet]

Regarding reading off files: Do you have a simple example file that fails? Mshr relies on CGAL for reading these files. I did experience once that CGAL was very picky about an extra lineshift in the file. We could then either report it to CGAL or write an OFF-file parser. The format is very simple, so will not be much work to write our own parser.

[Nico Schlömer]

Issue #45 was marked as a duplicate of this issue.

[Nico Schlömer]

Simply adding 100 to the geometry specs à la

import mshr
import dolfin

c = mshr.Circle(dolfin.Point(0., 0., 0.), 1, 100)
m = mshr.generate_mesh(c, 3)
dolfin.plot(m)
dolfin.interactive()

isn't really helping much in general:

If the number is chosen too small, the boundary vertices do not sit on the geometry; if too large, they mesh is too coarse at the boundary.

Somehow it must be kept in sync with the target refinement, but I don't know exactly how. As pointed out before, mshr should take care of this.

Anything I can do to help here?

[Nico Schlömer]

Another installment of this bug: This is a plot of the error of a finite volume solution of a PDE.

Exact Dirichlet boundary conditions are used, so you'd expect the error on the boundary to be of the order of machine precision. The contrary is true for every other point, and that's because the boundary points of the mesh don't sit exactly on, in this case, a circle.

[Benjamin Dam Kehlet]

This is on top of my todo-list for mshr now. This week is pretty busy, though.