Flag cells to be refined more than one levels?

[Greg Wang]

Dear all,

This might be a somewhat odd feature request but I was wondering whether the functionality described below can be readily achieved in deal.II already:

Say we are given a coarse mesh with every cell never refined, and a vector of non-negative integers associated with each cell. The integer is the level of isotropic refinement to be applied to each cell. The number can be 0, which means not doing anything here; it can be 1, which corresponds to what set_refine_flag() does. But it can also be any number bigger than 1. Say we have a quadrilateral cell and the number associated with it is 2, then the cell is to be refined isotropically once into four children and then every child is refined again, resulting in 16 grandchildren.

If we were to hard code this, would it be sound to simply rely on set_refine_flag() and execute the refinement one level after another for M levels, with M being the maximum number of the given vector? And perhaps to do this we would need iterators triangulation.begin(n), triangulation.end(n) and isotropic_child() involved somehow?

Any help and/or pro tips would be greatly appreciated!

Best regards,
Yuan

​

[Peter Munch]

Hi Greg,

I guess what you could do is to treat the vector as vector associated with a DoFHandler with FE_DGQ(0). During each refinement you "interpolate" the vector to the new mesh and decrease the value by one. This approach naturally extends to a parallel setting.

Alternately, you could create a single vector associated the coarse cells and while looping over active cells of the refined mesh you would recursively call the parent cell until you reach the coarsest cell whose id you can use to access the vector. In the parallel setting, the vector associated to the coarse cells would need to replicated between processes.

Hope this help,

PM

[Yuan Wang]

Hi Peter,

Thanks a lot for your help! I was trying to implement a MWE with the first method you suggested. However, the function VectorTools::interpolate_to_different_mesh() seems to be causing a linker error:

/home/usrname/Projects/hdg_dealii/adr/ip/multilevel_refinement/MultilevelRefinement.cc:66: error: undefined reference to 'void dealii::VectorTools::interpolate_to_different_mesh<2, 2, dealii::Vector<int> >(dealii::DoFHandler<2, 2> const&, dealii::Vector<int> const&, dealii::DoFHandler<2, 2> const&, dealii::Vector<int>&)' collect2: error: ld returned 1 exit status make[2]: *** [CMakeFiles/MultilevelRefinement.dir/build.make:133: MultilevelRefinement] Error 1 make[1]: *** [CMakeFiles/Makefile2:90: CMakeFiles/MultilevelRefinement.dir/all] Error 2 make: *** [Makefile:91: all] Error 2

Here is, I think, the relevant chunk of code (original code file is also attached):

// ... #include <deal.II/numerics/vector_tools.h> // ... Triangulation<dim> triangulation_old; triangulation_old.copy_triangulation(triangulation); DoFHandler<dim> dof_handler_old(triangulation_old); dof_handler_old.distribute_dofs(fe); Vector<int> ref_level_old(ref_level); triangulation.execute_coarsening_and_refinement(); dof_handler.reinit(triangulation); dof_handler.distribute_dofs(fe); ref_level.reinit(dof_handler.n_dofs()); VectorTools::interpolate_to_different_mesh(dof_handler_old, ref_level_old, dof_handler, ref_level); // ...

Can you help shed some light on this?

Best regards,

Yuan

​

[Yuan Wang]

Here is perhaps a minimaler (non)working example:

#include <deal.II/grid/tria.h> #include <deal.II/grid/grid_generator.h> #include <deal.II/grid/grid_refinement.h> #include <deal.II/dofs/dof_handler.h> #include <deal.II/dofs/dof_tools.h> #include <deal.II/fe/fe_dgq.h> #include <deal.II/lac/vector.h> #include <deal.II/numerics/vector_tools.h> using namespace dealii; int main() { FE_DGQ<2> fe(0); Triangulation<2> trig; Triangulation<2> trig_child; GridGenerator::hyper_cube(trig); trig.refine_global(1); trig_child.copy_triangulation(trig); trig_child.refine_global(1); DoFHandler<2> dof_handler(trig); dof_handler.distribute_dofs(fe); Vector<int> ref_level({1,2,3,4}); // ref_level is the vector on // the coarser mesh trig DoFHandler<2> dof_handler_child(trig_child); dof_handler_child.distribute_dofs(fe); Vector<int> ref_level_child(dof_handler_child.n_dofs()); // ref_level_child is the vector // on the finer mesh trig_child VectorTools::interpolate_to_different_mesh(dof_handler, ref_level, dof_handler_child, ref_level_child); }

And I just found that the bug comes from feeding Vector<int> to
interpolate_to_different_mesh() while int is not a valid
VectorType::value_type therein. Changing int to double fixes it.

Best,
Yuan

​

[Wolfgang Bangerth]

On 11/10/23 17:43, Yuan Wang wrote:
> And I just found that the bug comes from feeding Vector<int> to
> interpolate_to_different_mesh() while int is not a valid
> VectorType::value_type therein. Changing int to double fixes it.

Yes, that is right. We do not seem to enforce this, but we are generally
expecting Vector<T> to only be instantiated for T=some floating point type.

