Problems with gmsh 2D quad grids and deal.ii

[Johannes Reinhardt]

Hello everyone,

I encounter problems when trying to load 2D quad meshes generated by
gmsh into a deal.ii based code. I could roughly find out, what the
causes of the problems are, but I don't know whether it is a problem
with gmsh, my use of it, or a problem in deal.ii code. Therefore I have
sent this mail to the mailing lists of both gmsh and deal.ii. The rest
of this rather lengthy mail explains my findings in detail.

I need to calculate a certain quantity as a function of a change in the
geometry. For this reason I want to generate a large number of meshes
using gmsh, so I defined some functions in a common file, and then use
a python script to generate many short .geo files that just include the
common file and call the function with different parameters.

The meshes are 2D quad meshes, because my calculation, that is based on
deal.ii, which can only deal with quad meshes. The meshes are
relatively detailed already, because I do not use adaptive refinement
and need a good approximation of the curved boundaries of the various
surfaces.

However when I try to load the resulting .msh files in my calculation
(which uses deal.ii), some of the meshes fail to load. Of the
294 meshes I generated, 39 fail.

Of these 31 fail with the error message:

The violated condition was:
n_negative_cells==0 || n_negative_cells==cells.size()

and 8 fail with the error message:

The violated condition was:
needed_lines.find(std::make_pair(line_vertices.second,
line_vertices.first)) == needed_lines.end()

The first error seems to be connected to by a cell orientation
issue that leads to a invalid cell:
http://www.dealii.org/7.3.0/doxygen/deal.II/classGridReordering.html#ae3a479aa86fabf1f761ab05eb1e42838

The second seems to be due to a vertex ordering problem, I quote from
deal.II/source/grid/tria.cc around line 1640, where the error is thrown:

// assert that the line was
// not already inserted in
// reverse order. This
// happens in spite of the
// vertex rotation above,
// if the sense of the cell
// was incorrect.
//
// Here is what usually
// happened when this
// exception is thrown:
// consider these two cells
// and the vertices
// 3---4---5
// | | |
// 0---1---2
// If in the input vector
// the two cells are given
// with vertices <0 1 4 3>
// and <4 1 2 5>, in the
// first cell the middle
// line would have
// direction 1->4, while in
// the second it would be
// 4->1. This will cause
// the exception.

The first error seems to be rather common, when looking it up, I found
a tool to postprocess .msh files to resolve it:
https://code.google.com/p/tethex/

When using this tool with my meshes, the first error is indeed avoided,
and of my 294 meshes, but now 9 meshes fail with the second error, so
for one file the first error masked the second one.

I use deal.ii 7.3 and the Linux 64 bit release of gmsh 2.7.1 on Xubuntu
13.04 64bit.

I have attached some files from my investigations:

common_params.geo contains functions that are used to build the
geometry and common parameters.
generate_geo.py is used to generate individual .geos that include
common_params
make_meshes.sh is a shell script that sets up directories, runs the
python scripts, runs gmsh and runs tethex
mass_grid.cc is a small snippet that uses the deal.ii grid classes to
load meshes and logs what is going wrong
failure.log is a log about which .msh files can not be loaded with and
without using tethex and why.

I uploaded a selection of the msh files that should cover all the cases
here (~3MB):

http://ist-dein-freund.de/geometries.tar.gz

Thank you for your efforts in advance

Johannes

[Uwe Köcher]

Hi Johannes,

there is another amazing project tethex by Mikhael:

http://code.google.com/p/tethex/wiki/Tethex

I've tried this for a 2D circle grid (tet-grid from gmsh -> tethex -> deal.II) and this is working fine.

Maybe you try this first.

Best

  Uwe

[Mikhail Artemiev]

Hi,

I dug into Johannes' issue, and I can confirm that there is something strange here.
I used deal.II 7.2.0 and svn version (not the latest though) and gmsh 2.7.0 on Ubuntu 12.04.

For both versions (built in debug mode) the error was in function
    void dealii::internal::GridReordering2d::GridReordering::get_quads(std::vector<dealii::CellData<2> >&) const
The violated condition was:
    s[1] == s[3]
The name and call sequence of the exception was:
    ExcInternalError()
Additional Information:
(none)

