Triangulation of a cubic with inclusion

[Chenchen Liu]

Dear all,

I am dealing with a geometry that is a spare with a circular inclusion in 2D (see attached figure) and a cubic with a spherical inclusion in 3D. I have seen the function "hyper_cube_with_cylindrical_hole", but what I want is still solid inclusion rather than a hole. Is there anyway to build the mesh I want? Or I have to use the third party mesh software to do it? 

Thank you!

Chenchen

[Jean-Paul Pelteret]

Hi Chenchen,

You can indeed build a mesh for this geometry in deal.II. You just need to use multiple triangulations and some of the grid manipulation tools.

1. Build a triangulation (tria_1) with a custom set of subdivisions (3x3) using GridGenerator::subdivided_hyper_rectangle.

2. Build tria_2 that removes the central cell using GridGenerator::create_triangulation_with_removed_cells

 -- Note for 2-d only: An alternative to doing steps 1&2 would be to continue using GridGenerator ::hyper_cube_with_cylindrical_hole

3. Build tria_3 that represents the particle using GridGenerator::hyper_ball

4. Build the final tria that merges tria_2 and tria_3 using GridGenerator::merge_triangulations

5. Add a SphericalManifold manifold to the interface between the particle and surrounding matter. This ensures that upon grid refinement, the particle geometry is respected.

I hope that this helps!

J-P

[Sebastian.Gonzalez-Pintor]

Hi Chenchen

In addition to jean-Paul response, for 2d there is one example in the tutorial, but now I can not link to it because I cannot access the webpage. Nevertheless, you can do something as in the following link

https://groups.google.com/forum/#!topic/dealii/_uk1JDakBLc

You can refine this geometry, instead of using high order mappings, and you get what you are looking for. 

Sebas.

[Chenchen Liu]

Hi Jean-Paul,

Thank you for your reply. I tried the case for 2D, and what I have got is shown in the following,

The function "merge" can only merge two coarse meshes without refinement. It seems that tria_1 and tria_2 are not compatible. Will "manifold" work for this? I tried the following codes, but it does not work. In particular, I don't know how to set manifold at the interface of the two triangulations. Could you please take a look at it? Any suggestions would be helpful. Thanks!

void Step::make_grid ()  

{

    Point<2> center (0,0);

    GridGenerator::hyper_cube_with_cylindrical_hole ( tria_1, 0.25, 0.5);

    GridGenerator::hyper_ball (tria_2, center, 0.25);

    GridGenerator::merge_triangulations (tria_1, tria_2, triangulation);

    const SphericalManifold<2> boundary_description(center);

    triangulation.set_manifold (1, boundary_description);

    

    std::ofstream out ("grid-1.eps");   GridOut grid_out;   grid_out.write_eps (triangulation, out);

}

if I add the following codes, 

    Triangulation<2>::active_cell_iterator   cell = triangulation.begin_active(),  endc = triangulation.end();

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

         cell->set_all_manifold_ids (1);

I will get

Best,

Chenchen

在 2016年4月13日星期三 UTC-4上午2:50:14，Jean-Paul Pelteret写道：

[Jean-Paul Pelteret]

Ok, so GridGenerator::hyper_cube_with_cylindrical_hole can't be used in this case because of the differing discretisation at the coarsest level. Try the other approach where you remove a cell in the middle of a 3x3 grid. That definitely works.

[Chenchen Liu]

Hi Jean-Paul,

Thanks for your suggestion.

What I have obtained is strange for the outer meshes:

Still I am not sure about the manifold part. What I have done is to set manifold for boundary cells in the initial coarsest mesh. Please take a look at my code. Thanks!

void Step::make_grid ()  // Without Using GMESH to do it with inclusion.

{

    // Tria_1

    std::vector< unsigned int > repetitions(2);

    repetitions[0] = 3;  repetitions[1] = 3;

    GridGenerator::subdivided_hyper_rectangle (tria_1, repetitions,  Point<2>(-1.,-1.),Point<2>(1.,1.));

    

    // Tria_2

    std::set< typename Triangulation<2>::active_cell_iterator > cells_to_remove;

    Triangulation<2>::active_cell_iterator cell = tria_1.begin_active(), endc = tria_1.end();

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

    {

        if (cell->at_boundary ())

        {}

        else

            cells_to_remove.insert(cell);

    }

    GridGenerator::create_triangulation_with_removed_cells (tria_1, cells_to_remove, tria_2);

    

    // Tria_3

    Point<2> center (0,0);

    GridGenerator::hyper_ball (tria_3, center, 1./3.*sqrt(2.0));

    // Merge

    GridGenerator::merge_triangulations (tria_2, tria_3, triangulation);

  

    // Set manifold

    cell = triangulation.begin_active(), endc = triangulation.end();

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

        if (cell->at_boundary ())

            cell->set_all_manifold_ids (10);

    

    const SphericalManifold<2> boundary_description(center);

    triangulation.set_manifold (10, boundary_description);

    

    // Refine

    triangulation.refine_global (3);

    // Output 

    std::ofstream out ("grid-1.eps");

    GridOut grid_out;

    grid_out.write_eps (triangulation, out);

}

Best,

Chenchen

在 2016年4月13日星期三 UTC-4下午2:09:45，Jean-Paul Pelteret写道：

[Jean-Paul Pelteret]

Great, you've nearly got it! 

What's gone wrong here is that you've assigned the manifold geometry to the eight outer cells and all of their faces. This means that upon refinement all of the newly generated cell vertices are projected to the manifold. This is not quite what you want in this case. What you really want is just for the faces at the interface between the particle and the surroundings to conform to this geometry. So, instead of setting the manifolds of all the outer cells (and therefore their faces as well) with cell->set_all_manifold_ids(), what you need to do is only assign the manifold to the appropriate cell faces with cell->face(f)->set_all_manifold_ids().

[Chenchen Liu]

Hi Jean-Paul,

Thanks for your instruction. You are always helpful. I am now able to change the radius of the particle.

Best,

Chenchen

在 2016年4月14日星期四 UTC-4上午2:06:05，Jean-Paul Pelteret写道：