[deal.II] confusing result when read, refine, output, and re-read a mesh in gmsh format

[曾元圆]

hi all, I met a confusing problem when I try to create a mesh as follows:

1. I create a coarse mesh using gmsh. I used number 0, 1, 2  as 
boundary_indicators to denote three different boundary. The resulting 
square.msh looks good.

<https://lh3.googleusercontent.com/-eMdafmTXRrY/WubIVIjZ7_I/AAAAAAAAADM/JqHew64eO0AGaev3RvhmQPHCau3C0IlmACLcBGAs/s1600/square.msh.png>

2. Then I read the square.msh into dealii, and refine those cells with 
boundary_indicator 0. Then output again (I used GridOutFlags::Msh() flags 
to keep the boundary indicators  in the output file). The resulting 
square_refined0.msh still looks good.

<https://lh3.googleusercontent.com/-V7Rn2KaWX0E/WubIZyvg4-I/AAAAAAAAADQ/epK8TfiXuQ89bghrJVrTJ71Yh-P_HpifgCLcBGAs/s1600/square_refined0.msh.png>

3.  I want to refine the mesh more, so I re-read the square_refined0.msh 
again, run the same code (still refine those with boundary_indicator 0), 
and output the mesh again.  But this time, the resulting 
square_refined1.msh is no longer the same as I expected.

<https://lh3.googleusercontent.com/-KCj4yo49bo4/WubIg-_z_uI/AAAAAAAAADU/3CSNdmD3PwUKgybHXeyo338vWiXhmjlQwCLcBGAs/s1600/square_refined1.msh.png>

Apparantly this final mesh is not what I wanted, because I only want those 
active cells at boundary be refined.

I know in dealii, the default boundary_indicator of boundary faces is 0. 
But from this resulting mesh, it seems that the interior faces are also 
bound with indicators 0. I don't know why this happened. 

I do further experienments: I changed the boundary_indicator of the left 
and right side from 0 to 3, repeat above steps, and this time the resulting 
mesh is good:

<https://lh3.googleusercontent.com/-BInHxOcOgzA/WubKywGJ63I/AAAAAAAAADo/tRoGK6ByyegMG2r0cNX_ON0tXq3xo3AVgCLcBGAs/s1600/square_good_refined1.msh.png>

It seems that if I use 0 as boundary_indicator, then when dealii read in 
the .msh fille, it will lose some information  related to the bounday. So 
do I have to use non-zero boundary indicators? Or is this a bug? 

Someone explain it a little bit? 





My code is like:

void test_square(){

   Triangulation<2> triangulation;

   GridIn<2> gridin;

   gridin.attach_triangulation(triangulation);


   std::ifstream f("square_refined0.msh");

   gridin.read_msh(f);

   

   typename Triangulation<2>::active_cell_iterator 

            cell = triangulation.begin_active(),

            endc = triangulation.end();

   for(; cell!=endc; cell++){

      for(unsigned int face = 0; face<GeometryInfo<2>::faces_per_cell; 
face++){

         if(cell->face(face)->boundary_id() == 3)

            cell->set_refine_flag();

      }

   }

   triangulation.execute_coarsening_and_refinement();


   std::ofstream out("square_refined1.msh");    

   GridOut grid_out;

  

   grid_out.set_flags(GridOutFlags::Msh(true,false));     //write out the 
boundary indicator explicitly

   grid_out.write_msh(triangulation, out);

}

 



 

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

#include<deal.II/grid/tria.h>
#include<deal.II/grid/tria_accessor.h>
#include<deal.II/grid/tria_iterator.h>
#include<deal.II/grid/grid_generator.h>
#include<deal.II/grid/grid_tools.h>
#include<deal.II/grid/manifold_lib.h>
#include<deal.II/grid/grid_in.h>
#include<deal.II/grid/grid_out.h>

#include<iostream>
#include<fstream>

#include<map>

using namespace dealii;

void test_square(){
   Triangulation<2> triangulation;
   GridIn<2> gridin;
   gridin.attach_triangulation(triangulation);
   std::ifstream f("square.msh");
   gridin.read_msh(f);

   typename Triangulation<2>::active_cell_iterator
            cell = triangulation.begin_active(),
            endc = triangulation.end();
   for(; cell!=endc; cell++){
      for(unsigned int face = 0; face<GeometryInfo<2>::faces_per_cell; face++){
         if(cell->face(face)->boundary_id() == 0)
            cell->set_refine_flag();
      }
   }
   triangulation.execute_coarsening_and_refinement();

   std::ofstream out("square_refined0.msh");
   GridOut grid_out;

   grid_out.set_flags(GridOutFlags::Msh(true,false));
   grid_out.write_msh(triangulation, out);
}

int main(){
	try{
		   test_square();
	}
	catch(std::exception &exc){
		std::cerr<<std::endl<<std::endl<<"--------------"<<std::endl;
		std::cerr<<"exception on processing:"<<std::endl<<exc.what()<<std::endl;
		return 1;
	}
	catch(...){
		std::cerr<<std::endl<<"---------------"<<std::endl;
		std::cerr<<"unknown exception"<<std::endl;
		return 1;
	}
}

square.msh
Description: Mesh model