Mesh generation issue

[Peng Wei]

Hi all,

I was trying to adapt the mesh generation method introduced by Konstantin Ladutenko posted here: https://groups.google.com/forum/#!searchin/dealii/sphere$20mesh%7Csort:date/dealii/pvqNQUM3eyM/o0USYKbGBQAJ

But even though I ran exactly the same code, I got a messed-up mesh (see attached grid-0.eps). 

As far as my understanding (a beginner), the trouble-maker is the following piece of code:

// -------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 

 // skip the first and last elements in radii vector as far as they are already used.

            for (unsigned int i = 1; i < radii.size() - 1; ++i)

            {

                // Find tangental faces in cells and refine in tangental direction.

                for (auto cell : triangulation.active_cell_iterators()) {

                    if (grid_center.distance (cell->center())  < radii[i-1]) continue;

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

                    {

                        bool is_tangental = true;

                        const double dist = grid_center.distance (cell->face(f)->vertex(0));

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

                        {

                            double dist2 = grid_center.distance (cell->face(f)->vertex(v));

                            if (std::abs( dist2 - dist) > min_width)

                            {

                                is_tangental = false;

                                break;

                            }

                        }

                        if (is_tangental)

                        {

                            // Number of faces is two time larger than number of axis.

                            cell->set_refine_flag(RefinementCase<dim>::cut_axis(f/2));

                            break;

                        }

                    }  // end of for all faces

                }  // end of for all cells

                triangulation.execute_coarsening_and_refinement ();

                // Scale to current radius

                for (auto cell : triangulation.active_cell_iterators())

                {

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

                    {

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

                        if (dist > radii[i-1]*1.0001 && dist < outer_radius * 0.9999)

                            for (int j = 0; j < dim; ++j)

                                cell->vertex(v)(j) *= radii[i]/dist;

                    }  // end of for all vertices

                }  // end of for all cells

                Triangulation<dim>   tmp;

                GridGenerator::flatten_triangulation(triangulation,tmp);

                triangulation.clear();

                triangulation.copy_triangulation(tmp);

            }  // end of all radii

// -------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 

In particular, I tried with three layers (i.e., only three values in the input vector radii, std::vector<double> radii = {0.1, 0.9, 1.0};)

It seems that the upper-left cell is cut incorrectly, see also the attached grid.eps. The cut_axis is done with the following code:

                            cell->set_refine_flag(RefinementCase<dim>::cut_axis(f/2));

For that particular cell, f = 0. Thus I expect anisotropic cutting along x-axis. But in fact it is cut isotropically. 

I am using dealii version 8.5.0.

Many thanks for any advice.

Best,

Peng

[Jie Cheng]

Hi Peng

I've written something to generate sphere based on the documentation, which generates not as perfect, but usable mesh (see attached). If you are interested, click this.

Best
Jie