Best
W.

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

[Peter Munch]

Hi Greg,

good that you find a working solution! However, I would suggest against using `VectorTools::interpolate_to_different_mesd` but use `SolutionTransfer` or `parallel::distributed::SolutionTransfer` to move the refinement information from a coarse mesh to a fine mesh. You can take a look at the tutorial how these classes are used.

Hope this helps,
Peter

[Yuan Wang]

Hi Peter,

Thanks for the suggestion! I had an example using SolutionTransfer worked out
and only then realized the premise of this whole feature request missed the mark…

So let me restate what is desired: I basically want to record the refinement
(and coarsening) history of a Triangulation and reconstruct the mesh. There
are already available save_refine_flags and load_refine_flags functions
that are precisely for this but it’s not really suitable for my case because I
need to reconstruct the mesh at every time-step and it’s better to avoid I/O
operations since they tend to slow things down.

In theory, it should work by recording level and index numbers of all cells
that has_children() and then execute_coarsening_and_refinement level by
level on the same coarse mesh. However, oddities of cell indexing were
observed. See attached mesh as an example. On level 1, I was expecting to see
indices in the order of 0, 1, 2, 3, 6, 7, etc but instead saw the output below
(the third column is the barycenter coordinate of the cell). The cells are not
lexicographically ordered on level 1 despite being lexicographically refined
when the reconstruction part was carried out. This results in different meshes.

Level: 0; Index: 0; (-8 -8) Level: 0; Index: 1; (8 -8) Level: 0; Index: 2; (-8 8) Level: 0; Index: 3; (8 8) Level: 1; Index: 0; (-12 -12) Level: 1; Index: 1; (-4 -12) Level: 1; Index: 2; (-12 -4) Level: 1; Index: 3; (-4 -4) Level: 1; Index: 4; (-12 4) Level: 1; Index: 5; (-4 4) Level: 1; Index: 6; (-12 12) Level: 1; Index: 7; (-4 12) Level: 1; Index: 8; (4 4) Level: 1; Index: 9; (12 4) Level: 1; Index: 10; (4 12) Level: 1; Index: 11; (12 12) Level: 1; Index: 14; (4 -4) Level: 1; Index: 15; (12 -4) Level: 2; Index: 32; (2 2) ...

And here is the chunk of my code that does the recording:

std::map<int, std::vector<int>> cell_map; for (auto &cell : triangulation.cell_iterators()) { if (cell->has_children()) { cell_map[cell->level()].push_back(cell->index()); } }

…then the reconstruction:

for (auto cell_data = cell_map.begin(); cell_data != cell_map.end(); ++cell_data) { int cell_level = cell_data->first; std::vector<int> cell_indices = cell_data->second; for (auto cell_index : cell_indices) { CellAccessor<dim, dim> cell(&triangulation_copy, cell_level, cell_index); cell.set_refine_flag(); } triangulation_copy.execute_coarsening_and_refinement(); }

See also the picture of the reconstructed mesh. I believe that level 1 cell 12
was refined there due to 1-irregularity constraint.

Would appreciate your comments and help!

Best,
Yuan

​

[Wolfgang Bangerth]

On 11/21/23 12:59, Yuan Wang wrote:
> So let me restate what is desired: I basically want to record the refinement
> (and coarsening) history of a Triangulation and reconstruct the mesh.

The PersistentTriangulation is your class :-)

[Yuan Wang]

Hi Wolfgang,

This is a massive time-saver thanks so much!

Although I’m wondering how well, if at all, PersistentTriangulation would
work in the distributed parallel setting. The documentation and the inheritance
diagram of Triangulation are not looking very encouraging.

So I made a small example with the coarse_grid being a
parallel::distributed::Triangulation<dim>. And then the main work of the code
is done to PersistentTriangulation<dim> triangulation(coarse_grid). Things
worked out pretty fine actually! Although I didn’t do much other than moving
vertices around and setting refinement/coarsening flags - so there is no
dofhandler or finite element involved yet.

I think I’ll go ahead and apply this to my main project code and pray for the
best but I was wondering perhaps you would have a word on this matter?

Best regards,
Yuan

​

[Wolfgang Bangerth]

> Although I’m wondering how well, if at all, PersistentTriangulation would
> work in the distributed parallel setting. The documentation and the inheritance
> diagram of Triangulation are not looking very encouraging.

Yes, you are correct.

> So I made a small example with the coarse_grid being a
> parallel::distributed::Triangulation<dim>. And then the main work of the code
> is done to PersistentTriangulation<dim> triangulation(coarse_grid). Things
> worked out pretty fine actually!

Right. But the triangulation you restore this way is not parallel.

> I think I’ll go ahead and apply this to my main project code and pray for the
> best but I was wondering perhaps you would have a word on this matter?

If all you want is to store and re-create a triangulation, what you are
looking for is the concept of "serialization/deserialization". There is no
tutorial program (yet) that shows how this is supposed to work, but you can
search through the archives of the mailing list to find discussions about how
to do that in a number of posts.

You should also look at a draft tutorial program here:
https://github.com/dealii/dealii/pull/12429