Stacktrace:
-----------
#0  /home/artemiev/deal.II.dev.orig.new/lib/libdeal_II.g.so.8.0.pre: dealii::internal::GridReordering2d::GridReordering::get_quads(std::vector<dealii::CellData<2>, std::allocator<dealii::CellData<2> > >&) const
#1  /home/artemiev/deal.II.dev.orig.new/lib/libdeal_II.g.so.8.0.pre: dealii::internal::GridReordering2d::GridReordering::reorient(std::vector<dealii::CellData<2>, std::allocator<dealii::CellData<2> > >&)
#2  /home/artemiev/deal.II.dev.orig.new/lib/libdeal_II.g.so.8.0.pre: dealii::GridReordering<2, 2>::reorder_cells(std::vector<dealii::CellData<2>, std::allocator<dealii::CellData<2> > >&)
#3  /home/artemiev/deal.II.dev.orig.new/lib/libdeal_II.g.so.8.0.pre: dealii::GridIn<2, 2>::read_msh(std::istream&)
#4  ./proj: GridChecker::check_grid(std::string)
#5  ./proj: main
--------------------------------------------------------

@Johannes,
I got the same exception as you:

   The violated condition was:
      needed_lines.find(std::make_pair(line_vertices.second,

      line_vertices.first)) == needed_lines.end()

for 7.2.0 version built in release mode - did you also launch your program linking to the release version of deal.II ?

What is more interesting - for the case described in code:
 3---4---5
 |    |    |
 0---1---2
both versions works fine.
I checked it for the following .msh file:

$MeshFormat
2.0 0 8
$EndMeshFormat
$Nodes
6
1 0 0 0
2 1 0 0
3 2 0 0
4 0 1 0
5 1 1 0
6 2 1 0
$EndNodes
$Elements
2
1 3 2 0 1 1 2 5 4
2 3 2 0 2 5 2 3 6
$EndElements

Mikhail

[Johannes Reinhardt]

Hi,

thanks for looking into this issue.

Indeed, I had linked to the release library of deal.ii 7.3, when I link
to the debug lib I get the internal error in get_quads that you quoted.

Greetings

Johannes

> 2>std::allocator<dealii::CellData<2> > >&)
> 2>#3 /home/artemiev/deal.II.dev.orig.new/lib/libdeal_II.g.so.8.0.pre:
> 2>dealii::GridIn<2, 2>::read_msh(std::istream&) #4 ./proj:
> 2>GridChecker::check_grid(std::string) #5 ./proj: main
> 2>--------------------------------------------------------

--

[Wolfgang Bangerth]

> Of these 31 fail with the error message:
>
> The violated condition was:
> n_negative_cells==0 || n_negative_cells==cells.size()

This typically happens if you have inverted cells or cells that are not convex
in your mesh. This would be a result of what gmsh outputs.

(The test verifies that either no cells are inverted, or all are. If all are,
deal.II inverts all cells.)

> and 8 fail with the error message:
>
> The violated condition was:
> needed_lines.find(std::make_pair(line_vertices.second,
> line_vertices.first)) == needed_lines.end()

I don't know about this one (or the variant Mikhail posted -- thanks!). It
would be useful to have an input file that triggers this.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Johannes Reinhardt]

Hi,

> > Of these 31 fail with the error message:
> >
> > The violated condition was:
> > n_negative_cells==0 || n_negative_cells==cells.size()
>
> This typically happens if you have inverted cells or cells that are
> not convex in your mesh. This would be a result of what gmsh outputs.
>
> (The test verifies that either no cells are inverted, or all are. If
> all are, deal.II inverts all cells.)
>
>
> > and 8 fail with the error message:
> >
> > The violated condition was:
> > needed_lines.find(std::make_pair(line_vertices.second,
> > line_vertices.first)) == needed_lines.end()
>
> I don't know about this one (or the variant Mikhail posted --
> thanks!). It would be useful to have an input file that triggers this.

I uploaded a selection of the msh files that should cover all the cases
here (~3MB):

http://ist-dein-freund.de/geometries.tar.gz

Thank you

Johannes


--

[Wolfgang Bangerth]

>>> The violated condition was:
>>> needed_lines.find(std::make_pair(line_vertices.second,
>>> line_vertices.first)) == needed_lines.end()
>>
>> I don't know about this one (or the variant Mikhail posted --
>> thanks!). It would be useful to have an input file that triggers this.
>
> I uploaded a selection of the msh files that should cover all the cases
> here (~3MB):
>
> http://ist-dein-freund.de/geometries.tar.gz

That's too unspecific. As explained, the majority of the error results from
gmsh generating non-convex cells -- there is nothing deal.II can do about it.
Please post a file that shows the other kind of error. (In fact, even better:
open a deal.II bug report about it at
https://code.google.com/p/dealii/issues/list
and attach the mesh file, possibly along with a picture of the mesh for easier
orientation.)

Thanks

[Mikhail Artemiev]

Hello there!

@Johannes,
What I can offer to you is to use tethex for simplicial mesh.
I mean you can build triangular mesh (coarser than you used to build), and after that convert triangles into quadrangles (it makes the mesh finer).
I tested such technique for your files and scripts several times, and I didn't face the problem we are discussing here.
Perhaps, it could help as a transient solution.

Mikhail