mesh generation in step-6 in parallel::distributed::Triangulation and Trilinos

[Yiyang Zhang]

Hello guys.

I am trying to use the Grid generation mentioned in Step-6. It turns out that if we try to use it with parallel::distributed::Triangulation<>, things might be tricky.

So far, I notice that (i) I need to add a signal.post_refinement_connect to update the boundary indicators; all the MPI processes need to be notified when local cells are moved.

I have attached the relevant code below. Here are the questions:

i) the signal.post_refinement_connect() need to be set only when the grid is generated initially, then throughout the problem (even if we have AMR in the code), we don't need to set it again. Right?

ii) If I set manifold description for the mesh, do I need to set a similar signal.post_refinement_connect() to refresh the manifold description?

iii) the communicate_locally_moved_cells() is implemented in sync_moved_cells() function, is this implementation correct? I don't know how to modify the locally_owned_vertices vector (for example, the ith-component in this vector corresponds to which vertice), so I just plug all the locally_owned_vertices in.

iv) Is there anything missing? I am using Trilinos, do I need to make further modifications for this code?

v) A further question, according to parallel::distributed::Triangulation< dim, spacedim > Class Template Reference, the set_boundary_ids() is set to const. I am wondering why we can add the const-qualifier, since setting boundary indicator is changing tria (is it? which is a member of the class).

Thank you!

template<int dim>

    void MyClass<dim>::make_ball_mesh_2d(const Real radius, const int refine){

        //generate a hyper_ball tria with center at (0,0), and radius = radius

        const Real isCenter = 1e-3; //some small value

        const types::manifold_id MFD_ID_SPH = 1;

        GridGenerator::hyper_ball(tria, Point<dim>(), radius);

        //set up post_refinement signal

        //such that boundary indicators are refreshed after every refinement.

        //Ref. parallel::distributed::Triangulation<> class

        tria.signals.post_refinement.connect(std::bind(&MyClass<dim>::set_boundary_ids,

                                                       std::cref(*this),

                                                       std::ref(tria) ) );

        //set manifold id to get spherical manifold description

        tria.set_all_manifold_ids_on_boundary(MFD_ID_SPH);

        tria.set_manifold(MFD_ID_SPH, spherical_mfd);

        //set boundary ID

        set_boundary_ids(tria);

        const Point<dim> mesh_center;

        const Real core_radius = 1.0 / 100.0 * radius;

        const Real inner_radius = 1.0 / 50.0 * radius;

        // Step 1: Shrink the inner cell

        for (typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active();

             cell != tria.end(); ++cell) {

            if(cell->is_locally_owned()){

                if (mesh_center.distance(cell->center()) < isCenter)

                    //cell center almost on the mesh_center

                {

                    for (unsigned int v = 0;

                         v < GeometryInfo<dim>::vertices_per_cell;

                         ++v)

                        cell->vertex(v) *= core_radius / mesh_center.distance(cell->vertex(v));

                    //the vertices of this cell now have a radius of core_radius.

                }

            }

        }

        //when cells are moved locally, all the MPI processes need to be notified

        //by using void Triangulation< dim, spacedim >::communicate_locally_moved_vertices().

        sync_moved_cells();

        // Step 2: Refine all cells except the central one

        for (typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active();

             cell != tria.end(); ++cell) {

            if(cell->is_locally_owned()){

                if (mesh_center.distance(cell->center()) >= isCenter)

                    cell->set_refine_flag();

            }

        }

        tria.execute_coarsening_and_refinement();

        // Step 3: Resize the inner children of the outer cells

        for (typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active();

             cell != tria.end(); ++cell) {

            if(cell->is_locally_owned()){

                for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)

                {

                    const double dist = mesh_center.distance(cell->vertex(v));

                    if (dist > core_radius*1.0001 && dist < 0.9999*radius)

                        cell->vertex(v) *= inner_radius / dist;

                }

            }

        }

        sync_moved_cells();

        // Step 4: Apply curved manifold description

        for (typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active();

             cell != tria.end(); ++cell) {

            if(cell->is_locally_owned()){

                bool is_in_inner_circle = false;

                for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)

                    if (mesh_center.distance(cell->vertex(v)) < inner_radius)

                    {

                        is_in_inner_circle = true;

                        break;

                    }

                if (is_in_inner_circle == false)

                    cell->set_all_manifold_ids(MFD_ID_SPH);

            }

        }

        tria.refine_global(refine);

        return;

  }

    template<int dim>

   void MyClass<dim>::set_boundary_ids(parallel::distributed::Triangulation<dim>& tria_) const {

       for (typename Triangulation<dim>::active_cell_iterator

            cell = tria_.begin_active();

            cell != tria_.end();

            ++cell){

           if(cell->is_locally_owned()){

               for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)

                   if (cell->face(f)->at_boundary())

                       cell->face(f)->set_all_boundary_ids(BD_INFTY);

           }

       }

       return;

   }

   template<int dim>

   void MyClass<dim>::sync_moved_cells(){

       std::vector<bool> locally_owned_vertices = GridTools::get_locally_owned_vertices(tria);

       tria.communicate_locally_moved_vertices(locally_owned_vertices);

       return;

  }

