Change in the behavior of Boundary Id interpretation from GMSH 2 to GMSH 4?

[Bruno Blais]

Hello,
I have been having some issues when trying to run a "toy" code (a simple Steady Navier-Stokes solver) on a case I used to run.
The boundary group Id 0, which is manually set in my GMSH script (vonKarman.geo), is detected perfectly as both a BC and a manifold when I export to GMSH 2 and I get the expected results.
However, if I export the same GMSH file to V4 the same code does not run anymore. I guess there is a change in behavior I did not notice in how dealii interprets the new format? What am I doing wrong?
The error I get is (when running in debug mode):

--------------------------------------------------------
An error occurred in line <10387> of file </home/blaisb/codes/dealii/dealii/source/grid/tria.cc> in function
    void dealii::Triangulation<dim, spacedim>::set_all_manifold_ids_on_boundary(dealii::types::boundary_id, dealii::types::manifold_id) [with int dim = 2; int spacedim = 2; dealii::types::boundary_id = unsigned int; dealii::types::manifold_id = unsigned int]
The violated condition was:
    boundary_found
Additional information:
    The given boundary_id 0 is not defined in this Triangulation!

Stacktrace:
-----------
#0  /home/blaisb/codes/dealii/build/lib/libdeal_II.g.so.9.1.0-pre: dealii::Triangulation<2, 2>::set_all_manifold_ids_on_boundary(unsigned int, unsigned int)
#1  ./schursteadynavierstokes: SchurNavierStokesSolver<2>::runVonKarman()
#2  ./schursteadynavierstokes: main
--------------------------------------------------------

However, BC ID 0 as a physical group is defined explicitely in the script. I have joined the .h and .cc of my code ( I am sorry it is ugly as hell, it was just made for a quick test).

The GMSH file is as follow (also in the tar archive):

y0=0;
y1=1;

x0 =0.;
x1 =10;

xc=1;
yc=0.5;
r=0.025;

lo = 2.0e-1;
lc = 1.0e-1;

Point(0) = {x0, y0, 0, lo};
Point(1) = {x0, y1, 0, lo};
Point(2) = {x1, y0, 0, lo};
Point(3) = {x1, y1, 0, lo};

Point(4) = {xc, yc, 0, lc};
Point(5) = {xc-r, yc, 0, lc};
Point(6) = {xc, yc+r, 0, lc};
Point(7) = {xc+r, yc, 0, lc};
Point(8) = {xc, yc-r, 0, lc};


// Contour Box
Line(1) = {0,1};
Line(2) = {1,3};
Line(3) = {3,2};
Line(4) = {2,0};

// Inner Circle
Circle(5)={5,4,6};
Circle(6)={6,4,7};
Circle(7)={7,4,8};
Circle(8)={8,4,5};

Line Loop(1) = {1,2,3,4};
Line Loop(2) = {5,6,7,8};

// Creates the physical entities
Physical Line(0)={5,6,7,8};
Physical Line(1)={1};
Physical Line(2)={2,4};
Physical Line(3)={3};

Plane Surface(1) = {1,2};
Physical Surface(2) = {1};
Recombine Surface{1};

The DEALII lines regarding the manifold are :

    grid_in.read_msh(input_file);
    Point<dim,double> circleCenter(1.0,0.5);
    static const SphericalManifold<dim> boundary(circleCenter);
    triangulation.set_all_manifold_ids_on_boundary(0,0);
    triangulation.set_manifold (0, boundary);

And the lines linked to the BC are :

