problem with CylindricalManifold<3> after removing cells from a triangulation

[Simon]

Dear all,

 

I created a mesh using GridGenerator::hyper_cube_with_cylindrical_hole and removed after that half of the total number of the cells in order to make use of symmetry.

Then I set a SphericalManifold<2> (2D) respectively CylindricalManifold<3> (3D) in order to get a true circle (2D) respectively a true cylinder (3D) after global and/or adaptive refining.

Without removing some cells, i.e. just using the GridGenerator::hyper_cube_with_cylindrical_hole without using symmetry, both 2D and 3D variants work fine.

With removing some cells, i.e. using symmetry, only in 2D I get a true circle as result. In 3D however, the SphericalManifold<3> does not work anymore. (I can adaptively refine the inner hole, but the result is not a true cylinder anymore which was the case without removing the cells).

In the following is a snippet of my code. It´s hard for me to find out what the problem actually is, since in 2D the manifold description worked after removing some cells.

Thanks for helping!

Best

Simon

const double outer_radius = 1.0;
const double inner_radius = 0.5;
const Point<dim> center;

GridGenerator::hyper_cube_with_cylindrical_hole(triangulation_tmp,

                                                    inner_radius,
                                                     outer_radius,
                                                    0.5,
                                                    1,
                                                     false /*boundary_id_inner_hole is set to 1*/);

std::set<typename Triangulation<dim>::active_cell_iterator> cells_to_remove;
typename Triangulation<dim>::active_cell_iterator cell =triangulation_tmp.begin_active(),
        endc = triangulation_tmp.end();

for(; cell!=endc; cell++)
 {
      if(cell->center()[1]<0)
       {
            cells_to_remove.insert(cell);
        }

  }

GridGenerator::create_triangulation_with_removed_cells(triangulation_tmp, cells_to_remove, triangulation);

std::unique_ptr<Manifold<dim>> ptr_manifold=nullptr;

    
    if(dim==2)
    {
        ptr_manifold = std::make_unique<SphericalManifold<dim>>(center);
    }
    else if(dim==3)
    {
        ptr_manifold = std::make_unique<CylindricalManifold<dim>>(dim-1);
    }
    else
    {
        throw std::runtime_error("only allowed for dim == 2 or dim == 3");
    }
types::boundary_id boundary_id_inner_hole=1;
types::manifold_id manifold_id_inner_hole=1;
   triangulation.set_all_manifold_ids_on_boundary(boundary_id_inner_hole,manifold_id_inner_hole);

triangulation.set_manifold (manifold_id_inner_hole, *ptr_manifold);
    

..........

[Simon]

By the way, I forgot to precise my question in more detail:

In the code snippet above I set the corresponding Manifold Descriptions for dim=2 (spherical) and dim=3 (cylindrical) on the inner boundary created by  GridGenerator::hyper_cube_with_cylindrical_hole and removed after that some cells by calling GridGenerator::create_triangulation_with_removed_cells .

However, for the 3D-Case, after adaptively refining the inner hole there is still a FlatManifold attached to it and not a CylindricalManifold. For dim=2, I get the desired SphericalManifold and everything is fine.

I am pretty sure that I set the manifold_id on the inner boundary correctly, because after visualizing the mesh the cells are actually adaptively refined around the inner boundary, but as I said a FlatManifold is used in the 3D-Case.

-So does GridGenerator::hyper_cube_with_cylindrical_hole attach some important "information" in terms of id_s,... to my triangulation which my new triangulation created by  GridGenerator::create_triangulation_with_removed_cell does not have?

-Or does the CylindricalManifold only work for closed boundaries? Because without calling GridGenerator::create_triangulation_with_removed_cells also in the 3D-Case I get the desired CylindricalManifold.

Best

Simon

[luca.heltai]

There may be a problem with the way boundary ids are set.

Can you try the following?

after creating the grid with removed cells, loop over all cells and all faces, and if at boundary with boundary id == 1, then call

cell->face(f)->set_all_manifold_ids(1);

notice the “_all_”, i.e., it is not “set_manifold_id(1)”. This ensures you set manifold ids also on edges.

It may be that boundary ids are not set correctly on edges when creating the triangulation with removed cells. In this case this is a bug (simple to fix, btw), that you only see in 3d, because refinement starts from edges.

L.

> --
> 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.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/d326f1ec-97db-419c-a055-f727ebe55454n%40googlegroups.com.

[Simon]

Hello Luca,

your suggestion fixed my problem, so thank you very much!

Just to make sure that I understand it correctly:

In my code snippet from above I call the member function "triangulation.set_all_manifold_ids_on_boundary(1,1).

I actually thought that this call does the same thing as you suggested, i.e. cell->face(f)->set_all_manifold_ids(1).

But I guess the problem is that in my version I only set the manifold ids of the face and not its four edges as well, is it?

Best

Simon

[Luca Heltai]

The issue is with the function that removes the cells. It creates the grid correctly but forgets to set the ids of the edges. It is not with your code. For boundary ids this is usually not an issue, since interpolation and boundary conditions are set using the face ids (not the edges), but for manifold ids this is a problem.  

I'll open an issue, and use your test to check things. 

Thanks for coming up with such a minimal working example. :-)

Luca

Il giorno 26 apr 2021, alle ore 10:44, Simon <simon.w...@gmail.com> ha scritto:



To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/bf4e5915-84e6-4f29-89e4-5a29eef84652n%40googlegroups.com.