[Bruno Turcksin]

Hi,

On Wednesday, October 18, 2017 at 4:42:03 PM UTC-4, Yiyang Zhang wrote:

i) the signal.post_refinement_connect() need to be set only when the grid is generated initially, then throughout the problem (even if we have AMR in the code), we don't need to set it again. Right?

No you still need it when you do refinement see here. The problem is that cells can be moved to a processor where you the boundary id was not set. If the boundary ids never change in your problem, you can set them everywhere at the beginning and then you don't need to do anything special when you refine. By everywhere, I mean not only cells that are not owned locally but also the cells that are not active.

ii) If I set manifold description for the mesh, do I need to set a similar signal.post_refinement_connect() to refresh the manifold description?

I am not sure but I would think so.

 

iii) the communicate_locally_moved_cells() is implemented in sync_moved_cells() function, is this implementation correct? I don't know how to modify the locally_owned_vertices vector (for example, the ith-component in this vector corresponds to which vertice), so I just plug all the locally_owned_vertices in.

I think it's ok. Does it look right?

iv) Is there anything missing? I am using Trilinos, do I need to make further modifications for this code?

Trilinos is used for the matrices, vectors, and the solver. The generation of the mesh is totally independent of these objects. You could use the same code with petsc.

v) A further question, according to parallel::distributed::Triangulation< dim, spacedim > Class Template Reference, the set_boundary_ids() is set to const. I am wondering why we can add the const-qualifier, since setting boundary indicator is changing tria (is it? which is a member of the class).

You are not changing the Triangulation but the cells that the Triangulation does not own but knows how to manipulate.

Best,

Bruno

[Yiyang Zhang]

Hello Bruno,

Thank you for your reply!

i) I checked this page. I am a little confused. As soon as signal.post_refinement.connect() is called, I assume set_boundary_ids() is always attached to the Triangulation object, and it will always refresh the boundary IDs when a refinement is done. Right?

Are you saying we need to call tria.signal.post_refinement.connect() every time before we call execute_mesh_refinement()? seems to be weird...

iii) I think it's ok. But should it be more efficient if we identify which cells are moved, and send a modified locally_owned_vertices vector? I don't know the implementation, so I don't know if that will change the efficiency.

Best,

Yiyang

[Bruno Turcksin]

Yiang,

2017-10-19 13:53 GMT-04:00 Yiyang Zhang <zyyce...@gmail.com>:

i) I checked this page. I am a little confused. As soon as signal.post_refinement.connect() is called, I assume set_boundary_ids() is always attached to the Triangulation object, and it will always refresh the boundary IDs when a refinement is done. Right?

Sorry I misread what you wrote. You are correct you only set it once but it is called many times.

 

iii) I think it's ok. But should it be more efficient if we identify which cells are moved, and send a modified locally_owned_vertices vector? I don't know the implementation, so I don't know if that will change the efficiency.

This is something I wouldn't bother with. Most of the time of your program is going to be spent in your solver. Doing the optimization that your are talking about will probably not make a big difference unless your mesh is very fine. Anyway if you want to look how to do it right, you can take a look at the tests (in particular this one show a way to fill the vector).

Best,

Bruno

[Yiyang Zhang]

Hello Bruno,

Thank you for your answer! Also thank you for the URL, those tests are really a good place to learn deal.II.

Best,

Yiyang