emplate <int dim>
void SchurNavierStokesSolver<dim>::setup_dofs ()
{
    system_matrix.clear();
    pressure_mass_matrix.clear();

    dof_handler.distribute_dofs(fe);

    std::vector<unsigned int> block_component(dim+1, 0);
    block_component[dim] = 1;
    DoFRenumbering::component_wise (dof_handler, block_component);
    dofs_per_block.resize (2);
    DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
    unsigned int dof_u = dofs_per_block[0];
    unsigned int dof_p = dofs_per_block[1];

    FEValuesExtractors::Vector velocities(0);
    {
      nonzero_constraints.clear();
      DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints);
      VectorTools::interpolate_boundary_values(dof_handler,
                                               0,
                                               ZeroFunction<dim>(dim+1),
                                               nonzero_constraints,
                                               fe.component_mask(velocities));

      if (simulationCase_==TaylorCouette)
      {
          VectorTools::interpolate_boundary_values(dof_handler,
                                                   1,
                                                   RotatingWall<dim>(),
                                                   nonzero_constraints,
                                                   fe.component_mask(velocities));
      }
      if (simulationCase_==BackwardStep)
      {
          VectorTools::interpolate_boundary_values(dof_handler,
                                                   1,
                                                   PoiseuilleInlet<dim>(),
                                                   nonzero_constraints,
                                                   fe.component_mask(velocities));
      }



      if (simulationCase_==vonKarman)
      {
          std::set<types::boundary_id> no_normal_flux_boundaries;
          no_normal_flux_boundaries.insert (2);
          VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
                                                           no_normal_flux_boundaries,
                                                           nonzero_constraints
                                                           );
      }

      if (simulationCase_==vonKarman)
      {
          VectorTools::interpolate_boundary_values(dof_handler,
                                                   1,
                                                   ConstantXInlet<dim>(),
                                                   nonzero_constraints,
                                                   fe.component_mask(velocities));
      }
    }
    nonzero_constraints.close();

    {
      zero_constraints.clear();
      DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints);
      VectorTools::interpolate_boundary_values(dof_handler,
                                               0,
                                               ZeroFunction<dim>(dim+1),
                                               zero_constraints,
                                               fe.component_mask(velocities));

      if (simulationCase_==TaylorCouette || simulationCase_==BackwardStep || simulationCase_==vonKarman)
      {
          VectorTools::interpolate_boundary_values(dof_handler,
                                               1,
                                               ZeroFunction<dim>(dim+1),
                                               zero_constraints,
                                               fe.component_mask(velocities));
      }

      if (simulationCase_==vonKarman)
      {
          std::set<types::boundary_id> no_normal_flux_boundaries;
          no_normal_flux_boundaries.insert (2);
          VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
                                                           no_normal_flux_boundaries,
                                                           zero_constraints
                                                           );
      }


    }
    zero_constraints.close();
    std::cout << "   Number of active cells: "
              << triangulation.n_active_cells()
              << std::endl
              << "   Number of degrees of freedom: "
              << dof_handler.n_dofs()
              << " (" << dof_u << '+' << dof_p << ')'
              << std::endl;

}

[Wolfgang Bangerth]

Bruno,

I put this piece of code into your source file right after reading in the mesh:
|
for (const auto &cell : triangulation.active_cell_iterators())
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary())
std::cout << "cell=" << cell << ", face=" << f
<< ", boundary_id=" << cell->face(f)->boundary_id()
<< std::endl;
|
The output shows that all boundary faces have boundary ids between 1 and 8.

Of course you already know that from the error message. The question is how
these get there. I have to say that I don't understand the gmsh file format
version 4 -- it's quite complicated with the entities. But do you think you
could try to debug this problem by stepping into the GridIn and seeing where
and why the SubCellData objects are created that assign nonzero boundary_ids
to the boundary faces? The gmsh format 4 reader is new (just a few weeks) and
it would not be a surprise if it still had bugs in it. We would definitely
appreciate any help with this!

A good step would also be to create a smaller mesh that illustrates the issue.
The mesh you have has 400+ cells and 17 entities. It is difficult to read
through the mesh file and understand what each part does.

For all further work, the following, much-reduced program also shows the issue:
|
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_in.h>


int main ()
{
using namespace dealii;

const unsigned int dim = 2;

Triangulation<2> triangulation;

GridIn<dim> grid_in;
grid_in.attach_triangulation (triangulation);
std::ifstream input_file("vonKarman.msh");
grid_in.read_msh(input_file);

for (const auto &cell : triangulation.active_cell_iterators())
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary())
std::cout << "cell=" << cell << ", face=" << f
<< ", boundary_id=" << cell->face(f)->boundary_id()
<< std::endl;


triangulation.set_all_manifold_ids_on_boundary(0,0);
}
|

Best
W.

