Skip to content

Commit

Permalink
Provide internal function for faster access of cell(3d)->line_index
Browse files Browse the repository at this point in the history
  • Loading branch information
@kronbichler
kronbichler committed Jun 8, 2022
1 parent f262a26 commit d4d2cd4
Showing 1 changed file with 211 additions and 109 deletions.
320 changes: 211 additions & 109 deletions source/grid/tria.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1807,6 +1807,65 @@ namespace internal
*/
struct Implementation
{
/**
* A helper function to provide faster access to cell->line_index() in 3D
*/
template <int dim, int spacedim>
static std::array<unsigned int, 12>
get_line_indices_of_cell(
const typename Triangulation<dim, spacedim>::cell_iterator &)
{
Assert(false, ExcImpossibleInDim(dim));
return {};
}

/**
* A helper function to provide faster access to cell->line_index() in
* 3D, given that the regular function is pretty slow.
*/
static std::array<unsigned int, 12>
get_line_indices_of_cell(const Triangulation<3, 3>::cell_iterator &cell)
{
std::array<unsigned int, 12> line_indices;
// Can apply a more efficient function for hexes
const auto ref_cell = cell->reference_cell();
if (ref_cell == ReferenceCells::Hexahedron)
{
for (unsigned int f = 4; f < 6; ++f)
{
const unsigned char orientation =
cell->get_triangulation()
.levels[cell->level()]
->face_orientations[cell->index() * 6 + f];
const std::array<unsigned int, 4> my_indices{
{ref_cell.standard_to_real_face_line(0, f, orientation),
ref_cell.standard_to_real_face_line(1, f, orientation),
ref_cell.standard_to_real_face_line(2, f, orientation),
ref_cell.standard_to_real_face_line(3, f, orientation)}};
for (unsigned int l = 0; l < 4; ++l)
line_indices[4 * (f - 4) + l] =
cell->quad(f)->line_index(my_indices[l]);
}
for (unsigned int f = 0; f < 2; ++f)
{
const unsigned char orientation =
cell->get_triangulation()
.levels[cell->level()]
->face_orientations[cell->index() * 6 + f];
const std::array<unsigned int, 2> my_indices{
{ref_cell.standard_to_real_face_line(0, f, orientation),
ref_cell.standard_to_real_face_line(1, f, orientation)}};
line_indices[8 + f] = cell->quad(f)->line_index(my_indices[0]);
line_indices[10 + f] = cell->quad(f)->line_index(my_indices[1]);
}
}
else
for (unsigned int l = 0; l < std::min(12U, cell->n_lines()); ++l)
line_indices[l] = cell->line_index(l);

return line_indices;
}

/**
* For a given Triangulation, update that part of the number
* cache that relates to lines. For 1d, we have to deal with the
Expand Down Expand Up @@ -10754,6 +10813,8 @@ namespace internal
Triangulation<3, spacedim> &triangulation)
{
const unsigned int dim = 3;
using line_iterator =
typename Triangulation<dim, spacedim>::raw_line_iterator;

// first clear flags on lines, since we need them to determine
// which lines will be refined
Expand Down Expand Up @@ -10782,24 +10843,32 @@ namespace internal
for (const auto &cell : triangulation.cell_iterators())
if (cell->refine_flag_set())
{
for (unsigned int line = 0; line < cell->n_lines(); ++line)
const std::array<unsigned int, 12> line_indices =
get_line_indices_of_cell(cell);
for (unsigned int l = 0; l < cell->n_lines(); ++l)
if (GeometryInfo<dim>::line_refinement_case(
cell->refine_flag_set(), line) ==
cell->refine_flag_set(), l) ==
RefinementCase<1>::cut_x)
// flag a line, that will be
// refined
cell->line(line)->set_user_flag();
{
// flag a line to be refined
line_iterator(&triangulation, 0, line_indices[l])
->set_user_flag();
}
}
else if (cell->has_children() &&
!cell->child(0)->coarsen_flag_set())
{
for (unsigned int line = 0; line < cell->n_lines(); ++line)
const std::array<unsigned int, 12> line_indices =
get_line_indices_of_cell(cell);
for (unsigned int l = 0; l < cell->n_lines(); ++l)
if (GeometryInfo<dim>::line_refinement_case(
cell->refinement_case(), line) ==
cell->refinement_case(), l) ==
RefinementCase<1>::cut_x)
// flag a line, that is refined
// and will stay so
cell->line(line)->set_user_flag();
{
// flag a line that is refined and will stay so
line_iterator(&triangulation, 0, line_indices[l])
->set_user_flag();
}
}
else if (cell->has_children() &&
cell->child(0)->coarsen_flag_set())
Expand All @@ -10814,75 +10883,87 @@ namespace internal
cell = triangulation.last_active();
cell != triangulation.end();
--cell)
for (unsigned int line = 0; line < cell->n_lines(); ++line)
{
if (cell->line(line)->has_children())
{
// if this line is refined, its children should
// not have further children
//
// however, if any of the children is flagged
// for further refinement, we need to refine
// this cell also (at least, if the cell is not
// already flagged)
bool offending_line_found = false;
{
const std::array<unsigned int, 12> line_indices =
get_line_indices_of_cell(cell);
for (unsigned int l = 0; l < cell->n_lines(); ++l)
{
line_iterator line(&triangulation, 0, line_indices[l]);
if (line->has_children())
{
// if this line is refined, its children should
// not have further children
//
// however, if any of the children is flagged
// for further refinement, we need to refine
// this cell also (at least, if the cell is not
// already flagged)
bool offending_line_found = false;

for (unsigned int c = 0; c < 2; ++c)
{
Assert(line->child(c)->has_children() == false,
ExcInternalError());

for (unsigned int c = 0; c < 2; ++c)
{
Assert(cell->line(line)->child(c)->has_children() ==
false,
ExcInternalError());

if (cell->line(line)->child(c)->user_flag_set() &&
(GeometryInfo<dim>::line_refinement_case(
cell->refine_flag_set(), line) ==
RefinementCase<1>::no_refinement))
{
// tag this cell for refinement
cell->clear_coarsen_flag();
// if anisotropic coarsening is allowed:
// extend the refine_flag in the needed
// direction, else set refine_flag
// (isotropic)
if (triangulation.smooth_grid &
Triangulation<dim, spacedim>::
allow_anisotropic_smoothing)
cell->flag_for_line_refinement(line);
else
cell->set_refine_flag();

for (unsigned int l = 0; l < cell->n_lines(); ++l)
if (GeometryInfo<dim>::line_refinement_case(
cell->refine_flag_set(), line) ==
RefinementCase<1>::cut_x)
// flag a line, that will be refined
cell->line(l)->set_user_flag();

// note that we have changed the grid
offending_line_found = true;

// it may save us several loop
// iterations if we flag all lines of
// this cell now (and not at the outset
// of the next iteration) for refinement
for (unsigned int l = 0; l < cell->n_lines(); ++l)
if (!cell->line(l)->has_children() &&
(GeometryInfo<dim>::line_refinement_case(
cell->refine_flag_set(), l) !=
RefinementCase<1>::no_refinement))
cell->line(l)->set_user_flag();

break;
}
}
if (line->child(c)->user_flag_set() &&
(GeometryInfo<dim>::line_refinement_case(
cell->refine_flag_set(), l) ==
RefinementCase<1>::no_refinement))
{
// tag this cell for refinement
cell->clear_coarsen_flag();
// if anisotropic coarsening is allowed:
// extend the refine_flag in the needed
// direction, else set refine_flag
// (isotropic)
if (triangulation.smooth_grid &
Triangulation<dim, spacedim>::
allow_anisotropic_smoothing)
cell->flag_for_line_refinement(l);
else
cell->set_refine_flag();

for (unsigned int k = 0; k < cell->n_lines();
++k)
if (GeometryInfo<dim>::line_refinement_case(
cell->refine_flag_set(), k) ==
RefinementCase<1>::cut_x)
// flag a line, that will be refined
line_iterator(&triangulation,
0,
line_indices[k])
->set_user_flag();

// note that we have changed the grid
offending_line_found = true;

// it may save us several loop
// iterations if we flag all lines of
// this cell now (and not at the outset
// of the next iteration) for refinement
for (unsigned int k = 0; k < cell->n_lines();
++k)
if (!cell->line(k)->has_children() &&
(GeometryInfo<dim>::line_refinement_case(
cell->refine_flag_set(), k) !=
RefinementCase<1>::no_refinement))
line_iterator(&triangulation,
0,
line_indices[k])
->set_user_flag();

break;
}
}

if (offending_line_found)
{
mesh_changed = true;
break;
}
}
}
if (offending_line_found)
{
mesh_changed = true;
break;
}
}
}
}


// there is another thing here: if any of the lines will
Expand All @@ -10899,27 +10980,33 @@ namespace internal
triangulation.last();
cell != triangulation.end();
--cell)
{
if (cell->user_flag_set())
for (unsigned int line = 0; line < cell->n_lines(); ++line)
if (cell->line(line)->has_children() &&
(cell->line(line)->child(0)->user_flag_set() ||
cell->line(line)->child(1)->user_flag_set()))
{
for (unsigned int c = 0; c < cell->n_children(); ++c)
cell->child(c)->clear_coarsen_flag();
cell->clear_user_flag();
for (unsigned int l = 0; l < cell->n_lines(); ++l)
if (GeometryInfo<dim>::line_refinement_case(
cell->refinement_case(), l) ==
RefinementCase<1>::cut_x)
// flag a line, that is refined
// and will stay so
cell->line(l)->set_user_flag();
mesh_changed = true;
break;
}
}
if (cell->user_flag_set())
{
const std::array<unsigned int, 12> line_indices =
get_line_indices_of_cell(cell);
for (unsigned int l = 0; l < cell->n_lines(); ++l)
{
line_iterator line(&triangulation, 0, line_indices[l]);
if (line->has_children() &&
(line->child(0)->user_flag_set() ||
line->child(1)->user_flag_set()))
{
for (unsigned int c = 0; c < cell->n_children(); ++c)
cell->child(c)->clear_coarsen_flag();
cell->clear_user_flag();
for (unsigned int k = 0; k < cell->n_lines(); ++k)
if (GeometryInfo<dim>::line_refinement_case(
cell->refinement_case(), k) ==
RefinementCase<1>::cut_x)
// flag a line that is refined and will stay
// so
line_iterator(&triangulation, 0, line_indices[k])
->set_user_flag();
mesh_changed = true;
break;
}
}
}
}
while (mesh_changed == true);
}
Expand Down Expand Up @@ -11463,15 +11550,23 @@ std::vector<types::manifold_id>
Triangulation<dim, spacedim>::get_manifold_ids() const
{
std::set<types::manifold_id> m_ids;
for (auto cell : active_cell_iterators())
for (const auto &cell : active_cell_iterators())
if (cell->is_locally_owned())
{
m_ids.insert(cell->manifold_id());
for (const auto &face : cell->face_iterators())
m_ids.insert(face->manifold_id());
if (dim == 3)
for (const unsigned int l : cell->line_indices())
m_ids.insert(cell->line(l)->manifold_id());
{
const std::array<unsigned int, 12> line_indices =
internal::TriangulationImplementation::Implementation::
get_line_indices_of_cell(cell);
for (unsigned int l = 0; l < cell->n_lines(); ++l)
{
raw_line_iterator line(this, 0, line_indices[l]);
m_ids.insert(line->manifold_id());
}
}
}
return {m_ids.begin(), m_ids.end()};
}
Expand Down Expand Up @@ -11831,15 +11926,22 @@ Triangulation<dim, spacedim>::create_triangulation(
{
while (cell_info->id != cell->id().template to_binary<dim>())
++cell;
if (dim > 1)
for (const auto face : cell->face_indices())
cell->face(face)->set_manifold_id(
cell_info->manifold_quad_ids[face]);

if (dim == 3)
for (const auto quad : cell->face_indices())
cell->quad(quad)->set_manifold_id(
cell_info->manifold_quad_ids[quad]);

if (dim >= 2)
for (const auto line : cell->line_indices())
cell->line(line)->set_manifold_id(
cell_info->manifold_line_ids[line]);
{
const std::array<unsigned int, 12> line_indices =
internal::TriangulationImplementation::Implementation::
get_line_indices_of_cell(cell);
for (unsigned int l = 0; l < cell->n_lines(); ++l)
{
raw_line_iterator line(this, 0, line_indices[l]);
line->set_manifold_id(cell_info->manifold_line_ids[l]);
}
}

cell->set_manifold_id(cell_info->manifold_id);
}
Expand Down

0 comments on commit d4d2cd4

Please sign in to comment.