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
--
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.
grid-0.eps
Description: PostScript document
grid.eps
Description: PostScript document