> Line(4)={2,0};

>
> // Inner Circle
> Circle(5)={5,4,6};
> Circle(6)={6,4,7};
> Circle(7)={7,4,8};
> Circle(8)={8,4,5};
>
> LineLoop(1)={1,2,3,4};

> LineLoop(2)={5,6,7,8};

>
> // Creates the physical entities
> PhysicalLine(0)={5,6,7,8};
> PhysicalLine(1)={1};
> PhysicalLine(2)={2,4};
> PhysicalLine(3)={3};
>
> PlaneSurface(1)={1,2};
> PhysicalSurface(2)={1};

> RecombineSurface{1};

>
> |
>
>
> The DEALII lines regarding the manifold are :
> |
>     grid_in.read_msh(input_file);
> Point<dim,double>circleCenter(1.0,0.5);
> staticconstSphericalManifold<dim>boundary(circleCenter);
>     triangulation.set_all_manifold_ids_on_boundary(0,0);
>     triangulation.set_manifold (0,boundary);
> |
>
> And the lines linked to the BC are :
> |
> emplate <intdim>

> voidSchurNavierStokesSolver<dim>::setup_dofs ()

> if(simulationCase_==TaylorCouette||simulationCase_==BackwardStep||simulationCase_==vonKarman)

> {
> VectorTools::interpolate_boundary_values(dof_handler,
> 1,
> ZeroFunction<dim>(dim+1),
>                                                zero_constraints,
>                                                fe.component_mask(velocities));
> }
>
> if(simulationCase_==vonKarman)
> {

>           std::set<types::boundary_id>no_normal_flux_boundaries;
>           no_normal_flux_boundaries.insert (2);
> VectorTools::compute_no_normal_flux_constraints (dof_handler,0,
>

>  no_normal_flux_boundaries,
>                                                            zero_constraints
> );
> }
>
>
> }
>     zero_constraints.close();
>     std::cout <<"   Number of active cells: "
> <<triangulation.n_active_cells()
> <<std::endl
> <<"   Number of degrees of freedom: "
> <<dof_handler.n_dofs()
> <<" ("<<dof_u <<'+'<<dof_p <<')'
> <<std::endl;
>
> }
>
> |
>
>

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.


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

[Bruno Blais]

Dear Wolfgang,

Thanks for looking into it.

I can't guarantee a timeline, but I'll try to take a look at it this weekend / early next week.

I'll use  your code to debug, it is indeed significantly smaller!. I'll try to see if I can come up with a very small mesh that breaks it down (like 10cells or so) in a similar fashion and try to work through it afterward.

Thanks for the help, i'll keep you posted :)!

Best,

Bruno

[Wolfgang Bangerth]

On 11/30/2018 02:06 PM, Bruno Blais wrote:
> Thanks for looking into it.
> I can't guarantee a timeline, but I'll try to take a look at it this
> weekend / early next week.
> I'll use  your code to debug, it is indeed significantly smaller!.

Yes, by a factor of 20 or 30 ;-))

> I'll
> try to see if I can come up with a very small mesh that breaks it down
> (like 10cells or so) in a similar fashion and try to work through it
> afterward.

I think that would be very useful. Let me know if you get stuck. I don't
have the time any more to debug all of these issues myself, but I'm
quite happy to help with advice! Together we'll figure this out!

Best
W.

[Daniel Arndt]

Bruno,

I created a corresponding issue in the GitHub repository at https://github.com/dealii/dealii/issues/7501.

Best,

Daniel

[Bruno Blais]

Daniel,

Thank you very much. I real the issue and saw that you guys found a smaller mesh to break the example.

I'll try to look further into it in the following weeks and see if I can come up with a nice solution.

Best

Bruno

[Daniel Arndt]

Bruno,

Can you try if the patch in https://github.com/dealii/dealii/pull/7637 fixes your problem?

Best,

Daniel

[Bruno Blais]

Dear Daniel,
Sorry I dropped the ball on this one.

I checked out the patch and tested it on two of my full cases and everything works well.
No differences between a GMSH ASCII 2 and GMSH ASCII 4, whereas before I had the previous issue.

From my point of view this fixes everything.

Thank you very much!