Reference documentation for deal.II version GIT 35969cdc9b 2023-12-09 01:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 #include <deal.II/base/mpi.h>
18 #include <deal.II/base/mpi.templates.h>
22 
23 #ifdef DEAL_II_WITH_ARBORX
26 #else
27 template <int dim, typename Number>
28 class BoundingBox;
29 
30 namespace ArborXWrappers
31 {
32  class DistributedTree
33  {
34  public:
35  template <int dim, typename Number>
36  DistributedTree(const MPI_Comm &,
37  const std::vector<BoundingBox<dim, Number>> &);
38  template <typename QueryType>
39  std::pair<std::vector<std::pair<int, int>>, std::vector<int>>
40  query(const QueryType &queries);
41  };
42  class BoundingBoxIntersectPredicate
43  {};
44 } // namespace ArborXWrappers
45 #endif
46 
47 #ifdef DEAL_II_WITH_CGAL
49 # include <deal.II/cgal/utilities.h>
50 #endif
51 
56 
59 #include <deal.II/dofs/dof_tools.h>
60 
61 #include <deal.II/fe/fe_nothing.h>
62 #include <deal.II/fe/fe_q.h>
63 #include <deal.II/fe/fe_values.h>
64 #include <deal.II/fe/mapping_fe.h>
65 #include <deal.II/fe/mapping_q.h>
66 
70 #include <deal.II/grid/manifold.h>
71 #include <deal.II/grid/tria.h>
74 
79 #include <deal.II/lac/solver_cg.h>
83 #include <deal.II/lac/vector.h>
85 
88 
90 
91 #include <boost/random/mersenne_twister.hpp>
92 #include <boost/random/uniform_real_distribution.hpp>
93 
94 #include <array>
95 #include <cmath>
96 #include <iostream>
97 #include <limits>
98 #include <list>
99 #include <numeric>
100 #include <set>
101 #include <tuple>
102 #include <unordered_map>
103 
105 
106 
107 namespace GridTools
108 {
109  // define some transformations
110  namespace internal
111  {
112  template <int spacedim>
113  class Shift
114  {
115  public:
116  explicit Shift(const Tensor<1, spacedim> &shift)
117  : shift(shift)
118  {}
121  {
122  return p + shift;
123  }
124 
125  private:
127  };
128 
129 
130  // Transformation to rotate around one of the cartesian z-axis in 2d.
131  class Rotate2d
132  {
133  public:
134  explicit Rotate2d(const double angle)
135  : rotation_matrix(
136  Physics::Transformations::Rotations::rotation_matrix_2d(angle))
137  {}
138  Point<2>
139  operator()(const Point<2> &p) const
140  {
141  return static_cast<Point<2>>(rotation_matrix * p);
142  }
143 
144  private:
146  };
147 
148 
149  // Transformation to rotate around one of the cartesian axes.
150  class Rotate3d
151  {
152  public:
153  Rotate3d(const Tensor<1, 3, double> &axis, const double angle)
154  : rotation_matrix(
155  Physics::Transformations::Rotations::rotation_matrix_3d(axis,
156  angle))
157  {}
158 
159  Point<3>
160  operator()(const Point<3> &p) const
161  {
162  return static_cast<Point<3>>(rotation_matrix * p);
163  }
164 
165  private:
167  };
168 
169 
170  template <int spacedim>
171  class Scale
172  {
173  public:
174  explicit Scale(const double factor)
175  : factor(factor)
176  {}
179  {
180  return p * factor;
181  }
182 
183  private:
184  const double factor;
185  };
186  } // namespace internal
187 
188 
189  template <int dim, int spacedim>
190  void
191  shift(const Tensor<1, spacedim> &shift_vector,
193  {
195  }
196 
197 
198 
199  template <int dim, int spacedim>
200  void
202  {
203  (void)angle;
204  (void)triangulation;
205 
206  AssertThrow(false,
207  ExcMessage(
208  "GridTools::rotate() is only available for spacedim = 2."));
209  }
210 
211 
212 
213  template <>
214  void
216  {
218  }
219 
220 
221 
222  template <>
223  void
225  {
227  }
228 
229 
230  template <int dim>
231  void
233  const double angle,
235  {
237  }
238 
239 
240  template <int dim>
241  void
242  rotate(const double angle,
243  const unsigned int axis,
245  {
246  Assert(axis < 3, ExcMessage("Invalid axis given!"));
247 
248  Tensor<1, 3, double> vector;
249  vector[axis] = 1.;
250 
252  }
253 
254 
255  template <int dim, int spacedim>
256  void
257  scale(const double scaling_factor,
259  {
260  Assert(scaling_factor > 0, ExcScalingFactorNotPositive(scaling_factor));
262  }
263 
264 
265  namespace internal
266  {
272  inline void
274  const AffineConstraints<double> &constraints,
275  Vector<double> &u)
276  {
277  const unsigned int n_dofs = S.n();
278  const auto op = linear_operator(S);
279  const auto SF = constrained_linear_operator(constraints, op);
281  prec.initialize(S, 1.2);
282 
283  SolverControl control(n_dofs, 1.e-10, false, false);
285  SolverCG<Vector<double>> solver(control, mem);
286 
287  Vector<double> f(n_dofs);
288 
289  const auto constrained_rhs =
290  constrained_right_hand_side(constraints, op, f);
291  solver.solve(SF, u, constrained_rhs, prec);
292 
293  constraints.distribute(u);
294  }
295  } // namespace internal
296 
297 
298  // Implementation for dimensions except 1
299  template <int dim>
300  void
301  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
303  const Function<dim> *coefficient,
304  const bool solve_for_absolute_positions)
305  {
306  if (dim == 1)
307  Assert(false, ExcNotImplemented());
308 
309  // first provide everything that is needed for solving a Laplace
310  // equation.
311  FE_Q<dim> q1(1);
312 
313  DoFHandler<dim> dof_handler(triangulation);
314  dof_handler.distribute_dofs(q1);
315 
316  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
317  DoFTools::make_sparsity_pattern(dof_handler, dsp);
318  dsp.compress();
319 
320  SparsityPattern sparsity_pattern;
321  sparsity_pattern.copy_from(dsp);
322  sparsity_pattern.compress();
323 
324  SparseMatrix<double> S(sparsity_pattern);
325 
326  QGauss<dim> quadrature(4);
327 
328  Assert(triangulation.all_reference_cells_are_hyper_cube(),
330  const auto reference_cell = ReferenceCells::get_hypercube<dim>();
332  reference_cell.template get_default_linear_mapping<dim, dim>(),
333  dof_handler,
334  quadrature,
335  S,
336  coefficient);
337 
338  // set up the boundary values for the laplace problem
339  std::array<AffineConstraints<double>, dim> constraints;
340  typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
341  new_points.end();
342 
343  // Fill these maps using the data given by new_points
344  for (const auto &cell : dof_handler.active_cell_iterators())
345  {
346  // Loop over all vertices of the cell and see if it is listed in the map
347  // given as first argument of the function. We visit vertices multiple
348  // times, so also check that if we have already added a constraint, we
349  // don't do it a second time again.
350  for (const unsigned int vertex_no : cell->vertex_indices())
351  {
352  const unsigned int vertex_index = cell->vertex_index(vertex_no);
353  const Point<dim> &vertex_point = cell->vertex(vertex_no);
354 
355  const typename std::map<unsigned int, Point<dim>>::const_iterator
356  map_iter = new_points.find(vertex_index);
357 
358  if (map_iter != map_end)
359  for (unsigned int i = 0; i < dim; ++i)
360  if (constraints[i].is_constrained(
361  cell->vertex_dof_index(vertex_no, 0)) == false)
362  {
363  constraints[i].add_constraint(
364  cell->vertex_dof_index(vertex_no, 0),
365  {},
366  (solve_for_absolute_positions ?
367  map_iter->second(i) :
368  map_iter->second(i) - vertex_point[i]));
369  }
370  }
371  }
372 
373  for (unsigned int i = 0; i < dim; ++i)
374  constraints[i].close();
375 
376  // solve the dim problems with different right hand sides.
377  Vector<double> us[dim];
378  for (unsigned int i = 0; i < dim; ++i)
379  us[i].reinit(dof_handler.n_dofs());
380 
381  // solve linear systems in parallel
382  Threads::TaskGroup<> tasks;
383  for (unsigned int i = 0; i < dim; ++i)
384  tasks +=
385  Threads::new_task(&internal::laplace_solve, S, constraints[i], us[i]);
386  tasks.join_all();
387 
388  // change the coordinates of the points of the triangulation
389  // according to the computed values
390  std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
391  for (const auto &cell : dof_handler.active_cell_iterators())
392  for (const unsigned int vertex_no : cell->vertex_indices())
393  if (vertex_touched[cell->vertex_index(vertex_no)] == false)
394  {
395  Point<dim> &v = cell->vertex(vertex_no);
396 
397  const types::global_dof_index dof_index =
398  cell->vertex_dof_index(vertex_no, 0);
399  for (unsigned int i = 0; i < dim; ++i)
400  if (solve_for_absolute_positions)
401  v(i) = us[i](dof_index);
402  else
403  v(i) += us[i](dof_index);
404 
405  vertex_touched[cell->vertex_index(vertex_no)] = true;
406  }
407  }
408 
413  template <int dim, int spacedim>
414  void
415  distort_random(const double factor,
417  const bool keep_boundary,
418  const unsigned int seed)
419  {
420  // if spacedim>dim we need to make sure that we perturb
421  // points but keep them on
422  // the manifold. however, this isn't implemented right now
423  Assert(spacedim == dim, ExcNotImplemented());
424 
425 
426  // find the smallest length of the
427  // lines adjacent to the
428  // vertex. take the initial value
429  // to be larger than anything that
430  // might be found: the diameter of
431  // the triangulation, here
432  // estimated by adding up the
433  // diameters of the coarse grid
434  // cells.
435  double almost_infinite_length = 0;
437  triangulation.begin(0);
438  cell != triangulation.end(0);
439  ++cell)
440  almost_infinite_length += cell->diameter();
441 
442  std::vector<double> minimal_length(triangulation.n_vertices(),
443  almost_infinite_length);
444 
445  // also note if a vertex is at the boundary
446  std::vector<bool> at_boundary(keep_boundary ? triangulation.n_vertices() :
447  0,
448  false);
449  // for parallel::shared::Triangulation we need to work on all vertices,
450  // not just the ones related to locally owned cells;
451  const bool is_parallel_shared =
453  &triangulation) != nullptr);
454  for (const auto &cell : triangulation.active_cell_iterators())
455  if (is_parallel_shared || cell->is_locally_owned())
456  {
457  if (dim > 1)
458  {
459  for (unsigned int i = 0; i < cell->n_lines(); ++i)
460  {
462  line = cell->line(i);
463 
464  if (keep_boundary && line->at_boundary())
465  {
466  at_boundary[line->vertex_index(0)] = true;
467  at_boundary[line->vertex_index(1)] = true;
468  }
469 
470  minimal_length[line->vertex_index(0)] =
471  std::min(line->diameter(),
472  minimal_length[line->vertex_index(0)]);
473  minimal_length[line->vertex_index(1)] =
474  std::min(line->diameter(),
475  minimal_length[line->vertex_index(1)]);
476  }
477  }
478  else // dim==1
479  {
480  if (keep_boundary)
481  for (unsigned int vertex = 0; vertex < 2; ++vertex)
482  if (cell->at_boundary(vertex) == true)
483  at_boundary[cell->vertex_index(vertex)] = true;
484 
485  minimal_length[cell->vertex_index(0)] =
486  std::min(cell->diameter(),
487  minimal_length[cell->vertex_index(0)]);
488  minimal_length[cell->vertex_index(1)] =
489  std::min(cell->diameter(),
490  minimal_length[cell->vertex_index(1)]);
491  }
492  }
493 
494  // create a random number generator for the interval [-1,1]
495  boost::random::mt19937 rng(seed);
496  boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
497 
498  // If the triangulation is distributed, we need to
499  // exchange the moved vertices across mpi processes
500  if (auto distributed_triangulation =
502  &triangulation))
503  {
504  const std::vector<bool> locally_owned_vertices =
506  std::vector<bool> vertex_moved(triangulation.n_vertices(), false);
507 
508  // Next move vertices on locally owned cells
509  for (const auto &cell : triangulation.active_cell_iterators())
510  if (cell->is_locally_owned())
511  {
512  for (const unsigned int vertex_no : cell->vertex_indices())
513  {
514  const unsigned global_vertex_no =
515  cell->vertex_index(vertex_no);
516 
517  // ignore this vertex if we shall keep the boundary and
518  // this vertex *is* at the boundary, if it is already moved
519  // or if another process moves this vertex
520  if ((keep_boundary && at_boundary[global_vertex_no]) ||
521  vertex_moved[global_vertex_no] ||
522  !locally_owned_vertices[global_vertex_no])
523  continue;
524 
525  // first compute a random shift vector
526  Point<spacedim> shift_vector;
527  for (unsigned int d = 0; d < spacedim; ++d)
528  shift_vector(d) = uniform_distribution(rng);
529 
530  shift_vector *= factor * minimal_length[global_vertex_no] /
531  std::sqrt(shift_vector.square());
532 
533  // finally move the vertex
534  cell->vertex(vertex_no) += shift_vector;
535  vertex_moved[global_vertex_no] = true;
536  }
537  }
538 
539  distributed_triangulation->communicate_locally_moved_vertices(
540  locally_owned_vertices);
541  }
542  else
543  // if this is a sequential triangulation, we could in principle
544  // use the algorithm above, but we'll use an algorithm that we used
545  // before the parallel::distributed::Triangulation was introduced
546  // in order to preserve backward compatibility
547  {
548  // loop over all vertices and compute their new locations
549  const unsigned int n_vertices = triangulation.n_vertices();
550  std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
551  const std::vector<Point<spacedim>> &old_vertex_locations =
552  triangulation.get_vertices();
553 
554  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
555  {
556  // ignore this vertex if we will keep the boundary and
557  // this vertex *is* at the boundary
558  if (keep_boundary && at_boundary[vertex])
559  new_vertex_locations[vertex] = old_vertex_locations[vertex];
560  else
561  {
562  // compute a random shift vector
563  Point<spacedim> shift_vector;
564  for (unsigned int d = 0; d < spacedim; ++d)
565  shift_vector(d) = uniform_distribution(rng);
566 
567  shift_vector *= factor * minimal_length[vertex] /
568  std::sqrt(shift_vector.square());
569 
570  // record new vertex location
571  new_vertex_locations[vertex] =
572  old_vertex_locations[vertex] + shift_vector;
573  }
574  }
575 
576  // now do the actual move of the vertices
577  for (const auto &cell : triangulation.active_cell_iterators())
578  for (const unsigned int vertex_no : cell->vertex_indices())
579  cell->vertex(vertex_no) =
580  new_vertex_locations[cell->vertex_index(vertex_no)];
581  }
582 
583  // Correct hanging nodes if necessary
584  if (dim >= 2)
585  {
586  // We do the same as in GridTools::transform
587  //
588  // exclude hanging nodes at the boundaries of artificial cells:
589  // these may belong to ghost cells for which we know the exact
590  // location of vertices, whereas the artificial cell may or may
591  // not be further refined, and so we cannot know whether
592  // the location of the hanging node is correct or not
594  cell = triangulation.begin_active(),
595  endc = triangulation.end();
596  for (; cell != endc; ++cell)
597  if (!cell->is_artificial())
598  for (const unsigned int face : cell->face_indices())
599  if (cell->face(face)->has_children() &&
600  !cell->face(face)->at_boundary())
601  {
602  // this face has hanging nodes
603  if (dim == 2)
604  cell->face(face)->child(0)->vertex(1) =
605  (cell->face(face)->vertex(0) +
606  cell->face(face)->vertex(1)) /
607  2;
608  else if (dim == 3)
609  {
610  cell->face(face)->child(0)->vertex(1) =
611  .5 * (cell->face(face)->vertex(0) +
612  cell->face(face)->vertex(1));
613  cell->face(face)->child(0)->vertex(2) =
614  .5 * (cell->face(face)->vertex(0) +
615  cell->face(face)->vertex(2));
616  cell->face(face)->child(1)->vertex(3) =
617  .5 * (cell->face(face)->vertex(1) +
618  cell->face(face)->vertex(3));
619  cell->face(face)->child(2)->vertex(3) =
620  .5 * (cell->face(face)->vertex(2) +
621  cell->face(face)->vertex(3));
622 
623  // center of the face
624  cell->face(face)->child(0)->vertex(3) =
625  .25 * (cell->face(face)->vertex(0) +
626  cell->face(face)->vertex(1) +
627  cell->face(face)->vertex(2) +
628  cell->face(face)->vertex(3));
629  }
630  }
631  }
632  }
633 
634 
635 
636  template <int dim, template <int, int> class MeshType, int spacedim>
638  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
639  unsigned int find_closest_vertex(const MeshType<dim, spacedim> &mesh,
640  const Point<spacedim> &p,
641  const std::vector<bool> &marked_vertices)
642  {
643  // first get the underlying triangulation from the mesh and determine
644  // vertices and used vertices
646 
647  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
648 
649  Assert(tria.get_vertices().size() == marked_vertices.size() ||
650  marked_vertices.empty(),
652  marked_vertices.size()));
653 
654  // marked_vertices is expected to be a subset of used_vertices. Thus,
655  // comparing the range marked_vertices.begin() to marked_vertices.end() with
656  // the range used_vertices.begin() to used_vertices.end() the element in the
657  // second range must be valid if the element in the first range is valid.
658  Assert(
659  marked_vertices.empty() ||
660  std::equal(marked_vertices.begin(),
661  marked_vertices.end(),
662  tria.get_used_vertices().begin(),
663  [](bool p, bool q) { return !p || q; }),
664  ExcMessage(
665  "marked_vertices should be a subset of used vertices in the triangulation "
666  "but marked_vertices contains one or more vertices that are not used vertices!"));
667 
668  // If marked_indices is empty, consider all used_vertices for finding the
669  // closest vertex to the point. Otherwise, marked_indices is used.
670  const std::vector<bool> &vertices_to_use =
671  (marked_vertices.empty()) ? tria.get_used_vertices() : marked_vertices;
672 
673  // At the beginning, the first used vertex is considered to be the closest
674  // one.
675  std::vector<bool>::const_iterator first =
676  std::find(vertices_to_use.begin(), vertices_to_use.end(), true);
677 
678  // Assert that at least one vertex is actually used
679  Assert(first != vertices_to_use.end(), ExcInternalError());
680 
681  unsigned int best_vertex = std::distance(vertices_to_use.begin(), first);
682  double best_dist = (p - vertices[best_vertex]).norm_square();
683 
684  // For all remaining vertices, test
685  // whether they are any closer
686  for (unsigned int j = best_vertex + 1; j < vertices.size(); ++j)
687  if (vertices_to_use[j])
688  {
689  const double dist = (p - vertices[j]).norm_square();
690  if (dist < best_dist)
691  {
692  best_vertex = j;
693  best_dist = dist;
694  }
695  }
696 
697  return best_vertex;
698  }
699 
700 
701 
702  template <int dim, template <int, int> class MeshType, int spacedim>
704  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
705  unsigned int find_closest_vertex(const Mapping<dim, spacedim> &mapping,
706  const MeshType<dim, spacedim> &mesh,
707  const Point<spacedim> &p,
708  const std::vector<bool> &marked_vertices)
709  {
710  // Take a shortcut in the simple case.
711  if (mapping.preserves_vertex_locations() == true)
712  return find_closest_vertex(mesh, p, marked_vertices);
713 
714  // first get the underlying triangulation from the mesh and determine
715  // vertices and used vertices
717 
718  auto vertices = extract_used_vertices(tria, mapping);
719 
720  Assert(tria.get_vertices().size() == marked_vertices.size() ||
721  marked_vertices.empty(),
723  marked_vertices.size()));
724 
725  // marked_vertices is expected to be a subset of used_vertices. Thus,
726  // comparing the range marked_vertices.begin() to marked_vertices.end()
727  // with the range used_vertices.begin() to used_vertices.end() the element
728  // in the second range must be valid if the element in the first range is
729  // valid.
730  Assert(
731  marked_vertices.empty() ||
732  std::equal(marked_vertices.begin(),
733  marked_vertices.end(),
734  tria.get_used_vertices().begin(),
735  [](bool p, bool q) { return !p || q; }),
736  ExcMessage(
737  "marked_vertices should be a subset of used vertices in the triangulation "
738  "but marked_vertices contains one or more vertices that are not used vertices!"));
739 
740  // Remove from the map unwanted elements.
741  if (marked_vertices.size() != 0)
742  for (auto it = vertices.begin(); it != vertices.end();)
743  {
744  if (marked_vertices[it->first] == false)
745  {
746  it = vertices.erase(it);
747  }
748  else
749  {
750  ++it;
751  }
752  }
753 
754  return find_closest_vertex(vertices, p);
755  }
756 
757 
758 
759  template <int dim, int spacedim>
760  std::vector<std::vector<Tensor<1, spacedim>>>
762  const Triangulation<dim, spacedim> &mesh,
763  const std::vector<
765  &vertex_to_cells)
766  {
767  const std::vector<Point<spacedim>> &vertices = mesh.get_vertices();
768  const unsigned int n_vertices = vertex_to_cells.size();
769 
770  AssertDimension(vertices.size(), n_vertices);
771 
772 
773  std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
774  n_vertices);
775  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
776  if (mesh.vertex_used(vertex))
777  {
778  const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
779  vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
780 
781  typename std::set<typename Triangulation<dim, spacedim>::
782  active_cell_iterator>::iterator it =
783  vertex_to_cells[vertex].begin();
784  for (unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
785  {
786  vertex_to_cell_centers[vertex][cell] =
787  (*it)->center() - vertices[vertex];
788  vertex_to_cell_centers[vertex][cell] /=
789  vertex_to_cell_centers[vertex][cell].norm();
790  }
791  }
792  return vertex_to_cell_centers;
793  }
794 
795 
796  namespace internal
797  {
798  template <int spacedim>
799  bool
801  const unsigned int a,
802  const unsigned int b,
803  const Tensor<1, spacedim> &point_direction,
804  const std::vector<Tensor<1, spacedim>> &center_directions)
805  {
806  const double scalar_product_a = center_directions[a] * point_direction;
807  const double scalar_product_b = center_directions[b] * point_direction;
808 
809  // The function is supposed to return if a is before b. We are looking
810  // for the alignment of point direction and center direction, therefore
811  // return if the scalar product of a is larger.
812  return (scalar_product_a > scalar_product_b);
813  }
814  } // namespace internal
815 
816  template <int dim, template <int, int> class MeshType, int spacedim>
818  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
819 #ifndef _MSC_VER
820  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
821 #else
822  std::pair<typename ::internal::
823  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
824  Point<dim>>
825 #endif
827  const Mapping<dim, spacedim> &mapping,
828  const MeshType<dim, spacedim> &mesh,
829  const Point<spacedim> &p,
830  const std::vector<
831  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
832  &vertex_to_cells,
833  const std::vector<std::vector<Tensor<1, spacedim>>>
834  &vertex_to_cell_centers,
835  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
836  const std::vector<bool> &marked_vertices,
837  const RTree<std::pair<Point<spacedim>, unsigned int>>
838  &used_vertices_rtree,
839  const double tolerance,
840  const RTree<
841  std::pair<BoundingBox<spacedim>,
843  *relevant_cell_bounding_boxes_rtree)
844  {
845  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
846  Point<dim>>
847  cell_and_position;
848  cell_and_position.first = mesh.end();
849 
850  // To handle points at the border we keep track of points which are close to
851  // the unit cell:
852  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
853  Point<dim>>
854  cell_and_position_approx;
855 
856  if (relevant_cell_bounding_boxes_rtree != nullptr &&
857  !relevant_cell_bounding_boxes_rtree->empty())
858  {
859  // create a bounding box around point p with 2*tolerance as side length.
860  const auto bb = BoundingBox<spacedim>(p).create_extended(tolerance);
861 
862  if (relevant_cell_bounding_boxes_rtree->qbegin(
863  boost::geometry::index::intersects(bb)) ==
864  relevant_cell_bounding_boxes_rtree->qend())
865  return cell_and_position;
866  }
867 
868  bool found_cell = false;
869  bool approx_cell = false;
870 
871  unsigned int closest_vertex_index = 0;
872  // ensure closest vertex index is a marked one, otherwise cell (with vertex
873  // 0) might be found even though it is not marked. This is only relevant if
874  // searching with rtree, using find_closest_vertex already can manage not
875  // finding points
876  if (marked_vertices.size() && !used_vertices_rtree.empty())
877  {
878  const auto itr =
879  std::find(marked_vertices.begin(), marked_vertices.end(), true);
880  Assert(itr != marked_vertices.end(),
881  ::ExcMessage("No vertex has been marked!"));
882  closest_vertex_index = std::distance(marked_vertices.begin(), itr);
883  }
884 
885  Tensor<1, spacedim> vertex_to_point;
886  auto current_cell = cell_hint;
887 
888  // check whether cell has at least one marked vertex
889  const auto cell_marked = [&mesh, &marked_vertices](const auto &cell) {
890  if (marked_vertices.empty())
891  return true;
892 
893  if (cell != mesh.active_cell_iterators().end())
894  for (unsigned int i = 0; i < cell->n_vertices(); ++i)
895  if (marked_vertices[cell->vertex_index(i)])
896  return true;
897 
898  return false;
899  };
900 
901  // check whether any cell in collection is marked
902  const auto any_cell_marked = [&cell_marked](const auto &cells) {
903  return std::any_of(cells.begin(),
904  cells.end(),
905  [&cell_marked](const auto &cell) {
906  return cell_marked(cell);
907  });
908  };
909  (void)any_cell_marked;
910 
911  while (found_cell == false)
912  {
913  // First look at the vertices of the cell cell_hint. If it's an
914  // invalid cell, then query for the closest global vertex
915  if (current_cell.state() == IteratorState::valid &&
916  cell_marked(cell_hint))
917  {
918  const auto cell_vertices = mapping.get_vertices(current_cell);
919  const unsigned int closest_vertex =
920  find_closest_vertex_of_cell<dim, spacedim>(current_cell,
921  p,
922  mapping);
923  vertex_to_point = p - cell_vertices[closest_vertex];
924  closest_vertex_index = current_cell->vertex_index(closest_vertex);
925  }
926  else
927  {
928  // For some clang-based compilers and boost versions the call to
929  // RTree::query doesn't compile. Since using an rtree here is just a
930  // performance improvement disabling this branch is OK.
931  // This is fixed in boost in
932  // https://github.com/boostorg/numeric_conversion/commit/50a1eae942effb0a9b90724323ef8f2a67e7984a
933 #if defined(DEAL_II_WITH_BOOST_BUNDLED) || \
934  !(defined(__clang_major__) && __clang_major__ >= 16) || \
935  BOOST_VERSION >= 108100
936  if (!used_vertices_rtree.empty())
937  {
938  // If we have an rtree at our disposal, use it.
939  using ValueType = std::pair<Point<spacedim>, unsigned int>;
940  std::function<bool(const ValueType &)> marked;
941  if (marked_vertices.size() == mesh.n_vertices())
942  marked = [&marked_vertices](const ValueType &value) -> bool {
943  return marked_vertices[value.second];
944  };
945  else
946  marked = [](const ValueType &) -> bool { return true; };
947 
948  std::vector<std::pair<Point<spacedim>, unsigned int>> res;
949  used_vertices_rtree.query(
950  boost::geometry::index::nearest(p, 1) &&
951  boost::geometry::index::satisfies(marked),
952  std::back_inserter(res));
953 
954  // Searching for a point which is located outside the
955  // triangulation results in res.size() = 0
956  Assert(res.size() < 2,
957  ::ExcMessage("There can not be multiple results"));
958 
959  if (res.size() > 0)
960  if (any_cell_marked(vertex_to_cells[res[0].second]))
961  closest_vertex_index = res[0].second;
962  }
963  else
964 #endif
965  {
966  closest_vertex_index = GridTools::find_closest_vertex(
967  mapping, mesh, p, marked_vertices);
968  }
969  vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
970  }
971 
972 #ifdef DEBUG
973  {
974  // Double-check if found index is at marked cell
975  Assert(any_cell_marked(vertex_to_cells[closest_vertex_index]),
976  ::ExcMessage("Found non-marked vertex"));
977  }
978 #endif
979 
980  const double vertex_point_norm = vertex_to_point.norm();
981  if (vertex_point_norm > 0)
982  vertex_to_point /= vertex_point_norm;
983 
984  const unsigned int n_neighbor_cells =
985  vertex_to_cells[closest_vertex_index].size();
986 
987  // Create a corresponding map of vectors from vertex to cell center
988  std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
989 
990  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
991  neighbor_permutation[i] = i;
992 
993  auto comp = [&](const unsigned int a, const unsigned int b) -> bool {
994  return internal::compare_point_association<spacedim>(
995  a,
996  b,
997  vertex_to_point,
998  vertex_to_cell_centers[closest_vertex_index]);
999  };
1000 
1001  std::sort(neighbor_permutation.begin(),
1002  neighbor_permutation.end(),
1003  comp);
1004  // It is possible the vertex is close
1005  // to an edge, thus we add a tolerance
1006  // to keep also the "best" cell
1007  double best_distance = tolerance;
1008 
1009  // Search all of the cells adjacent to the closest vertex of the cell
1010  // hint. Most likely we will find the point in them.
1011  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
1012  {
1013  try
1014  {
1015  auto cell = vertex_to_cells[closest_vertex_index].begin();
1016  std::advance(cell, neighbor_permutation[i]);
1017 
1018  if (!(*cell)->is_artificial())
1019  {
1020  const Point<dim> p_unit =
1021  mapping.transform_real_to_unit_cell(*cell, p);
1022  if ((*cell)->reference_cell().contains_point(p_unit,
1023  tolerance))
1024  {
1025  cell_and_position.first = *cell;
1026  cell_and_position.second = p_unit;
1027  found_cell = true;
1028  approx_cell = false;
1029  break;
1030  }
1031  // The point is not inside this cell: checking how far
1032  // outside it is and whether we want to use this cell as a
1033  // backup if we can't find a cell within which the point
1034  // lies.
1035  const double dist = p_unit.distance(
1036  (*cell)->reference_cell().closest_point(p_unit));
1037  if (dist < best_distance)
1038  {
1039  best_distance = dist;
1040  cell_and_position_approx.first = *cell;
1041  cell_and_position_approx.second = p_unit;
1042  approx_cell = true;
1043  }
1044  }
1045  }
1046  catch (typename Mapping<dim>::ExcTransformationFailed &)
1047  {}
1048  }
1049 
1050  if (found_cell == true)
1051  return cell_and_position;
1052  else if (approx_cell == true)
1053  return cell_and_position_approx;
1054 
1055  // The first time around, we check for vertices in the hint_cell. If
1056  // that does not work, we set the cell iterator to an invalid one, and
1057  // look for a global vertex close to the point. If that does not work,
1058  // we are in trouble, and just throw an exception.
1059  //
1060  // If we got here, then we did not find the point. If the
1061  // current_cell.state() here is not IteratorState::valid, it means that
1062  // the user did not provide a hint_cell, and at the beginning of the
1063  // while loop we performed an actual global search on the mesh
1064  // vertices. Not finding the point then means the point is outside the
1065  // domain, or that we've had problems with the algorithm above. Try as a
1066  // last resort the other (simpler) algorithm.
1067  if (current_cell.state() != IteratorState::valid)
1069  mapping, mesh, p, marked_vertices, tolerance);
1070 
1071  current_cell = typename MeshType<dim, spacedim>::active_cell_iterator();
1072  }
1073  return cell_and_position;
1074  }
1075 
1076 
1077 
1078  template <int dim, int spacedim>
1079  unsigned int
1082  const Point<spacedim> &position,
1083  const Mapping<dim, spacedim> &mapping)
1084  {
1085  const auto vertices = mapping.get_vertices(cell);
1086  double minimum_distance = position.distance_square(vertices[0]);
1087  unsigned int closest_vertex = 0;
1088 
1089  for (unsigned int v = 1; v < cell->n_vertices(); ++v)
1090  {
1091  const double vertex_distance = position.distance_square(vertices[v]);
1092  if (vertex_distance < minimum_distance)
1093  {
1094  closest_vertex = v;
1095  minimum_distance = vertex_distance;
1096  }
1097  }
1098  return closest_vertex;
1099  }
1100 
1101 
1102 
1103  namespace internal
1104  {
1105  namespace BoundingBoxPredicate
1106  {
1107  template <typename MeshType>
1109  concepts::is_triangulation_or_dof_handler<MeshType>)
1110  std::tuple<
1112  bool> compute_cell_predicate_bounding_box(const typename MeshType::
1113  cell_iterator &parent_cell,
1114  const std::function<bool(
1115  const typename MeshType::
1116  active_cell_iterator &)>
1117  &predicate)
1118  {
1119  bool has_predicate =
1120  false; // Start assuming there's no cells with predicate inside
1121  std::vector<typename MeshType::active_cell_iterator> active_cells;
1122  if (parent_cell->is_active())
1123  active_cells = {parent_cell};
1124  else
1125  // Finding all active cells descendants of the current one (or the
1126  // current one if it is active)
1127  active_cells = get_active_child_cells<MeshType>(parent_cell);
1128 
1129  const unsigned int spacedim = MeshType::space_dimension;
1130 
1131  // Looking for the first active cell which has the property predicate
1132  unsigned int i = 0;
1133  while (i < active_cells.size() && !predicate(active_cells[i]))
1134  ++i;
1135 
1136  // No active cells or no active cells with property
1137  if (active_cells.empty() || i == active_cells.size())
1138  {
1139  BoundingBox<spacedim> bbox;
1140  return std::make_tuple(bbox, has_predicate);
1141  }
1142 
1143  // The two boundary points defining the boundary box
1144  Point<spacedim> maxp = active_cells[i]->vertex(0);
1145  Point<spacedim> minp = active_cells[i]->vertex(0);
1146 
1147  for (; i < active_cells.size(); ++i)
1148  if (predicate(active_cells[i]))
1149  for (const unsigned int v : active_cells[i]->vertex_indices())
1150  for (unsigned int d = 0; d < spacedim; ++d)
1151  {
1152  minp[d] = std::min(minp[d], active_cells[i]->vertex(v)[d]);
1153  maxp[d] = std::max(maxp[d], active_cells[i]->vertex(v)[d]);
1154  }
1155 
1156  has_predicate = true;
1157  BoundingBox<spacedim> bbox(std::make_pair(minp, maxp));
1158  return std::make_tuple(bbox, has_predicate);
1159  }
1160  } // namespace BoundingBoxPredicate
1161  } // namespace internal
1162 
1163 
1164 
1165  template <typename MeshType>
1166  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1167  std::
1168  vector<BoundingBox<MeshType::space_dimension>> compute_mesh_predicate_bounding_box(
1169  const MeshType &mesh,
1170  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1171  &predicate,
1172  const unsigned int refinement_level,
1173  const bool allow_merge,
1174  const unsigned int max_boxes)
1175  {
1176  // Algorithm brief description: begin with creating bounding boxes of all
1177  // cells at refinement_level (and coarser levels if there are active cells)
1178  // which have the predicate property. These are then merged
1179 
1180  Assert(
1181  refinement_level <= mesh.n_levels(),
1182  ExcMessage(
1183  "Error: refinement level is higher then total levels in the triangulation!"));
1184 
1185  const unsigned int spacedim = MeshType::space_dimension;
1186  std::vector<BoundingBox<spacedim>> bounding_boxes;
1187 
1188  // Creating a bounding box for all active cell on coarser level
1189 
1190  for (unsigned int i = 0; i < refinement_level; ++i)
1191  for (const typename MeshType::cell_iterator &cell :
1192  mesh.active_cell_iterators_on_level(i))
1193  {
1194  bool has_predicate = false;
1195  BoundingBox<spacedim> bbox;
1196  std::tie(bbox, has_predicate) =
1198  MeshType>(cell, predicate);
1199  if (has_predicate)
1200  bounding_boxes.push_back(bbox);
1201  }
1202 
1203  // Creating a Bounding Box for all cells on the chosen refinement_level
1204  for (const typename MeshType::cell_iterator &cell :
1205  mesh.cell_iterators_on_level(refinement_level))
1206  {
1207  bool has_predicate = false;
1208  BoundingBox<spacedim> bbox;
1209  std::tie(bbox, has_predicate) =
1211  MeshType>(cell, predicate);
1212  if (has_predicate)
1213  bounding_boxes.push_back(bbox);
1214  }
1215 
1216  if (!allow_merge)
1217  // If merging is not requested return the created bounding_boxes
1218  return bounding_boxes;
1219  else
1220  {
1221  // Merging part of the algorithm
1222  // Part 1: merging neighbors
1223  // This array stores the indices of arrays we have already merged
1224  std::vector<unsigned int> merged_boxes_idx;
1225  bool found_neighbors = true;
1226 
1227  // We merge only neighbors which can be expressed by a single bounding
1228  // box e.g. in 1d [0,1] and [1,2] can be described with [0,2] without
1229  // losing anything
1230  while (found_neighbors)
1231  {
1232  found_neighbors = false;
1233  for (unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
1234  {
1235  if (std::find(merged_boxes_idx.begin(),
1236  merged_boxes_idx.end(),
1237  i) == merged_boxes_idx.end())
1238  for (unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
1239  if (std::find(merged_boxes_idx.begin(),
1240  merged_boxes_idx.end(),
1241  j) == merged_boxes_idx.end() &&
1242  bounding_boxes[i].get_neighbor_type(
1243  bounding_boxes[j]) ==
1245  {
1246  bounding_boxes[i].merge_with(bounding_boxes[j]);
1247  merged_boxes_idx.push_back(j);
1248  found_neighbors = true;
1249  }
1250  }
1251  }
1252 
1253  // Copying the merged boxes into merged_b_boxes
1254  std::vector<BoundingBox<spacedim>> merged_b_boxes;
1255  for (unsigned int i = 0; i < bounding_boxes.size(); ++i)
1256  if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
1257  merged_boxes_idx.end())
1258  merged_b_boxes.push_back(bounding_boxes[i]);
1259 
1260  // Part 2: if there are too many bounding boxes, merging smaller boxes
1261  // This has sense only in dimension 2 or greater, since in dimension 1,
1262  // neighboring intervals can always be merged without problems
1263  if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
1264  {
1265  std::vector<double> volumes;
1266  for (unsigned int i = 0; i < merged_b_boxes.size(); ++i)
1267  volumes.push_back(merged_b_boxes[i].volume());
1268 
1269  while (merged_b_boxes.size() > max_boxes)
1270  {
1271  unsigned int min_idx =
1272  std::min_element(volumes.begin(), volumes.end()) -
1273  volumes.begin();
1274  volumes.erase(volumes.begin() + min_idx);
1275  // Finding a neighbor
1276  bool not_removed = true;
1277  for (unsigned int i = 0;
1278  i < merged_b_boxes.size() && not_removed;
1279  ++i)
1280  // We merge boxes if we have "attached" or "mergeable"
1281  // neighbors, even though mergeable should be dealt with in
1282  // Part 1
1283  if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
1284  merged_b_boxes[min_idx]) ==
1286  merged_b_boxes[i].get_neighbor_type(
1287  merged_b_boxes[min_idx]) ==
1289  {
1290  merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
1291  merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
1292  not_removed = false;
1293  }
1294  Assert(!not_removed,
1295  ExcMessage("Error: couldn't merge bounding boxes!"));
1296  }
1297  }
1298  Assert(merged_b_boxes.size() <= max_boxes,
1299  ExcMessage(
1300  "Error: couldn't reach target number of bounding boxes!"));
1301  return merged_b_boxes;
1302  }
1303  }
1304 
1305 
1306 
1307  template <int spacedim>
1308 #ifndef DOXYGEN
1309  std::tuple<std::vector<std::vector<unsigned int>>,
1310  std::map<unsigned int, unsigned int>,
1311  std::map<unsigned int, std::vector<unsigned int>>>
1312 #else
1313  return_type
1314 #endif
1316  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1317  const std::vector<Point<spacedim>> &points)
1318  {
1319  unsigned int n_procs = global_bboxes.size();
1320  std::vector<std::vector<unsigned int>> point_owners(n_procs);
1321  std::map<unsigned int, unsigned int> map_owners_found;
1322  std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
1323 
1324  unsigned int n_points = points.size();
1325  for (unsigned int pt = 0; pt < n_points; ++pt)
1326  {
1327  // Keep track of how many processes we guess to own the point
1328  std::vector<unsigned int> owners_found;
1329  // Check in which other processes the point might be
1330  for (unsigned int rk = 0; rk < n_procs; ++rk)
1331  {
1332  for (const BoundingBox<spacedim> &bbox : global_bboxes[rk])
1333  if (bbox.point_inside(points[pt]))
1334  {
1335  point_owners[rk].emplace_back(pt);
1336  owners_found.emplace_back(rk);
1337  break; // We can check now the next process
1338  }
1339  }
1340  Assert(owners_found.size() > 0,
1341  ExcMessage("No owners found for the point " +
1342  std::to_string(pt)));
1343  if (owners_found.size() == 1)
1344  map_owners_found[pt] = owners_found[0];
1345  else
1346  // Multiple owners
1347  map_owners_guessed[pt] = owners_found;
1348  }
1349 
1350  return std::make_tuple(std::move(point_owners),
1351  std::move(map_owners_found),
1352  std::move(map_owners_guessed));
1353  }
1354 
1355  template <int spacedim>
1356 #ifndef DOXYGEN
1357  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1358  std::map<unsigned int, unsigned int>,
1359  std::map<unsigned int, std::vector<unsigned int>>>
1360 #else
1361  return_type
1362 #endif
1364  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1365  const std::vector<Point<spacedim>> &points)
1366  {
1367  std::map<unsigned int, std::vector<unsigned int>> point_owners;
1368  std::map<unsigned int, unsigned int> map_owners_found;
1369  std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
1370  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> search_result;
1371 
1372  unsigned int n_points = points.size();
1373  for (unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
1374  {
1375  search_result.clear(); // clearing last output
1376 
1377  // Running tree search
1378  covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
1379  std::back_inserter(search_result));
1380 
1381  // Keep track of how many processes we guess to own the point
1382  std::set<unsigned int> owners_found;
1383  // Check in which other processes the point might be
1384  for (const auto &rank_bbox : search_result)
1385  {
1386  // Try to add the owner to the owners found,
1387  // and check if it was already present
1388  const bool pt_inserted = owners_found.insert(pt_n).second;
1389  if (pt_inserted)
1390  point_owners[rank_bbox.second].emplace_back(pt_n);
1391  }
1392  Assert(owners_found.size() > 0,
1393  ExcMessage("No owners found for the point " +
1394  std::to_string(pt_n)));
1395  if (owners_found.size() == 1)
1396  map_owners_found[pt_n] = *owners_found.begin();
1397  else
1398  // Multiple owners
1399  std::copy(owners_found.begin(),
1400  owners_found.end(),
1401  std::back_inserter(map_owners_guessed[pt_n]));
1402  }
1403 
1404  return std::make_tuple(std::move(point_owners),
1405  std::move(map_owners_found),
1406  std::move(map_owners_guessed));
1407  }
1408 
1409 
1410 
1411  template <int dim, int spacedim>
1412  std::map<unsigned int, types::global_vertex_index>
1415  {
1416  std::map<unsigned int, types::global_vertex_index>
1417  local_to_global_vertex_index;
1418 
1419 #ifndef DEAL_II_WITH_MPI
1420 
1421  // without MPI, this function doesn't make sense because on cannot
1422  // use parallel::distributed::Triangulation in any meaningful
1423  // way
1424  (void)triangulation;
1425  Assert(false,
1426  ExcMessage("This function does not make any sense "
1427  "for parallel::distributed::Triangulation "
1428  "objects if you do not have MPI enabled."));
1429 
1430 #else
1431 
1432  using active_cell_iterator =
1434  const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
1436 
1437  // Create a local index for the locally "owned" vertices
1438  types::global_vertex_index next_index = 0;
1439  unsigned int max_cellid_size = 0;
1440  std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
1441  vertices_added;
1442  std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
1443  std::map<types::subdomain_id,
1444  std::vector<std::tuple<types::global_vertex_index,
1446  std::string>>>
1447  vertices_to_send;
1448  std::set<active_cell_iterator> missing_vert_cells;
1449  std::set<unsigned int> used_vertex_index;
1450  for (const auto &cell : triangulation.active_cell_iterators())
1451  {
1452  if (cell->is_locally_owned())
1453  {
1454  for (const unsigned int i : cell->vertex_indices())
1455  {
1456  types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
1457  for (const auto &adjacent_cell :
1458  vertex_to_cell[cell->vertex_index(i)])
1459  lowest_subdomain_id = std::min(lowest_subdomain_id,
1460  adjacent_cell->subdomain_id());
1461 
1462  // See if this process "owns" this vertex
1463  if (lowest_subdomain_id == cell->subdomain_id())
1464  {
1465  // Check that the vertex we are working on is a vertex that
1466  // has not been dealt with yet
1467  if (used_vertex_index.find(cell->vertex_index(i)) ==
1468  used_vertex_index.end())
1469  {
1470  // Set the local index
1471  local_to_global_vertex_index[cell->vertex_index(i)] =
1472  next_index++;
1473 
1474  // Store the information that will be sent to the
1475  // adjacent cells on other subdomains
1476  for (const auto &adjacent_cell :
1477  vertex_to_cell[cell->vertex_index(i)])
1478  if (adjacent_cell->subdomain_id() !=
1479  cell->subdomain_id())
1480  {
1481  std::pair<types::subdomain_id,
1483  tmp(adjacent_cell->subdomain_id(),
1484  cell->vertex_index(i));
1485  if (vertices_added.find(tmp) ==
1486  vertices_added.end())
1487  {
1488  vertices_to_send[adjacent_cell
1489  ->subdomain_id()]
1490  .emplace_back(i,
1491  cell->vertex_index(i),
1492  cell->id().to_string());
1493  if (cell->id().to_string().size() >
1494  max_cellid_size)
1495  max_cellid_size =
1496  cell->id().to_string().size();
1497  vertices_added.insert(tmp);
1498  }
1499  }
1500  used_vertex_index.insert(cell->vertex_index(i));
1501  }
1502  }
1503  else
1504  {
1505  // We don't own the vertex so we will receive its global
1506  // index
1507  vertices_to_recv[lowest_subdomain_id].insert(
1508  cell->vertex_index(i));
1509  missing_vert_cells.insert(cell);
1510  }
1511  }
1512  }
1513 
1514  // Some hanging nodes are vertices of ghost cells. They need to be
1515  // received.
1516  if (cell->is_ghost())
1517  {
1518  for (const unsigned int i : cell->face_indices())
1519  {
1520  if (cell->at_boundary(i) == false)
1521  {
1522  if (cell->neighbor(i)->is_active())
1523  {
1524  typename Triangulation<dim,
1525  spacedim>::active_cell_iterator
1526  adjacent_cell = cell->neighbor(i);
1527  if ((adjacent_cell->is_locally_owned()))
1528  {
1529  types::subdomain_id adj_subdomain_id =
1530  adjacent_cell->subdomain_id();
1531  if (cell->subdomain_id() < adj_subdomain_id)
1532  for (unsigned int j = 0;
1533  j < cell->face(i)->n_vertices();
1534  ++j)
1535  {
1536  vertices_to_recv[cell->subdomain_id()].insert(
1537  cell->face(i)->vertex_index(j));
1538  missing_vert_cells.insert(cell);
1539  }
1540  }
1541  }
1542  }
1543  }
1544  }
1545  }
1546 
1547  // Get the size of the largest CellID string
1548  max_cellid_size =
1549  Utilities::MPI::max(max_cellid_size, triangulation.get_communicator());
1550 
1551  // Make indices global by getting the number of vertices owned by each
1552  // processors and shifting the indices accordingly
1554  int ierr = MPI_Exscan(&next_index,
1555  &shift,
1556  1,
1558  MPI_SUM,
1559  triangulation.get_communicator());
1560  AssertThrowMPI(ierr);
1561 
1562  for (auto &global_index_it : local_to_global_vertex_index)
1563  global_index_it.second += shift;
1564 
1565 
1566  const int mpi_tag = Utilities::MPI::internal::Tags::
1568  const int mpi_tag2 = Utilities::MPI::internal::Tags::
1570 
1571 
1572  // In a first message, send the global ID of the vertices and the local
1573  // positions in the cells. In a second messages, send the cell ID as a
1574  // resize string. This is done in two messages so that types are not mixed
1575 
1576  // Send the first message
1577  std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
1578  vertices_to_send.size());
1579  std::vector<MPI_Request> first_requests(vertices_to_send.size());
1580  typename std::map<types::subdomain_id,
1581  std::vector<std::tuple<types::global_vertex_index,
1583  std::string>>>::iterator
1584  vert_to_send_it = vertices_to_send.begin(),
1585  vert_to_send_end = vertices_to_send.end();
1586  for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
1587  ++vert_to_send_it, ++i)
1588  {
1589  int destination = vert_to_send_it->first;
1590  const unsigned int n_vertices = vert_to_send_it->second.size();
1591  const int buffer_size = 2 * n_vertices;
1592  vertices_send_buffers[i].resize(buffer_size);
1593 
1594  // fill the buffer
1595  for (unsigned int j = 0; j < n_vertices; ++j)
1596  {
1597  vertices_send_buffers[i][2 * j] =
1598  std::get<0>(vert_to_send_it->second[j]);
1599  vertices_send_buffers[i][2 * j + 1] =
1600  local_to_global_vertex_index[std::get<1>(
1601  vert_to_send_it->second[j])];
1602  }
1603 
1604  // Send the message
1605  ierr = MPI_Isend(vertices_send_buffers[i].data(),
1606  buffer_size,
1608  destination,
1609  mpi_tag,
1610  triangulation.get_communicator(),
1611  &first_requests[i]);
1612  AssertThrowMPI(ierr);
1613  }
1614 
1615  // Receive the first message
1616  std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
1617  vertices_to_recv.size());
1618  typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
1619  vert_to_recv_it = vertices_to_recv.begin(),
1620  vert_to_recv_end = vertices_to_recv.end();
1621  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1622  ++vert_to_recv_it, ++i)
1623  {
1624  int source = vert_to_recv_it->first;
1625  const unsigned int n_vertices = vert_to_recv_it->second.size();
1626  const int buffer_size = 2 * n_vertices;
1627  vertices_recv_buffers[i].resize(buffer_size);
1628 
1629  // Receive the message
1630  ierr = MPI_Recv(vertices_recv_buffers[i].data(),
1631  buffer_size,
1633  source,
1634  mpi_tag,
1635  triangulation.get_communicator(),
1636  MPI_STATUS_IGNORE);
1637  AssertThrowMPI(ierr);
1638  }
1639 
1640  // At this point, wait for all of the isend operations to finish:
1641  MPI_Waitall(first_requests.size(),
1642  first_requests.data(),
1643  MPI_STATUSES_IGNORE);
1644 
1645 
1646  // Send second message
1647  std::vector<std::vector<char>> cellids_send_buffers(
1648  vertices_to_send.size());
1649  std::vector<MPI_Request> second_requests(vertices_to_send.size());
1650  vert_to_send_it = vertices_to_send.begin();
1651  for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
1652  ++vert_to_send_it, ++i)
1653  {
1654  int destination = vert_to_send_it->first;
1655  const unsigned int n_vertices = vert_to_send_it->second.size();
1656  const int buffer_size = max_cellid_size * n_vertices;
1657  cellids_send_buffers[i].resize(buffer_size);
1658 
1659  // fill the buffer
1660  unsigned int pos = 0;
1661  for (unsigned int j = 0; j < n_vertices; ++j)
1662  {
1663  std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
1664  for (unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
1665  {
1666  if (k < cell_id.size())
1667  cellids_send_buffers[i][pos] = cell_id[k];
1668  // if necessary fill up the reserved part of the buffer with an
1669  // invalid value
1670  else
1671  cellids_send_buffers[i][pos] = '-';
1672  }
1673  }
1674 
1675  // Send the message
1676  ierr = MPI_Isend(cellids_send_buffers[i].data(),
1677  buffer_size,
1678  MPI_CHAR,
1679  destination,
1680  mpi_tag2,
1681  triangulation.get_communicator(),
1682  &second_requests[i]);
1683  AssertThrowMPI(ierr);
1684  }
1685 
1686  // Receive the second message
1687  std::vector<std::vector<char>> cellids_recv_buffers(
1688  vertices_to_recv.size());
1689  vert_to_recv_it = vertices_to_recv.begin();
1690  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1691  ++vert_to_recv_it, ++i)
1692  {
1693  int source = vert_to_recv_it->first;
1694  const unsigned int n_vertices = vert_to_recv_it->second.size();
1695  const int buffer_size = max_cellid_size * n_vertices;
1696  cellids_recv_buffers[i].resize(buffer_size);
1697 
1698  // Receive the message
1699  ierr = MPI_Recv(cellids_recv_buffers[i].data(),
1700  buffer_size,
1701  MPI_CHAR,
1702  source,
1703  mpi_tag2,
1704  triangulation.get_communicator(),
1705  MPI_STATUS_IGNORE);
1706  AssertThrowMPI(ierr);
1707  }
1708 
1709 
1710  // Match the data received with the required vertices
1711  vert_to_recv_it = vertices_to_recv.begin();
1712  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1713  ++i, ++vert_to_recv_it)
1714  {
1715  for (unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
1716  {
1717  const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
1718  const types::global_vertex_index global_id_recv =
1719  vertices_recv_buffers[i][2 * j + 1];
1720  const std::string cellid_recv(
1721  &cellids_recv_buffers[i][max_cellid_size * j],
1722  &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
1723  bool found = false;
1724  typename std::set<active_cell_iterator>::iterator
1725  cell_set_it = missing_vert_cells.begin(),
1726  end_cell_set = missing_vert_cells.end();
1727  for (; (found == false) && (cell_set_it != end_cell_set);
1728  ++cell_set_it)
1729  {
1730  typename std::set<active_cell_iterator>::iterator
1731  candidate_cell =
1732  vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
1733  end_cell =
1734  vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
1735  for (; candidate_cell != end_cell; ++candidate_cell)
1736  {
1737  std::string current_cellid =
1738  (*candidate_cell)->id().to_string();
1739  current_cellid.resize(max_cellid_size, '-');
1740  if (current_cellid.compare(cellid_recv) == 0)
1741  {
1742  local_to_global_vertex_index
1743  [(*candidate_cell)->vertex_index(local_pos_recv)] =
1744  global_id_recv;
1745  found = true;
1746 
1747  break;
1748  }
1749  }
1750  }
1751  }
1752  }
1753 
1754  // At this point, wait for all of the isend operations of the second round
1755  // to finish:
1756  MPI_Waitall(second_requests.size(),
1757  second_requests.data(),
1758  MPI_STATUSES_IGNORE);
1759 #endif
1760 
1761  return local_to_global_vertex_index;
1762  }
1763 
1764 
1765 
1766  template <int dim, int spacedim>
1767  void
1768  partition_triangulation(const unsigned int n_partitions,
1770  const SparsityTools::Partitioner partitioner)
1771  {
1773  &triangulation) == nullptr),
1774  ExcMessage("Objects of type parallel::distributed::Triangulation "
1775  "are already partitioned implicitly and can not be "
1776  "partitioned again explicitly."));
1777 
1778  std::vector<unsigned int> cell_weights;
1779 
1780  // Get cell weighting if a signal has been attached to the triangulation
1781  if (!triangulation.signals.weight.empty())
1782  {
1783  cell_weights.resize(triangulation.n_active_cells(), 0U);
1784 
1785  // In a first step, obtain the weights of the locally owned
1786  // cells. For all others, the weight remains at the zero the
1787  // vector was initialized with above.
1788  for (const auto &cell : triangulation.active_cell_iterators())
1789  if (cell->is_locally_owned())
1790  cell_weights[cell->active_cell_index()] =
1791  triangulation.signals.weight(cell, CellStatus::cell_will_persist);
1792 
1793  // If this is a parallel triangulation, we then need to also
1794  // get the weights for all other cells. We have asserted above
1795  // that this function can't be used for
1796  // parallel::distributed::Triangulation objects, so the only
1797  // ones we have to worry about here are
1798  // parallel::shared::Triangulation
1799  if (const auto shared_tria =
1801  &triangulation))
1802  Utilities::MPI::sum(cell_weights,
1803  shared_tria->get_communicator(),
1804  cell_weights);
1805 
1806  // verify that the global sum of weights is larger than 0
1807  Assert(std::accumulate(cell_weights.begin(),
1808  cell_weights.end(),
1809  std::uint64_t(0)) > 0,
1810  ExcMessage("The global sum of weights over all active cells "
1811  "is zero. Please verify how you generate weights."));
1812  }
1813 
1814  // Call the other more general function
1815  partition_triangulation(n_partitions,
1816  cell_weights,
1817  triangulation,
1818  partitioner);
1819  }
1820 
1821 
1822 
1823  template <int dim, int spacedim>
1824  void
1825  partition_triangulation(const unsigned int n_partitions,
1826  const std::vector<unsigned int> &cell_weights,
1828  const SparsityTools::Partitioner partitioner)
1829  {
1831  &triangulation) == nullptr),
1832  ExcMessage("Objects of type parallel::distributed::Triangulation "
1833  "are already partitioned implicitly and can not be "
1834  "partitioned again explicitly."));
1835  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
1836 
1837  // check for an easy return
1838  if (n_partitions == 1)
1839  {
1840  for (const auto &cell : triangulation.active_cell_iterators())
1841  cell->set_subdomain_id(0);
1842  return;
1843  }
1844 
1845  // we decompose the domain by first
1846  // generating the connection graph of all
1847  // cells with their neighbors, and then
1848  // passing this graph off to METIS.
1849  // finally defer to the other function for
1850  // partitioning and assigning subdomain ids
1851  DynamicSparsityPattern cell_connectivity;
1852  get_face_connectivity_of_cells(triangulation, cell_connectivity);
1853 
1854  SparsityPattern sp_cell_connectivity;
1855  sp_cell_connectivity.copy_from(cell_connectivity);
1856  partition_triangulation(n_partitions,
1857  cell_weights,
1858  sp_cell_connectivity,
1859  triangulation,
1860  partitioner);
1861  }
1862 
1863 
1864 
1865  template <int dim, int spacedim>
1866  void
1867  partition_triangulation(const unsigned int n_partitions,
1868  const SparsityPattern &cell_connection_graph,
1870  const SparsityTools::Partitioner partitioner)
1871  {
1873  &triangulation) == nullptr),
1874  ExcMessage("Objects of type parallel::distributed::Triangulation "
1875  "are already partitioned implicitly and can not be "
1876  "partitioned again explicitly."));
1877 
1878  std::vector<unsigned int> cell_weights;
1879 
1880  // Get cell weighting if a signal has been attached to the triangulation
1881  if (!triangulation.signals.weight.empty())
1882  {
1883  cell_weights.resize(triangulation.n_active_cells(), 0U);
1884 
1885  // In a first step, obtain the weights of the locally owned
1886  // cells. For all others, the weight remains at the zero the
1887  // vector was initialized with above.
1888  for (const auto &cell : triangulation.active_cell_iterators() |
1890  cell_weights[cell->active_cell_index()] =
1891  triangulation.signals.weight(cell, CellStatus::cell_will_persist);
1892 
1893  // If this is a parallel triangulation, we then need to also
1894  // get the weights for all other cells. We have asserted above
1895  // that this function can't be used for
1896  // parallel::distribute::Triangulation objects, so the only
1897  // ones we have to worry about here are
1898  // parallel::shared::Triangulation
1899  if (const auto shared_tria =
1901  &triangulation))
1902  Utilities::MPI::sum(cell_weights,
1903  shared_tria->get_communicator(),
1904  cell_weights);
1905 
1906  // verify that the global sum of weights is larger than 0
1907  Assert(std::accumulate(cell_weights.begin(),
1908  cell_weights.end(),
1909  std::uint64_t(0)) > 0,
1910  ExcMessage("The global sum of weights over all active cells "
1911  "is zero. Please verify how you generate weights."));
1912  }
1913 
1914  // Call the other more general function
1915  partition_triangulation(n_partitions,
1916  cell_weights,
1917  cell_connection_graph,
1918  triangulation,
1919  partitioner);
1920  }
1921 
1922 
1923 
1924  template <int dim, int spacedim>
1925  void
1926  partition_triangulation(const unsigned int n_partitions,
1927  const std::vector<unsigned int> &cell_weights,
1928  const SparsityPattern &cell_connection_graph,
1930  const SparsityTools::Partitioner partitioner)
1931  {
1933  &triangulation) == nullptr),
1934  ExcMessage("Objects of type parallel::distributed::Triangulation "
1935  "are already partitioned implicitly and can not be "
1936  "partitioned again explicitly."));
1937  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
1938  Assert(cell_connection_graph.n_rows() == triangulation.n_active_cells(),
1939  ExcMessage("Connectivity graph has wrong size"));
1940  Assert(cell_connection_graph.n_cols() == triangulation.n_active_cells(),
1941  ExcMessage("Connectivity graph has wrong size"));
1942 
1943  // signal that partitioning is going to happen
1944  triangulation.signals.pre_partition();
1945 
1946  // check for an easy return
1947  if (n_partitions == 1)
1948  {
1949  for (const auto &cell : triangulation.active_cell_iterators())
1950  cell->set_subdomain_id(0);
1951  return;
1952  }
1953 
1954  // partition this connection graph and get
1955  // back a vector of indices, one per degree
1956  // of freedom (which is associated with a
1957  // cell)
1958  std::vector<unsigned int> partition_indices(triangulation.n_active_cells());
1959  SparsityTools::partition(cell_connection_graph,
1960  cell_weights,
1961  n_partitions,
1962  partition_indices,
1963  partitioner);
1964 
1965  // finally loop over all cells and set the subdomain ids
1966  for (const auto &cell : triangulation.active_cell_iterators())
1967  cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
1968  }
1969 
1970 
1971  namespace internal
1972  {
1976  template <class IT>
1977  void
1979  unsigned int &current_proc_idx,
1980  unsigned int &current_cell_idx,
1981  const unsigned int n_active_cells,
1982  const unsigned int n_partitions)
1983  {
1984  if (cell->is_active())
1985  {
1986  while (current_cell_idx >=
1987  std::floor(static_cast<uint_least64_t>(n_active_cells) *
1988  (current_proc_idx + 1) / n_partitions))
1989  ++current_proc_idx;
1990  cell->set_subdomain_id(current_proc_idx);
1991  ++current_cell_idx;
1992  }
1993  else
1994  {
1995  for (unsigned int n = 0; n < cell->n_children(); ++n)
1997  current_proc_idx,
1998  current_cell_idx,
2000  n_partitions);
2001  }
2002  }
2003  } // namespace internal
2004 
2005  template <int dim, int spacedim>
2006  void
2007  partition_triangulation_zorder(const unsigned int n_partitions,
2009  const bool group_siblings)
2010  {
2012  &triangulation) == nullptr),
2013  ExcMessage("Objects of type parallel::distributed::Triangulation "
2014  "are already partitioned implicitly and can not be "
2015  "partitioned again explicitly."));
2016  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
2017  Assert(triangulation.signals.weight.empty(), ExcNotImplemented());
2018 
2019  // signal that partitioning is going to happen
2020  triangulation.signals.pre_partition();
2021 
2022  // check for an easy return
2023  if (n_partitions == 1)
2024  {
2025  for (const auto &cell : triangulation.active_cell_iterators())
2026  cell->set_subdomain_id(0);
2027  return;
2028  }
2029 
2030  // Duplicate the coarse cell reordoring
2031  // as done in p4est
2032  std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
2033  std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
2034 
2035  DynamicSparsityPattern cell_connectivity;
2037  0,
2038  cell_connectivity);
2039  coarse_cell_to_p4est_tree_permutation.resize(triangulation.n_cells(0));
2040  SparsityTools::reorder_hierarchical(cell_connectivity,
2041  coarse_cell_to_p4est_tree_permutation);
2042 
2043  p4est_tree_to_coarse_cell_permutation =
2044  Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
2045 
2046  unsigned int current_proc_idx = 0;
2047  unsigned int current_cell_idx = 0;
2048  const unsigned int n_active_cells = triangulation.n_active_cells();
2049 
2050  // set subdomain id for active cell descendants
2051  // of each coarse cell in permuted order
2052  for (unsigned int idx = 0; idx < triangulation.n_cells(0); ++idx)
2053  {
2054  const unsigned int coarse_cell_idx =
2055  p4est_tree_to_coarse_cell_permutation[idx];
2056  typename Triangulation<dim, spacedim>::cell_iterator coarse_cell(
2057  &triangulation, 0, coarse_cell_idx);
2058 
2060  current_proc_idx,
2061  current_cell_idx,
2063  n_partitions);
2064  }
2065 
2066  // if all children of a cell are active (e.g. we
2067  // have a cell that is refined once and no part
2068  // is refined further), p4est places all of them
2069  // on the same processor. The new owner will be
2070  // the processor with the largest number of children
2071  // (ties are broken by picking the lower rank).
2072  // Duplicate this logic here.
2073  if (group_siblings)
2074  {
2076  cell = triangulation.begin(),
2077  endc = triangulation.end();
2078  for (; cell != endc; ++cell)
2079  {
2080  if (cell->is_active())
2081  continue;
2082  bool all_children_active = true;
2083  std::map<unsigned int, unsigned int> map_cpu_n_cells;
2084  for (unsigned int n = 0; n < cell->n_children(); ++n)
2085  if (!cell->child(n)->is_active())
2086  {
2087  all_children_active = false;
2088  break;
2089  }
2090  else
2091  ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
2092 
2093  if (!all_children_active)
2094  continue;
2095 
2096  unsigned int new_owner = cell->child(0)->subdomain_id();
2097  for (std::map<unsigned int, unsigned int>::iterator it =
2098  map_cpu_n_cells.begin();
2099  it != map_cpu_n_cells.end();
2100  ++it)
2101  if (it->second > map_cpu_n_cells[new_owner])
2102  new_owner = it->first;
2103 
2104  for (unsigned int n = 0; n < cell->n_children(); ++n)
2105  cell->child(n)->set_subdomain_id(new_owner);
2106  }
2107  }
2108  }
2109 
2110 
2111  template <int dim, int spacedim>
2112  void
2114  {
2115  unsigned int n_levels = triangulation.n_levels();
2116  for (int lvl = n_levels - 1; lvl >= 0; --lvl)
2117  {
2118  for (const auto &cell : triangulation.cell_iterators_on_level(lvl))
2119  {
2120  if (cell->is_active())
2121  cell->set_level_subdomain_id(cell->subdomain_id());
2122  else
2123  {
2124  Assert(cell->child(0)->level_subdomain_id() !=
2126  ExcInternalError());
2127  cell->set_level_subdomain_id(
2128  cell->child(0)->level_subdomain_id());
2129  }
2130  }
2131  }
2132  }
2133 
2134  namespace internal
2135  {
2136  namespace
2137  {
2138  // Split get_subdomain_association() for p::d::T since we want to compile
2139  // it in 1d but none of the p4est stuff is available in 1d.
2140  template <int dim, int spacedim>
2141  void
2144  &triangulation,
2145  const std::vector<CellId> &cell_ids,
2146  std::vector<types::subdomain_id> &subdomain_ids)
2147  {
2148 #ifndef DEAL_II_WITH_P4EST
2149  (void)triangulation;
2150  (void)cell_ids;
2151  (void)subdomain_ids;
2152  Assert(
2153  false,
2154  ExcMessage(
2155  "You are attempting to use a functionality that is only available "
2156  "if deal.II was configured to use p4est, but cmake did not find a "
2157  "valid p4est library."));
2158 #else
2159  // for parallel distributed triangulations, we will ask the p4est oracle
2160  // about the global partitioning of active cells since this information
2161  // is stored on every process
2162  for (const auto &cell_id : cell_ids)
2163  {
2164  // find descendent from coarse quadrant
2165  typename ::internal::p4est::types<dim>::quadrant p4est_cell,
2167 
2168  ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
2169  for (const auto &child_index : cell_id.get_child_indices())
2170  {
2171  ::internal::p4est::init_quadrant_children<dim>(
2172  p4est_cell, p4est_children);
2173  p4est_cell =
2174  p4est_children[static_cast<unsigned int>(child_index)];
2175  }
2176 
2177  // find owning process, i.e., the subdomain id
2178  const int owner =
2180  const_cast<typename ::internal::p4est::types<dim>::forest
2181  *>(triangulation.get_p4est()),
2182  cell_id.get_coarse_cell_id(),
2183  &p4est_cell,
2185  triangulation.get_communicator()));
2186 
2187  Assert(owner >= 0, ExcMessage("p4est should know the owner."));
2188 
2189  subdomain_ids.push_back(owner);
2190  }
2191 #endif
2192  }
2193 
2194 
2195 
2196  template <int spacedim>
2197  void
2200  const std::vector<CellId> &,
2201  std::vector<types::subdomain_id> &)
2202  {
2203  Assert(false, ExcNotImplemented());
2204  }
2205  } // anonymous namespace
2206  } // namespace internal
2207 
2208 
2209 
2210  template <int dim, int spacedim>
2211  std::vector<types::subdomain_id>
2213  const std::vector<CellId> &cell_ids)
2214  {
2215  std::vector<types::subdomain_id> subdomain_ids;
2216  subdomain_ids.reserve(cell_ids.size());
2217 
2218  if (dynamic_cast<
2220  &triangulation) != nullptr)
2221  {
2222  Assert(false, ExcNotImplemented());
2223  }
2225  *parallel_tria = dynamic_cast<
2227  &triangulation))
2228  {
2229  internal::get_subdomain_association(*parallel_tria,
2230  cell_ids,
2231  subdomain_ids);
2232  }
2233  else if (const parallel::shared::Triangulation<dim, spacedim> *shared_tria =
2235  *>(&triangulation))
2236  {
2237  // for parallel shared triangulations, we need to access true subdomain
2238  // ids which are also valid for artificial cells
2239  const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
2240  shared_tria->get_true_subdomain_ids_of_cells();
2241 
2242  for (const auto &cell_id : cell_ids)
2243  {
2244  const unsigned int active_cell_index =
2245  shared_tria->create_cell_iterator(cell_id)->active_cell_index();
2246  subdomain_ids.push_back(
2247  true_subdomain_ids_of_cells[active_cell_index]);
2248  }
2249  }
2250  else
2251  {
2252  // the most general type of triangulation is the serial one. here, all
2253  // subdomain information is directly available
2254  for (const auto &cell_id : cell_ids)
2255  {
2256  subdomain_ids.push_back(
2257  triangulation.create_cell_iterator(cell_id)->subdomain_id());
2258  }
2259  }
2260 
2261  return subdomain_ids;
2262  }
2263 
2264 
2265 
2266  template <int dim, int spacedim>
2267  void
2269  std::vector<types::subdomain_id> &subdomain)
2270  {
2271  Assert(subdomain.size() == triangulation.n_active_cells(),
2272  ExcDimensionMismatch(subdomain.size(),
2273  triangulation.n_active_cells()));
2274  for (const auto &cell : triangulation.active_cell_iterators())
2275  subdomain[cell->active_cell_index()] = cell->subdomain_id();
2276  }
2277 
2278 
2279 
2280  template <int dim, int spacedim>
2281  unsigned int
2284  const types::subdomain_id subdomain)
2285  {
2286  unsigned int count = 0;
2287  for (const auto &cell : triangulation.active_cell_iterators())
2288  if (cell->subdomain_id() == subdomain)
2289  ++count;
2290 
2291  return count;
2292  }
2293 
2294 
2295 
2296  template <int dim, int spacedim>
2297  std::vector<bool>
2299  {
2300  // start with all vertices
2301  std::vector<bool> locally_owned_vertices =
2302  triangulation.get_used_vertices();
2303 
2304  // if the triangulation is distributed, eliminate those that
2305  // are owned by other processors -- either because the vertex is
2306  // on an artificial cell, or because it is on a ghost cell with
2307  // a smaller subdomain
2308  if (const auto *tr = dynamic_cast<
2310  &triangulation))
2311  for (const auto &cell : triangulation.active_cell_iterators())
2312  if (cell->is_artificial() ||
2313  (cell->is_ghost() &&
2314  (cell->subdomain_id() < tr->locally_owned_subdomain())))
2315  for (const unsigned int v : cell->vertex_indices())
2316  locally_owned_vertices[cell->vertex_index(v)] = false;
2317 
2318  return locally_owned_vertices;
2319  }
2320 
2321 
2322 
2323  namespace internal
2324  {
2325  namespace FixUpDistortedChildCells
2326  {
2327  // compute the mean square
2328  // deviation of the alternating
2329  // forms of the children of the
2330  // given object from that of
2331  // the object itself. for
2332  // objects with
2333  // structdim==spacedim, the
2334  // alternating form is the
2335  // determinant of the jacobian,
2336  // whereas for faces with
2337  // structdim==spacedim-1, the
2338  // alternating form is the
2339  // (signed and scaled) normal
2340  // vector
2341  //
2342  // this average square
2343  // deviation is computed for an
2344  // object where the center node
2345  // has been replaced by the
2346  // second argument to this
2347  // function
2348  template <typename Iterator, int spacedim>
2349  double
2350  objective_function(const Iterator &object,
2351  const Point<spacedim> &object_mid_point)
2352  {
2353  const unsigned int structdim =
2354  Iterator::AccessorType::structure_dimension;
2355  Assert(spacedim == Iterator::AccessorType::dimension,
2356  ExcInternalError());
2357 
2358  // everything below is wrong
2359  // if not for the following
2360  // condition
2361  Assert(object->refinement_case() ==
2363  ExcNotImplemented());
2364  // first calculate the
2365  // average alternating form
2366  // for the parent cell/face
2369  Tensor<spacedim - structdim, spacedim>
2370  parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
2371 
2372  for (const unsigned int i : object->vertex_indices())
2373  parent_vertices[i] = object->vertex(i);
2374 
2376  parent_vertices, parent_alternating_forms);
2377 
2378  const Tensor<spacedim - structdim, spacedim>
2379  average_parent_alternating_form =
2380  std::accumulate(parent_alternating_forms,
2381  parent_alternating_forms +
2384 
2385  // now do the same
2386  // computation for the
2387  // children where we use the
2388  // given location for the
2389  // object mid point instead of
2390  // the one the triangulation
2391  // currently reports
2395  Tensor<spacedim - structdim, spacedim> child_alternating_forms
2398 
2399  for (unsigned int c = 0; c < object->n_children(); ++c)
2400  for (const unsigned int i : object->child(c)->vertex_indices())
2401  child_vertices[c][i] = object->child(c)->vertex(i);
2402 
2403  // replace mid-object
2404  // vertex. note that for
2405  // child i, the mid-object
2406  // vertex happens to have the
2407  // number
2408  // max_children_per_cell-i
2409  for (unsigned int c = 0; c < object->n_children(); ++c)
2410  child_vertices[c][GeometryInfo<structdim>::max_children_per_cell - c -
2411  1] = object_mid_point;
2412 
2413  for (unsigned int c = 0; c < object->n_children(); ++c)
2415  child_vertices[c], child_alternating_forms[c]);
2416 
2417  // on a uniformly refined
2418  // hypercube object, the child
2419  // alternating forms should
2420  // all be smaller by a factor
2421  // of 2^structdim than the
2422  // ones of the parent. as a
2423  // consequence, we'll use the
2424  // squared deviation from
2425  // this ideal value as an
2426  // objective function
2427  double objective = 0;
2428  for (unsigned int c = 0; c < object->n_children(); ++c)
2429  for (const unsigned int i : object->child(c)->vertex_indices())
2430  objective +=
2431  (child_alternating_forms[c][i] -
2432  average_parent_alternating_form / std::pow(2., 1. * structdim))
2433  .norm_square();
2434 
2435  return objective;
2436  }
2437 
2438 
2444  template <typename Iterator>
2446  get_face_midpoint(const Iterator &object,
2447  const unsigned int f,
2448  std::integral_constant<int, 1>)
2449  {
2450  return object->vertex(f);
2451  }
2452 
2453 
2454 
2460  template <typename Iterator>
2462  get_face_midpoint(const Iterator &object,
2463  const unsigned int f,
2464  std::integral_constant<int, 2>)
2465  {
2466  return object->line(f)->center();
2467  }
2468 
2469 
2470 
2476  template <typename Iterator>
2478  get_face_midpoint(const Iterator &object,
2479  const unsigned int f,
2480  std::integral_constant<int, 3>)
2481  {
2482  return object->face(f)->center();
2483  }
2484 
2485 
2486 
2509  template <typename Iterator>
2510  double
2511  minimal_diameter(const Iterator &object)
2512  {
2513  const unsigned int structdim =
2514  Iterator::AccessorType::structure_dimension;
2515 
2516  double diameter = object->diameter();
2517  for (const unsigned int f : object->face_indices())
2518  for (unsigned int e = f + 1; e < object->n_faces(); ++e)
2519  diameter = std::min(
2520  diameter,
2521  get_face_midpoint(object,
2522  f,
2523  std::integral_constant<int, structdim>())
2524  .distance(get_face_midpoint(
2525  object, e, std::integral_constant<int, structdim>())));
2526 
2527  return diameter;
2528  }
2529 
2530 
2531 
2536  template <typename Iterator>
2537  bool
2538  fix_up_object(const Iterator &object)
2539  {
2540  const unsigned int structdim =
2541  Iterator::AccessorType::structure_dimension;
2542  const unsigned int spacedim = Iterator::AccessorType::space_dimension;
2543 
2544  // right now we can only deal with cells that have been refined
2545  // isotropically because that is the only case where we have a cell
2546  // mid-point that can be moved around without having to consider
2547  // boundary information
2548  Assert(object->has_children(), ExcInternalError());
2549  Assert(object->refinement_case() ==
2551  ExcNotImplemented());
2552 
2553  // get the current location of the object mid-vertex:
2554  Point<spacedim> object_mid_point = object->child(0)->vertex(
2556 
2557  // now do a few steepest descent steps to reduce the objective
2558  // function. compute the diameter in the helper function above
2559  unsigned int iteration = 0;
2560  const double diameter = minimal_diameter(object);
2561 
2562  // current value of objective function and initial delta
2563  double current_value = objective_function(object, object_mid_point);
2564  double initial_delta = 0;
2565 
2566  do
2567  {
2568  // choose a step length that is initially 1/4 of the child
2569  // objects' diameter, and a sequence whose sum does not converge
2570  // (to avoid premature termination of the iteration)
2571  const double step_length = diameter / 4 / (iteration + 1);
2572 
2573  // compute the objective function's derivative using a two-sided
2574  // difference formula with eps=step_length/10
2575  Tensor<1, spacedim> gradient;
2576  for (unsigned int d = 0; d < spacedim; ++d)
2577  {
2578  const double eps = step_length / 10;
2579 
2581  h[d] = eps / 2;
2582 
2583  gradient[d] =
2585  object, project_to_object(object, object_mid_point + h)) -
2587  object, project_to_object(object, object_mid_point - h))) /
2588  eps;
2589  }
2590 
2591  // there is nowhere to go
2592  if (gradient.norm() == 0)
2593  break;
2594 
2595  // We need to go in direction -gradient. the optimal value of the
2596  // objective function is zero, so assuming that the model is
2597  // quadratic we would have to go -2*val/||gradient|| in this
2598  // direction, make sure we go at most step_length into this
2599  // direction
2600  object_mid_point -=
2601  std::min(2 * current_value / (gradient * gradient),
2602  step_length / gradient.norm()) *
2603  gradient;
2604  object_mid_point = project_to_object(object, object_mid_point);
2605 
2606  // compute current value of the objective function
2607  const double previous_value = current_value;
2608  current_value = objective_function(object, object_mid_point);
2609 
2610  if (iteration == 0)
2611  initial_delta = (previous_value - current_value);
2612 
2613  // stop if we aren't moving much any more
2614  if ((iteration >= 1) &&
2615  ((previous_value - current_value < 0) ||
2616  (std::fabs(previous_value - current_value) <
2617  0.001 * initial_delta)))
2618  break;
2619 
2620  ++iteration;
2621  }
2622  while (iteration < 20);
2623 
2624  // verify that the new
2625  // location is indeed better
2626  // than the one before. check
2627  // this by comparing whether
2628  // the minimum value of the
2629  // products of parent and
2630  // child alternating forms is
2631  // positive. for cells this
2632  // means that the
2633  // determinants have the same
2634  // sign, for faces that the
2635  // face normals of parent and
2636  // children point in the same
2637  // general direction
2638  double old_min_product, new_min_product;
2639 
2642  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
2643  parent_vertices[i] = object->vertex(i);
2644 
2645  Tensor<spacedim - structdim, spacedim>
2646  parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
2648  parent_vertices, parent_alternating_forms);
2649 
2653 
2654  for (unsigned int c = 0; c < object->n_children(); ++c)
2655  for (const unsigned int i : object->child(c)->vertex_indices())
2656  child_vertices[c][i] = object->child(c)->vertex(i);
2657 
2658  Tensor<spacedim - structdim, spacedim> child_alternating_forms
2661 
2662  for (unsigned int c = 0; c < object->n_children(); ++c)
2664  child_vertices[c], child_alternating_forms[c]);
2665 
2666  old_min_product =
2667  child_alternating_forms[0][0] * parent_alternating_forms[0];
2668  for (unsigned int c = 0; c < object->n_children(); ++c)
2669  for (const unsigned int i : object->child(c)->vertex_indices())
2670  for (const unsigned int j : object->vertex_indices())
2671  old_min_product = std::min<double>(old_min_product,
2672  child_alternating_forms[c][i] *
2673  parent_alternating_forms[j]);
2674 
2675  // for the new minimum value,
2676  // replace mid-object
2677  // vertex. note that for child
2678  // i, the mid-object vertex
2679  // happens to have the number
2680  // max_children_per_cell-i
2681  for (unsigned int c = 0; c < object->n_children(); ++c)
2682  child_vertices[c][GeometryInfo<structdim>::max_children_per_cell - c -
2683  1] = object_mid_point;
2684 
2685  for (unsigned int c = 0; c < object->n_children(); ++c)
2687  child_vertices[c], child_alternating_forms[c]);
2688 
2689  new_min_product =
2690  child_alternating_forms[0][0] * parent_alternating_forms[0];
2691  for (unsigned int c = 0; c < object->n_children(); ++c)
2692  for (const unsigned int i : object->child(c)->vertex_indices())
2693  for (const unsigned int j : object->vertex_indices())
2694  new_min_product = std::min<double>(new_min_product,
2695  child_alternating_forms[c][i] *
2696  parent_alternating_forms[j]);
2697 
2698  // if new minimum value is
2699  // better than before, then set the
2700  // new mid point. otherwise
2701  // return this object as one of
2702  // those that can't apparently
2703  // be fixed
2704  if (new_min_product >= old_min_product)
2705  object->child(0)->vertex(
2707  object_mid_point;
2708 
2709  // return whether after this
2710  // operation we have an object that
2711  // is well oriented
2712  return (std::max(new_min_product, old_min_product) > 0);
2713  }
2714 
2715 
2716 
2717  // possibly fix up the faces of a cell by moving around its mid-points
2718  template <int dim, int spacedim>
2719  void
2721  const typename ::Triangulation<dim, spacedim>::cell_iterator
2722  &cell,
2723  std::integral_constant<int, dim>,
2724  std::integral_constant<int, spacedim>)
2725  {
2726  // see if we first can fix up some of the faces of this object. We can
2727  // mess with faces if and only if the neighboring cell is not even
2728  // more refined than we are (since in that case the sub-faces have
2729  // themselves children that we can't move around any more). however,
2730  // the latter case shouldn't happen anyway: if the current face is
2731  // distorted but the neighbor is even more refined, then the face had
2732  // been deformed before already, and had been ignored at the time; we
2733  // should then also be able to ignore it this time as well
2734  for (auto f : cell->face_indices())
2735  {
2736  Assert(cell->face(f)->has_children(), ExcInternalError());
2737  Assert(cell->face(f)->refinement_case() ==
2739  ExcInternalError());
2740 
2741  bool subface_is_more_refined = false;
2742  for (unsigned int g = 0;
2743  g < GeometryInfo<dim>::max_children_per_face;
2744  ++g)
2745  if (cell->face(f)->child(g)->has_children())
2746  {
2747  subface_is_more_refined = true;
2748  break;
2749  }
2750 
2751  if (subface_is_more_refined == true)
2752  continue;
2753 
2754  // we finally know that we can do something about this face
2755  fix_up_object(cell->face(f));
2756  }
2757  }
2758  } /* namespace FixUpDistortedChildCells */
2759  } /* namespace internal */
2760 
2761 
2762  template <int dim, int spacedim>
2766  &distorted_cells,
2767  Triangulation<dim, spacedim> & /*triangulation*/)
2768  {
2769  static_assert(
2770  dim != 1 && spacedim != 1,
2771  "This function is only valid when dim != 1 or spacedim != 1.");
2772  typename Triangulation<dim, spacedim>::DistortedCellList unfixable_subset;
2773 
2774  // loop over all cells that we have to fix up
2775  for (typename std::list<
2776  typename Triangulation<dim, spacedim>::cell_iterator>::const_iterator
2777  cell_ptr = distorted_cells.distorted_cells.begin();
2778  cell_ptr != distorted_cells.distorted_cells.end();
2779  ++cell_ptr)
2780  {
2781  const typename Triangulation<dim, spacedim>::cell_iterator &cell =
2782  *cell_ptr;
2783 
2784  Assert(!cell->is_active(),
2785  ExcMessage(
2786  "This function is only valid for a list of cells that "
2787  "have children (i.e., no cell in the list may be active)."));
2788 
2790  cell,
2791  std::integral_constant<int, dim>(),
2792  std::integral_constant<int, spacedim>());
2793 
2794  // If possible, fix up the object.
2796  unfixable_subset.distorted_cells.push_back(cell);
2797  }
2798 
2799  return unfixable_subset;
2800  }
2801 
2802 
2803 
2804  template <int dim, int spacedim>
2805  void
2807  const bool reset_boundary_ids)
2808  {
2809  const auto src_boundary_ids = tria.get_boundary_ids();
2810  std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
2811  auto m_it = dst_manifold_ids.begin();
2812  for (const auto b : src_boundary_ids)
2813  {
2814  *m_it = static_cast<types::manifold_id>(b);
2815  ++m_it;
2816  }
2817  const std::vector<types::boundary_id> reset_boundary_id =
2818  reset_boundary_ids ?
2819  std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
2820  src_boundary_ids;
2821  map_boundary_to_manifold_ids(src_boundary_ids,
2822  dst_manifold_ids,
2823  tria,
2824  reset_boundary_id);
2825  }
2826 
2827 
2828 
2829  template <int dim, int spacedim>
2830  void
2832  const std::vector<types::boundary_id> &src_boundary_ids,
2833  const std::vector<types::manifold_id> &dst_manifold_ids,
2835  const std::vector<types::boundary_id> &reset_boundary_ids_)
2836  {
2837  AssertDimension(src_boundary_ids.size(), dst_manifold_ids.size());
2838  const auto reset_boundary_ids =
2839  reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
2840  AssertDimension(reset_boundary_ids.size(), src_boundary_ids.size());
2841 
2842  // in 3d, we not only have to copy boundary ids of faces, but also of edges
2843  // because we see them twice (once from each adjacent boundary face),
2844  // we cannot immediately reset their boundary ids. thus, copy first
2845  // and reset later
2846  if (dim >= 3)
2847  for (const auto &cell : tria.active_cell_iterators())
2848  for (auto f : cell->face_indices())
2849  if (cell->face(f)->at_boundary())
2850  for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
2851  {
2852  const auto bid = cell->face(f)->line(e)->boundary_id();
2853  const unsigned int ind = std::find(src_boundary_ids.begin(),
2854  src_boundary_ids.end(),
2855  bid) -
2856  src_boundary_ids.begin();
2857  if (ind < src_boundary_ids.size())
2858  cell->face(f)->line(e)->set_manifold_id(
2859  dst_manifold_ids[ind]);
2860  }
2861 
2862  // now do cells
2863  for (const auto &cell : tria.active_cell_iterators())
2864  for (auto f : cell->face_indices())
2865  if (cell->face(f)->at_boundary())
2866  {
2867  const auto bid = cell->face(f)->boundary_id();
2868  const unsigned int ind =
2869  std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
2870  src_boundary_ids.begin();
2871 
2872  if (ind < src_boundary_ids.size())
2873  {
2874  // assign the manifold id
2875  cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
2876  // then reset boundary id
2877  cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
2878  }
2879 
2880  if (dim >= 3)
2881  for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
2882  {
2883  const auto bid = cell->face(f)->line(e)->boundary_id();
2884  const unsigned int ind = std::find(src_boundary_ids.begin(),
2885  src_boundary_ids.end(),
2886  bid) -
2887  src_boundary_ids.begin();
2888  if (ind < src_boundary_ids.size())
2889  cell->face(f)->line(e)->set_boundary_id(
2890  reset_boundary_ids[ind]);
2891  }
2892  }
2893  }
2894 
2895 
2896  template <int dim, int spacedim>
2897  void
2899  const bool compute_face_ids)
2900  {
2902  cell = tria.begin_active(),
2903  endc = tria.end();
2904 
2905  for (; cell != endc; ++cell)
2906  {
2907  cell->set_manifold_id(cell->material_id());
2908  if (compute_face_ids == true)
2909  {
2910  for (auto f : cell->face_indices())
2911  {
2912  if (cell->at_boundary(f) == false)
2913  cell->face(f)->set_manifold_id(
2914  std::min(cell->material_id(),
2915  cell->neighbor(f)->material_id()));
2916  else
2917  cell->face(f)->set_manifold_id(cell->material_id());
2918  }
2919  }
2920  }
2921  }
2922 
2923 
2924  template <int dim, int spacedim>
2925  void
2928  const std::function<types::manifold_id(
2929  const std::set<types::manifold_id> &)> &disambiguation_function,
2930  bool overwrite_only_flat_manifold_ids)
2931  {
2932  // Easy case first:
2933  if (dim == 1)
2934  return;
2935  const unsigned int n_subobjects =
2936  dim == 2 ? tria.n_lines() : tria.n_lines() + tria.n_quads();
2937 
2938  // If user index is zero, then it has not been set.
2939  std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
2940  std::vector<unsigned int> backup;
2941  tria.save_user_indices(backup);
2943 
2944  unsigned next_index = 1;
2945  for (auto &cell : tria.active_cell_iterators())
2946  {
2947  if (dim > 1)
2948  for (unsigned int l = 0; l < cell->n_lines(); ++l)
2949  {
2950  if (cell->line(l)->user_index() == 0)
2951  {
2952  AssertIndexRange(next_index, n_subobjects + 1);
2953  manifold_ids[next_index].insert(cell->manifold_id());
2954  cell->line(l)->set_user_index(next_index++);
2955  }
2956  else
2957  manifold_ids[cell->line(l)->user_index()].insert(
2958  cell->manifold_id());
2959  }
2960  if (dim > 2)
2961  for (unsigned int l = 0; l < cell->n_faces(); ++l)
2962  {
2963  if (cell->quad(l)->user_index() == 0)
2964  {
2965  AssertIndexRange(next_index, n_subobjects + 1);
2966  manifold_ids[next_index].insert(cell->manifold_id());
2967  cell->quad(l)->set_user_index(next_index++);
2968  }
2969  else
2970  manifold_ids[cell->quad(l)->user_index()].insert(
2971  cell->manifold_id());
2972  }
2973  }
2974  for (auto &cell : tria.active_cell_iterators())
2975  {
2976  if (dim > 1)
2977  for (unsigned int l = 0; l < cell->n_lines(); ++l)
2978  {
2979  const auto id = cell->line(l)->user_index();
2980  // Make sure we change the manifold indicator only once
2981  if (id != 0)
2982  {
2983  if (cell->line(l)->manifold_id() ==
2985  overwrite_only_flat_manifold_ids == false)
2986  cell->line(l)->set_manifold_id(
2987  disambiguation_function(manifold_ids[id]));
2988  cell->line(l)->set_user_index(0);
2989  }
2990  }
2991  if (dim > 2)
2992  for (unsigned int l = 0; l < cell->n_faces(); ++l)
2993  {
2994  const auto id = cell->quad(l)->user_index();
2995  // Make sure we change the manifold indicator only once
2996  if (id != 0)
2997  {
2998  if (cell->quad(l)->manifold_id() ==
3000  overwrite_only_flat_manifold_ids == false)
3001  cell->quad(l)->set_manifold_id(
3002  disambiguation_function(manifold_ids[id]));
3003  cell->quad(l)->set_user_index(0);
3004  }
3005  }
3006  }
3007  tria.load_user_indices(backup);
3008  }
3009 
3010 
3011 
3012  template <int dim, int spacedim>
3013  void
3015  const double limit_angle_fraction)
3016  {
3017  if (dim == 1)
3018  return; // Nothing to do
3019 
3020  // Check that we don't have hanging nodes
3022  ExcMessage("The input Triangulation cannot "
3023  "have hanging nodes."));
3024 
3026 
3027  bool has_cells_with_more_than_dim_faces_on_boundary = true;
3028  bool has_cells_with_dim_faces_on_boundary = false;
3029 
3030  unsigned int refinement_cycles = 0;
3031 
3032  while (has_cells_with_more_than_dim_faces_on_boundary)
3033  {
3034  has_cells_with_more_than_dim_faces_on_boundary = false;
3035 
3036  for (const auto &cell : tria.active_cell_iterators())
3037  {
3038  unsigned int boundary_face_counter = 0;
3039  for (auto f : cell->face_indices())
3040  if (cell->face(f)->at_boundary())
3041  boundary_face_counter++;
3042  if (boundary_face_counter > dim)
3043  {
3044  has_cells_with_more_than_dim_faces_on_boundary = true;
3045  break;
3046  }
3047  else if (boundary_face_counter == dim)
3048  has_cells_with_dim_faces_on_boundary = true;
3049  }
3050  if (has_cells_with_more_than_dim_faces_on_boundary)
3051  {
3052  tria.refine_global(1);
3053  refinement_cycles++;
3054  }
3055  }
3056 
3057  if (has_cells_with_dim_faces_on_boundary)
3058  {
3059  tria.refine_global(1);
3060  refinement_cycles++;
3061  }
3062  else
3063  {
3064  while (refinement_cycles > 0)
3065  {
3066  for (const auto &cell : tria.active_cell_iterators())
3067  cell->set_coarsen_flag();
3069  refinement_cycles--;
3070  }
3071  return;
3072  }
3073 
3074  std::vector<bool> cells_to_remove(tria.n_active_cells(), false);
3075  std::vector<Point<spacedim>> vertices = tria.get_vertices();
3076 
3077  std::vector<bool> faces_to_remove(tria.n_raw_faces(), false);
3078 
3079  std::vector<CellData<dim>> cells_to_add;
3080  SubCellData subcelldata_to_add;
3081 
3082  // Trick compiler for dimension independent things
3083  const unsigned int v0 = 0, v1 = 1, v2 = (dim > 1 ? 2 : 0),
3084  v3 = (dim > 1 ? 3 : 0);
3085 
3086  for (const auto &cell : tria.active_cell_iterators())
3087  {
3088  double angle_fraction = 0;
3089  unsigned int vertex_at_corner = numbers::invalid_unsigned_int;
3090 
3091  if (dim == 2)
3092  {
3094  p0[spacedim > 1 ? 1 : 0] = 1;
3096  p1[0] = 1;
3097 
3098  if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
3099  {
3100  p0 = cell->vertex(v0) - cell->vertex(v2);
3101  p1 = cell->vertex(v3) - cell->vertex(v2);
3102  vertex_at_corner = v2;
3103  }
3104  else if (cell->face(v3)->at_boundary() &&
3105  cell->face(v1)->at_boundary())
3106  {
3107  p0 = cell->vertex(v2) - cell->vertex(v3);
3108  p1 = cell->vertex(v1) - cell->vertex(v3);
3109  vertex_at_corner = v3;
3110  }
3111  else if (cell->face(1)->at_boundary() &&
3112  cell->face(2)->at_boundary())
3113  {
3114  p0 = cell->vertex(v0) - cell->vertex(v1);
3115  p1 = cell->vertex(v3) - cell->vertex(v1);
3116  vertex_at_corner = v1;
3117  }
3118  else if (cell->face(2)->at_boundary() &&
3119  cell->face(0)->at_boundary())
3120  {
3121  p0 = cell->vertex(v2) - cell->vertex(v0);
3122  p1 = cell->vertex(v1) - cell->vertex(v0);
3123  vertex_at_corner = v0;
3124  }
3125  p0 /= p0.norm();
3126  p1 /= p1.norm();
3127  angle_fraction = std::acos(p0 * p1) / numbers::PI;
3128  }
3129  else
3130  {
3131  Assert(false, ExcNotImplemented());
3132  }
3133 
3134  if (angle_fraction > limit_angle_fraction)
3135  {
3136  auto flags_removal = [&](unsigned int f1,
3137  unsigned int f2,
3138  unsigned int n1,
3139  unsigned int n2) -> void {
3140  cells_to_remove[cell->active_cell_index()] = true;
3141  cells_to_remove[cell->neighbor(n1)->active_cell_index()] = true;
3142  cells_to_remove[cell->neighbor(n2)->active_cell_index()] = true;
3143 
3144  faces_to_remove[cell->face(f1)->index()] = true;
3145  faces_to_remove[cell->face(f2)->index()] = true;
3146 
3147  faces_to_remove[cell->neighbor(n1)->face(f1)->index()] = true;
3148  faces_to_remove[cell->neighbor(n2)->face(f2)->index()] = true;
3149  };
3150 
3151  auto cell_creation = [&](const unsigned int vv0,
3152  const unsigned int vv1,
3153  const unsigned int f0,
3154  const unsigned int f1,
3155 
3156  const unsigned int n0,
3157  const unsigned int v0n0,
3158  const unsigned int v1n0,
3159 
3160  const unsigned int n1,
3161  const unsigned int v0n1,
3162  const unsigned int v1n1) {
3163  CellData<dim> c1, c2;
3164  CellData<1> l1, l2;
3165 
3166  c1.vertices[v0] = cell->vertex_index(vv0);
3167  c1.vertices[v1] = cell->vertex_index(vv1);
3168  c1.vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
3169  c1.vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
3170 
3171  c1.manifold_id = cell->manifold_id();
3172  c1.material_id = cell->material_id();
3173 
3174  c2.vertices[v0] = cell->vertex_index(vv0);
3175  c2.vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
3176  c2.vertices[v2] = cell->vertex_index(vv1);
3177  c2.vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
3178 
3179  c2.manifold_id = cell->manifold_id();
3180  c2.material_id = cell->material_id();
3181 
3182  l1.vertices[0] = cell->vertex_index(vv0);
3183  l1.vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
3184 
3185  l1.boundary_id = cell->line(f0)->boundary_id();
3186  l1.manifold_id = cell->line(f0)->manifold_id();
3187  subcelldata_to_add.boundary_lines.push_back(l1);
3188 
3189  l2.vertices[0] = cell->vertex_index(vv0);
3190  l2.vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
3191 
3192  l2.boundary_id = cell->line(f1)->boundary_id();
3193  l2.manifold_id = cell->line(f1)->manifold_id();
3194  subcelldata_to_add.boundary_lines.push_back(l2);
3195 
3196  cells_to_add.push_back(c1);
3197  cells_to_add.push_back(c2);
3198  };
3199 
3200  if (dim == 2)
3201  {
3202  switch (vertex_at_corner)
3203  {
3204  case 0:
3205  flags_removal(0, 2, 3, 1);
3206  cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
3207  break;
3208  case 1:
3209  flags_removal(1, 2, 3, 0);
3210  cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
3211  break;
3212  case 2:
3213  flags_removal(3, 0, 1, 2);
3214  cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
3215  break;
3216  case 3:
3217  flags_removal(3, 1, 0, 2);
3218  cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
3219  break;
3220  }
3221  }
3222  else
3223  {
3224  Assert(false, ExcNotImplemented());
3225  }
3226  }
3227  }
3228 
3229  // if no cells need to be added, then no regularization is necessary.
3230  // Restore things as they were before this function was called.
3231  if (cells_to_add.empty())
3232  {
3233  while (refinement_cycles > 0)
3234  {
3235  for (const auto &cell : tria.active_cell_iterators())
3236  cell->set_coarsen_flag();
3238  refinement_cycles--;
3239  }
3240  return;
3241  }
3242 
3243  // add the cells that were not marked as skipped
3244  for (const auto &cell : tria.active_cell_iterators())
3245  {
3246  if (cells_to_remove[cell->active_cell_index()] == false)
3247  {
3248  CellData<dim> c(cell->n_vertices());
3249  for (const unsigned int v : cell->vertex_indices())
3250  c.vertices[v] = cell->vertex_index(v);
3251  c.manifold_id = cell->manifold_id();
3252  c.material_id = cell->material_id();
3253  cells_to_add.push_back(c);
3254  }
3255  }
3256 
3257  // Face counter for both dim == 2 and dim == 3
3259  face = tria.begin_active_face(),
3260  endf = tria.end_face();
3261  for (; face != endf; ++face)
3262  if ((face->at_boundary() ||
3263  face->manifold_id() != numbers::flat_manifold_id) &&
3264  faces_to_remove[face->index()] == false)
3265  {
3266  for (unsigned int l = 0; l < face->n_lines(); ++l)
3267  {
3268  CellData<1> line;
3269  if (dim == 2)
3270  {
3271  for (const unsigned int v : face->vertex_indices())
3272  line.vertices[v] = face->vertex_index(v);
3273  line.boundary_id = face->boundary_id();
3274  line.manifold_id = face->manifold_id();
3275  }
3276  else
3277  {
3278  for (const unsigned int v : face->line(l)->vertex_indices())
3279  line.vertices[v] = face->line(l)->vertex_index(v);
3280  line.boundary_id = face->line(l)->boundary_id();
3281  line.manifold_id = face->line(l)->manifold_id();
3282  }
3283  subcelldata_to_add.boundary_lines.push_back(line);
3284  }
3285  if (dim == 3)
3286  {
3287  CellData<2> quad(face->n_vertices());
3288  for (const unsigned int v : face->vertex_indices())
3289  quad.vertices[v] = face->vertex_index(v);
3290  quad.boundary_id = face->boundary_id();
3291  quad.manifold_id = face->manifold_id();
3292  subcelldata_to_add.boundary_quads.push_back(quad);
3293  }
3294  }
3296  cells_to_add,
3297  subcelldata_to_add);
3299 
3300  // Save manifolds
3301  auto manifold_ids = tria.get_manifold_ids();
3302  std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
3303  manifolds;
3304  // Set manifolds in new Triangulation
3305  for (const auto manifold_id : manifold_ids)
3307  manifolds[manifold_id] = tria.get_manifold(manifold_id).clone();
3308 
3309  tria.clear();
3310 
3311  tria.create_triangulation(vertices, cells_to_add, subcelldata_to_add);
3312 
3313  // Restore manifolds
3314  for (const auto manifold_id : manifold_ids)
3316  tria.set_manifold(manifold_id, *manifolds[manifold_id]);
3317  }
3318 
3319 
3320 
3321  template <int dim, int spacedim>
3322 #ifndef DOXYGEN
3323  std::tuple<
3324  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3325  std::vector<std::vector<Point<dim>>>,
3326  std::vector<std::vector<unsigned int>>>
3327 #else
3328  return_type
3329 #endif
3331  const Cache<dim, spacedim> &cache,
3332  const std::vector<Point<spacedim>> &points,
3334  &cell_hint)
3335  {
3336  const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
3337  // Splitting the tuple's components
3338  auto &cells = std::get<0>(cqmp);
3339  auto &qpoints = std::get<1>(cqmp);
3340  auto &maps = std::get<2>(cqmp);
3341 
3342  return std::make_tuple(std::move(cells),
3343  std::move(qpoints),
3344  std::move(maps));
3345  }
3346 
3347 
3348 
3349  template <int dim, int spacedim>
3350 #ifndef DOXYGEN
3351  std::tuple<
3352  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3353  std::vector<std::vector<Point<dim>>>,
3354  std::vector<std::vector<unsigned int>>,
3355  std::vector<unsigned int>>
3356 #else
3357  return_type
3358 #endif
3360  const Cache<dim, spacedim> &cache,
3361  const std::vector<Point<spacedim>> &points,
3363  &cell_hint)
3364  {
3365  Assert((dim == spacedim),
3366  ExcMessage("Only implemented for dim==spacedim."));
3367 
3368  // Alias
3369  namespace bgi = boost::geometry::index;
3370 
3371  // Get the mapping
3372  const auto &mapping = cache.get_mapping();
3373 
3374  // How many points are here?
3375  const unsigned int np = points.size();
3376 
3377  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
3378  cells_out;
3379  std::vector<std::vector<Point<dim>>> qpoints_out;
3380  std::vector<std::vector<unsigned int>> maps_out;
3381  std::vector<unsigned int> missing_points_out;
3382 
3383  // Now the easy case.
3384  if (np == 0)
3385  return std::make_tuple(std::move(cells_out),
3386  std::move(qpoints_out),
3387  std::move(maps_out),
3388  std::move(missing_points_out));
3389 
3390  // For the search we shall use the following tree
3391  const auto &b_tree = cache.get_cell_bounding_boxes_rtree();
3392 
3393  // Now make a tree of indices for the points
3394  // [TODO] This would work better with pack_rtree_of_indices, but
3395  // windows does not like it. Build a tree with pairs of point and id
3396  std::vector<std::pair<Point<spacedim>, unsigned int>> points_and_ids(np);
3397  for (unsigned int i = 0; i < np; ++i)
3398  points_and_ids[i] = std::make_pair(points[i], i);
3399  const auto p_tree = pack_rtree(points_and_ids);
3400 
3401  // Keep track of all found points
3402  std::vector<bool> found_points(points.size(), false);
3403 
3404  // Check if a point was found
3405  const auto already_found = [&found_points](const auto &id) {
3406  AssertIndexRange(id.second, found_points.size());
3407  return found_points[id.second];
3408  };
3409 
3410  // check if the given cell was already in the vector of cells before. If so,
3411  // insert in the corresponding vectors the reference point and the id.
3412  // Otherwise append a new entry to all vectors.
3413  const auto store_cell_point_and_id =
3414  [&](
3416  const Point<dim> &ref_point,
3417  const unsigned int &id) {
3418  const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
3419  if (it != cells_out.rend())
3420  {
3421  const auto cell_id =
3422  (cells_out.size() - 1 - (it - cells_out.rbegin()));
3423  qpoints_out[cell_id].emplace_back(ref_point);
3424  maps_out[cell_id].emplace_back(id);
3425  }
3426  else
3427  {
3428  cells_out.emplace_back(cell);
3429  qpoints_out.emplace_back(std::vector<Point<dim>>({ref_point}));
3430  maps_out.emplace_back(std::vector<unsigned int>({id}));
3431  }
3432  };
3433 
3434  // Check all points within a given pair of box and cell
3435  const auto check_all_points_within_box = [&](const auto &leaf) {
3436  const double relative_tolerance = 1e-12;
3437  const BoundingBox<spacedim> box =
3438  leaf.first.create_extended_relative(relative_tolerance);
3439  const auto &cell_hint = leaf.second;
3440 
3441  for (const auto &point_and_id :
3442  p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
3443  bgi::intersects(box)))
3444  {
3445  const auto id = point_and_id.second;
3446  const auto cell_and_ref =
3448  points[id],
3449  cell_hint);
3450  const auto &cell = cell_and_ref.first;
3451  const auto &ref_point = cell_and_ref.second;
3452 
3453  if (cell.state() == IteratorState::valid)
3454  store_cell_point_and_id(cell, ref_point, id);
3455  else
3456  missing_points_out.emplace_back(id);
3457 
3458  // Don't look anymore for this point
3459  found_points[id] = true;
3460  }
3461  };
3462 
3463  // If a hint cell was given, use it
3464  if (cell_hint.state() == IteratorState::valid)
3465  check_all_points_within_box(
3466  std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
3467 
3468  // Now loop over all points that have not been found yet
3469  for (unsigned int i = 0; i < np; ++i)
3470  if (found_points[i] == false)
3471  {
3472  // Get the closest cell to this point
3473  const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
3474  // Now checks all points that fall within this box
3475  if (leaf != b_tree.qend())
3476  check_all_points_within_box(*leaf);
3477  else
3478  {
3479  // We should not get here. Throw an error.
3480  Assert(false, ExcInternalError());
3481  }
3482  }
3483  // Now make sure we send out the rest of the points that we did not find.
3484  for (unsigned int i = 0; i < np; ++i)
3485  if (found_points[i] == false)
3486  missing_points_out.emplace_back(i);
3487 
3488  // Debug Checking
3489  AssertDimension(cells_out.size(), maps_out.size());
3490  AssertDimension(cells_out.size(), qpoints_out.size());
3491 
3492 #ifdef DEBUG
3493  unsigned int c = cells_out.size();
3494  unsigned int qps = 0;
3495  // The number of points in all
3496  // the cells must be the same as
3497  // the number of points we
3498  // started off from,
3499  // plus the points which were ignored
3500  for (unsigned int n = 0; n < c; ++n)
3501  {
3502  AssertDimension(qpoints_out[n].size(), maps_out[n].size());
3503  qps += qpoints_out[n].size();
3504  }
3505 
3506  Assert(qps + missing_points_out.size() == np,
3507  ExcDimensionMismatch(qps + missing_points_out.size(), np));
3508 #endif
3509 
3510  return std::make_tuple(std::move(cells_out),
3511  std::move(qpoints_out),
3512  std::move(maps_out),
3513  std::move(missing_points_out));
3514  }
3515 
3516 
3517 
3518  template <int dim, int spacedim>
3519 #ifndef DOXYGEN
3520  std::tuple<
3521  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3522  std::vector<std::vector<Point<dim>>>,
3523  std::vector<std::vector<unsigned int>>,
3524  std::vector<std::vector<Point<spacedim>>>,
3525  std::vector<std::vector<unsigned int>>>
3526 #else
3527  return_type
3528 #endif
3530  const GridTools::Cache<dim, spacedim> &cache,
3531  const std::vector<Point<spacedim>> &points,
3532  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3533  const double tolerance,
3534  const std::vector<bool> &marked_vertices,
3535  const bool enforce_unique_mapping)
3536  {
3537  // run internal function ...
3538  const auto all =
3540  points,
3541  global_bboxes,
3542  marked_vertices,
3543  tolerance,
3544  false,
3545  enforce_unique_mapping)
3546  .send_components;
3547 
3548  // ... and reshuffle the data
3549  std::tuple<
3550  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3551  std::vector<std::vector<Point<dim>>>,
3552  std::vector<std::vector<unsigned int>>,
3553  std::vector<std::vector<Point<spacedim>>>,
3554  std::vector<std::vector<unsigned int>>>
3555  result;
3556 
3557  std::pair<int, int> dummy{-1, -1};
3558 
3559  for (unsigned int i = 0; i < all.size(); ++i)
3560  {
3561  if (dummy != std::get<0>(all[i]))
3562  {
3563  std::get<0>(result).push_back(
3565  &cache.get_triangulation(),
3566  std::get<0>(all[i]).first,
3567  std::get<0>(all[i]).second});
3568 
3569  const unsigned int new_size = std::get<0>(result).size();
3570 
3571  std::get<1>(result).resize(new_size);
3572  std::get<2>(result).resize(new_size);
3573  std::get<3>(result).resize(new_size);
3574  std::get<4>(result).resize(new_size);
3575 
3576  dummy = std::get<0>(all[i]);
3577  }
3578 
3579  std::get<1>(result).back().push_back(
3580  std::get<3>(all[i])); // reference point
3581  std::get<2>(result).back().push_back(std::get<2>(all[i])); // index
3582  std::get<3>(result).back().push_back(std::get<4>(all[i])); // real point
3583  std::get<4>(result).back().push_back(std::get<1>(all[i])); // rank
3584  }
3585 
3586  return result;
3587  }
3588 
3589 
3590 
3591  namespace internal
3592  {
3599  template <int spacedim, typename T>
3600  std::tuple<std::vector<unsigned int>,
3601  std::vector<unsigned int>,
3602  std::vector<unsigned int>>
3604  const MPI_Comm comm,
3605  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3606  const std::vector<T> &entities,
3607  const double tolerance)
3608  {
3609  std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
3610  ranks_and_indices.reserve(entities.size());
3611 
3612 #if defined(DEAL_II_WITH_ARBORX)
3613  static constexpr bool use_arborx = true;
3614 #else
3615  static constexpr bool use_arborx = false;
3616 #endif
3617  // Lambda to process bboxes if global_bboxes.size()>1 or ArborX not
3618  // available
3619  const auto process_bboxes = [&]() -> void {
3620  std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes_temp;
3621  auto *global_bboxes_to_be_used = &global_bboxes;
3622 
3623  if (global_bboxes.size() == 1 && use_arborx == false)
3624  {
3625  global_bboxes_temp =
3626  Utilities::MPI::all_gather(comm, global_bboxes[0]);
3627  global_bboxes_to_be_used = &global_bboxes_temp;
3628  }
3629 
3630  // helper function to determine if a bounding box is valid
3631  const auto is_valid = [](const auto &bb) {
3632  for (unsigned int i = 0; i < spacedim; ++i)
3633  if (bb.get_boundary_points().first[i] >
3634  bb.get_boundary_points().second[i])
3635  return false;
3636 
3637  return true;
3638  };
3639 
3640  // linearize vector of vectors
3641  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
3642  boxes_and_ranks;
3643 
3644  for (unsigned rank = 0; rank < global_bboxes_to_be_used->size(); ++rank)
3645  for (const auto &box : (*global_bboxes_to_be_used)[rank])
3646  if (is_valid(box))
3647  boxes_and_ranks.emplace_back(box, rank);
3648 
3649  // pack boxes into r-tree
3650  const auto tree = pack_rtree(boxes_and_ranks);
3651 
3652  // loop over all entities
3653  for (unsigned int i = 0; i < entities.size(); ++i)
3654  {
3655  // create a bounding box with tolerance
3656  const auto bb =
3657  BoundingBox<spacedim>(entities[i]).create_extended(tolerance);
3658 
3659  // determine ranks potentially owning point/bounding box
3660  std::set<unsigned int> my_ranks;
3661 
3662  for (const auto &box_and_rank :
3663  tree | boost::geometry::index::adaptors::queried(
3664  boost::geometry::index::intersects(bb)))
3665  my_ranks.insert(box_and_rank.second);
3666 
3667  for (const auto rank : my_ranks)
3668  ranks_and_indices.emplace_back(rank, i);
3669  }
3670  };
3671 
3672  if constexpr (use_arborx)
3673  {
3674  if (global_bboxes.size() == 1)
3675  {
3676  ArborXWrappers::DistributedTree distributed_tree(
3677  comm, global_bboxes[0]);
3678  std::vector<BoundingBox<spacedim>> query_bounding_boxes;
3679  for (const auto &entity : entities)
3680  query_bounding_boxes.emplace_back(
3681  BoundingBox<spacedim>(entity).create_extended(tolerance));
3682 
3684  query_bounding_boxes);
3685  const auto &[indices_ranks, offsets] =
3686  distributed_tree.query(bb_intersect);
3687 
3688  for (unsigned long int i = 0; i < offsets.size() - 1; ++i)
3689  {
3690  std::set<unsigned int> my_ranks;
3691  for (int j = offsets[i]; j < offsets[i + 1]; ++j)
3692  my_ranks.insert(indices_ranks[j].second);
3693 
3694  for (const auto rank : my_ranks)
3695  ranks_and_indices.emplace_back(rank, i);
3696  }
3697  }
3698  else
3699  {
3700  // global_bboxes.size()>1
3701  process_bboxes();
3702  }
3703  }
3704  else
3705  {
3706  // No ArborX
3707  process_bboxes();
3708  }
3709 
3710 
3711  // convert to CRS
3712  std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
3713 
3714  std::vector<unsigned int> ranks;
3715  std::vector<unsigned int> ptr;
3716  std::vector<unsigned int> indices;
3717 
3718  unsigned int current_rank = numbers::invalid_unsigned_int;
3719 
3720  for (const std::pair<unsigned int, unsigned int> &i : ranks_and_indices)
3721  {
3722  if (current_rank != i.first)
3723  {
3724  current_rank = i.first;
3725  ranks.push_back(current_rank);
3726  ptr.push_back(indices.size());
3727  }
3728 
3729  indices.push_back(i.second);
3730  }
3731  ptr.push_back(indices.size());
3732 
3733  return {std::move(ranks), std::move(ptr), std::move(indices)};
3734  }
3735 
3736 
3737 
3738  template <int dim, int spacedim>
3739  std::vector<
3740  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
3741  Point<dim>>>
3743  const Cache<dim, spacedim> &cache,
3744  const Point<spacedim> &point,
3746  const std::vector<bool> &marked_vertices,
3747  const double tolerance,
3748  const bool enforce_unique_mapping)
3749  {
3750  std::vector<
3751  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
3752  Point<dim>>>
3753  locally_owned_active_cells_around_point;
3754 
3755  const auto first_cell = GridTools::find_active_cell_around_point(
3756  cache.get_mapping(),
3757  cache.get_triangulation(),
3758  point,
3759  cache.get_vertex_to_cell_map(),
3761  cell_hint,
3762  marked_vertices,
3763  cache.get_used_vertices_rtree(),
3764  tolerance,
3766 
3767  const unsigned int my_rank = Utilities::MPI::this_mpi_process(
3769 
3770  cell_hint = first_cell.first;
3771  if (cell_hint.state() == IteratorState::valid)
3772  {
3773  const auto active_cells_around_point =
3775  cache.get_mapping(),
3776  cache.get_triangulation(),
3777  point,
3778  tolerance,
3779  first_cell,
3780  &cache.get_vertex_to_cell_map());
3781 
3782  if (enforce_unique_mapping)
3783  {
3784  // check if the rank of this process is the lowest of all cells
3785  // if not, the other process will handle this cell and we don't
3786  // have to do here anything in the case of unique mapping
3787  unsigned int lowes_rank = numbers::invalid_unsigned_int;
3788 
3789  for (const auto &cell : active_cells_around_point)
3790  lowes_rank = std::min(lowes_rank, cell.first->subdomain_id());
3791 
3792  if (lowes_rank != my_rank)
3793  return {};
3794  }
3795 
3796  locally_owned_active_cells_around_point.reserve(
3797  active_cells_around_point.size());
3798 
3799  for (const auto &cell : active_cells_around_point)
3800  if (cell.first->is_locally_owned())
3801  locally_owned_active_cells_around_point.push_back(cell);
3802  }
3803 
3804  std::sort(locally_owned_active_cells_around_point.begin(),
3805  locally_owned_active_cells_around_point.end(),
3806  [](const auto &a, const auto &b) { return a.first < b.first; });
3807 
3808  if (enforce_unique_mapping &&
3809  locally_owned_active_cells_around_point.size() > 1)
3810  // in the case of unique mapping, we only need a single cell
3811  return {locally_owned_active_cells_around_point.front()};
3812  else
3813  return locally_owned_active_cells_around_point;
3814  }
3815 
3816  template <int dim, int spacedim>
3819  : n_searched_points(numbers::invalid_unsigned_int)
3820  {}
3821 
3822  template <int dim, int spacedim>
3823  void
3825  {
3826  // before reshuffeling the data check if data.recv_components and
3827  // n_searched_points are in a valid state.
3828  Assert(n_searched_points != numbers::invalid_unsigned_int,
3829  ExcInternalError());
3830  Assert(recv_components.empty() ||
3831  std::get<1>(*std::max_element(recv_components.begin(),
3832  recv_components.end(),
3833  [](const auto &a, const auto &b) {
3834  return std::get<1>(a) <
3835  std::get<1>(b);
3836  })) < n_searched_points,
3837  ExcInternalError());
3838 
3839  send_ranks.clear();
3840  recv_ranks.clear();
3841  send_ptrs.clear();
3842  recv_ptrs.clear();
3843 
3844  if (true)
3845  {
3846  // sort according to rank (and point index and cell) -> make
3847  // deterministic
3848  std::sort(send_components.begin(),
3849  send_components.end(),
3850  [&](const auto &a, const auto &b) {
3851  if (std::get<1>(a) != std::get<1>(b)) // rank
3852  return std::get<1>(a) < std::get<1>(b);
3853 
3854  if (std::get<2>(a) != std::get<2>(b)) // point index
3855  return std::get<2>(a) < std::get<2>(b);
3856 
3857  return std::get<0>(a) < std::get<0>(b); // cell
3858  });
3859 
3860  // perform enumeration and extract rank information
3861  for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
3862  i < send_components.size();
3863  ++i)
3864  {
3865  std::get<5>(send_components[i]) = i;
3866 
3867  if (dummy != std::get<1>(send_components[i]))
3868  {
3869  dummy = std::get<1>(send_components[i]);
3870  send_ranks.push_back(dummy);
3871  send_ptrs.push_back(i);
3872  }
3873  }
3874  send_ptrs.push_back(send_components.size());
3875 
3876  // sort according to cell, rank, point index (while keeping
3877  // partial ordering)
3878  std::sort(send_components.begin(),
3879  send_components.end(),
3880  [&](const auto &a, const auto &b) {
3881  if (std::get<0>(a) != std::get<0>(b))
3882  return std::get<0>(a) < std::get<0>(b); // cell
3883 
3884  if (std::get<1>(a) != std::get<1>(b))
3885  return std::get<1>(a) < std::get<1>(b); // rank
3886 
3887  if (std::get<2>(a) != std::get<2>(b))
3888  return std::get<2>(a) < std::get<2>(b); // point index
3889 
3890  return std::get<5>(a) < std::get<5>(b); // enumeration
3891  });
3892  }
3893 
3894  if (recv_components.size() > 0)
3895  {
3896  // sort according to rank (and point index) -> make deterministic
3897  std::sort(recv_components.begin(),
3898  recv_components.end(),
3899  [&](const auto &a, const auto &b) {
3900  if (std::get<0>(a) != std::get<0>(b))
3901  return std::get<0>(a) < std::get<0>(b); // rank
3902 
3903  return std::get<1>(a) < std::get<1>(b); // point index
3904  });
3905 
3906  // perform enumeration and extract rank information
3907  for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
3908  i < recv_components.size();
3909  ++i)
3910  {
3911  std::get<2>(recv_components[i]) = i;
3912 
3913  if (dummy != std::get<0>(recv_components[i]))
3914  {
3915  dummy = std::get<0>(recv_components[i]);
3916  recv_ranks.push_back(dummy);
3917  recv_ptrs.push_back(i);
3918  }
3919  }
3920  recv_ptrs.push_back(recv_components.size());
3921 
3922  // sort according to point index and rank (while keeping partial
3923  // ordering)
3924  std::sort(recv_components.begin(),
3925  recv_components.end(),
3926  [&](const auto &a, const auto &b) {
3927  if (std::get<1>(a) != std::get<1>(b))
3928  return std::get<1>(a) < std::get<1>(b); // point index
3929 
3930  if (std::get<0>(a) != std::get<0>(b))
3931  return std::get<0>(a) < std::get<0>(b); // rank
3932 
3933  return std::get<2>(a) < std::get<2>(b); // enumeration
3934  });
3935  }
3936  }
3937 
3938 
3939 
3940  template <int dim, int spacedim>
3943  const GridTools::Cache<dim, spacedim> &cache,
3944  const std::vector<Point<spacedim>> &points,
3945  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3946  const std::vector<bool> &marked_vertices,
3947  const double tolerance,
3948  const bool perform_handshake,
3949  const bool enforce_unique_mapping)
3950  {
3952  result.n_searched_points = points.size();
3953 
3954  auto &send_components = result.send_components;
3955  auto &recv_components = result.recv_components;
3956 
3957  const auto comm = cache.get_triangulation().get_communicator();
3958 
3959  const auto potential_owners = internal::guess_owners_of_entities(
3960  comm, global_bboxes, points, tolerance);
3961 
3962  const auto &potential_owners_ranks = std::get<0>(potential_owners);
3963  const auto &potential_owners_ptrs = std::get<1>(potential_owners);
3964  const auto &potential_owners_indices = std::get<2>(potential_owners);
3965 
3966  auto cell_hint = cache.get_triangulation().begin_active();
3967 
3968  const auto translate = [&](const unsigned int other_rank) {
3969  const auto ptr = std::find(potential_owners_ranks.begin(),
3970  potential_owners_ranks.end(),
3971  other_rank);
3972 
3973  Assert(ptr != potential_owners_ranks.end(), ExcInternalError());
3974 
3975  const auto other_rank_index =
3976  std::distance(potential_owners_ranks.begin(), ptr);
3977 
3978  return other_rank_index;
3979  };
3980 
3981  Assert(
3982  (marked_vertices.empty()) ||
3983  (marked_vertices.size() == cache.get_triangulation().n_vertices()),
3984  ExcMessage(
3985  "The marked_vertices vector has to be either empty or its size has "
3986  "to equal the number of vertices of the triangulation."));
3987 
3988  using RequestType = std::vector<std::pair<unsigned int, Point<spacedim>>>;
3989  using AnswerType = std::vector<unsigned int>;
3990 
3991  // In the case that a marked_vertices vector has been given and none
3992  // of its entries is true, we know that this process does not own
3993  // any of the incoming points (and it will not send any data) so
3994  // that we can take a short cut.
3995  const bool has_relevant_vertices =
3996  (marked_vertices.empty()) ||
3997  (std::find(marked_vertices.begin(), marked_vertices.end(), true) !=
3998  marked_vertices.end());
3999 
4000  const auto create_request = [&](const unsigned int other_rank) {
4001  const auto other_rank_index = translate(other_rank);
4002 
4003  RequestType request;
4004  request.reserve(potential_owners_ptrs[other_rank_index + 1] -
4005  potential_owners_ptrs[other_rank_index]);
4006 
4007  for (unsigned int i = potential_owners_ptrs[other_rank_index];
4008  i < potential_owners_ptrs[other_rank_index + 1];
4009  ++i)
4010  request.emplace_back(potential_owners_indices[i],
4011  points[potential_owners_indices[i]]);
4012 
4013  return request;
4014  };
4015 
4016  const auto answer_request =
4017  [&](const unsigned int &other_rank,
4018  const RequestType &request) -> AnswerType {
4019  AnswerType answer(request.size(), 0);
4020 
4021  if (has_relevant_vertices)
4022  {
4023  cell_hint = cache.get_triangulation().begin_active();
4024 
4025  for (unsigned int i = 0; i < request.size(); ++i)
4026  {
4027  const auto &index_and_point = request[i];
4028 
4029  const auto cells_and_reference_positions =
4031  cache,
4032  index_and_point.second,
4033  cell_hint,
4034  marked_vertices,
4035  tolerance,
4036  enforce_unique_mapping);
4037 
4038  if (cell_hint.state() != IteratorState::valid)
4039  cell_hint = cache.get_triangulation().begin_active();
4040 
4041  for (const auto &cell_and_reference_position :
4042  cells_and_reference_positions)
4043  {
4044  const auto cell = cell_and_reference_position.first;
4045  auto reference_position =
4046  cell_and_reference_position.second;
4047 
4048  reference_position =
4049  cell->reference_cell().closest_point(reference_position);
4050 
4051  send_components.emplace_back(
4052  std::pair<int, int>(cell->level(), cell->index()),
4053  other_rank,
4054  index_and_point.first,
4055  reference_position,
4056  index_and_point.second,
4058  }
4059 
4060  answer[i] = cells_and_reference_positions.size();
4061  }
4062  }
4063 
4064  if (perform_handshake)
4065  return answer;
4066  else
4067  return {};
4068  };
4069 
4070  const auto process_answer = [&](const unsigned int other_rank,
4071  const AnswerType &answer) {
4072  if (perform_handshake)
4073  {
4074  const auto other_rank_index = translate(other_rank);
4075 
4076  for (unsigned int i = 0; i < answer.size(); ++i)
4077  for (unsigned int j = 0; j < answer[i]; ++j)
4078  recv_components.emplace_back(
4079  other_rank,
4080  potential_owners_indices
4081  [i + potential_owners_ptrs[other_rank_index]],
4083  }
4084  };
4085 
4086  Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
4087  potential_owners_ranks,
4088  create_request,
4089  answer_request,
4090  process_answer,
4091  comm);
4092 
4093  result.finalize_setup();
4094 
4095  return result;
4096  }
4097 
4098 
4099 
4100  template <int structdim, int spacedim>
4101  template <int dim>
4102  DistributedComputePointLocationsInternal<dim, spacedim>
4105  const unsigned int n_points_1D,
4107  const Mapping<dim, spacedim> &mapping,
4108  const bool consistent_numbering_of_sender_and_receiver) const
4109  {
4110  using CellIterator =
4112 
4114  spacedim>
4115  result;
4116 
4117  // We need quadrature rules for the intersections. We are using a
4118  // QGaussSimplex quadrature rule since CGAL always returns simplices
4119  // as intersections.
4120  const QGaussSimplex<structdim> quadrature(n_points_1D);
4121 
4122  // Resulting quadrature points get different indices. In the case the
4123  // requested intersections are unique also the resulting quadrature
4124  // points are unique and we can simply number the points in an
4125  // ascending way.
4126  for (const auto &recv_component : recv_components)
4127  {
4128  // dependent on the size of the intersection an empty quadrature
4129  // is returned. Therefore, we have to compute the quadrature also
4130  // here.
4131  const Quadrature<spacedim> &quad =
4132  quadrature.compute_affine_transformation(
4133  std::get<2>(recv_component));
4134 
4135  for (unsigned int i = 0; i < quad.size(); ++i)
4136  {
4137  // the third component of result.recv_components is not needed
4138  // before finalize_setup.
4139  result.recv_components.emplace_back(
4140  std::get<0>(recv_component),
4141  result.recv_components.size(), // number of point
4143  }
4144  }
4145 
4146  // since empty quadratures might be present we have to set the number
4147  // of searched points after inserting the point indices into
4148  // recv_components
4149  result.n_searched_points = result.recv_components.size();
4150 
4151  // send_ranks, counter, and indices_of_rank is only needed if
4152  // consistent_numbering_of_sender_and_receiver==true
4153  // indices_of_rank is always empty if deal.II is compiled without MPI
4154  std::map<unsigned int, std::vector<unsigned int>> indices_of_rank;
4155  std::map<unsigned int, unsigned int> counter;
4156  std::set<unsigned int> send_ranks;
4157  if (consistent_numbering_of_sender_and_receiver)
4158  {
4159  for (const auto &sc : send_components)
4160  send_ranks.insert(std::get<1>(sc));
4161 
4162  for (const auto rank : send_ranks)
4163  counter[rank] = 0;
4164 
4165  // indices assigned at recv side needed to fill send_components
4166  indices_of_rank = communicate_indices(result.recv_components,
4168  }
4169 
4170  for (const auto &send_component : send_components)
4171  {
4172  const CellIterator cell(&tria,
4173  std::get<0>(send_component).first,
4174  std::get<0>(send_component).second);
4175 
4176  const Quadrature<spacedim> &quad =
4177  quadrature.compute_affine_transformation(
4178  std::get<3>(send_component));
4179 
4180  const auto rank = std::get<1>(send_component);
4181 
4182  for (unsigned int q = 0; q < quad.size(); ++q)
4183  {
4184  // the fifth component of result.send_components is filled
4185  // during sorting the data and initializing the CRS structures
4186  result.send_components.emplace_back(std::make_tuple(
4187  std::get<0>(send_component),
4188  rank,
4189  indices_of_rank.empty() ?
4190  result.send_components.size() :
4191  indices_of_rank.at(rank)[counter.at(rank)],
4192  mapping.transform_real_to_unit_cell(cell, quad.point(q)),
4193  quad.point(q),
4195 
4196  if (!indices_of_rank.empty())
4197  ++counter[rank];
4198  }
4199  }
4200 
4201  result.finalize_setup();
4202 
4203  return result;
4204  }
4205 
4206 
4207 
4208  template <int structdim, int spacedim>
4209  std::map<unsigned int, std::vector<unsigned int>>
4212  const std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
4213  &point_recv_components,
4214  const MPI_Comm comm) const
4215  {
4216 #ifndef DEAL_II_WITH_MPI
4217  Assert(false, ExcNeedsMPI());
4218  (void)point_recv_components;
4219  (void)comm;
4220  return {};
4221 #else
4222  // since we are converting to DistributedComputePointLocationsInternal
4223  // we use the RPE tag
4224  const auto mpi_tag =
4226 
4227  const unsigned int my_rank = Utilities::MPI::this_mpi_process(comm);
4228 
4229  std::set<unsigned int> send_ranks;
4230  for (const auto &sc : send_components)
4231  send_ranks.insert(std::get<1>(sc));
4232  std::set<unsigned int> recv_ranks;
4233  for (const auto &rc : recv_components)
4234  recv_ranks.insert(std::get<0>(rc));
4235 
4236  std::vector<MPI_Request> requests;
4237  requests.reserve(send_ranks.size());
4238 
4239  // rank to used indices on the rank needed on sending side
4240  std::map<unsigned int, std::vector<unsigned int>> indices_of_rank;
4241  indices_of_rank[my_rank] = std::vector<unsigned int>();
4242 
4243  // rank to used indices on the rank known on recv side
4244  std::map<unsigned int, std::vector<unsigned int>> send_indices_of_rank;
4245  for (const auto rank : recv_ranks)
4246  if (rank != my_rank)
4247  send_indices_of_rank[rank] = std::vector<unsigned int>();
4248 
4249  // fill the maps
4250  for (const auto &point_recv_component : point_recv_components)
4251  {
4252  const auto rank = std::get<0>(point_recv_component);
4253  const auto idx = std::get<1>(point_recv_component);
4254 
4255  if (rank == my_rank)
4256  indices_of_rank[rank].emplace_back(idx);
4257  else
4258  send_indices_of_rank[rank].emplace_back(idx);
4259  }
4260 
4261  // send indices to the ranks we normally receive from
4262  for (const auto rank : recv_ranks)
4263  {
4264  if (rank == my_rank)
4265  continue;
4266 
4267  auto buffer = Utilities::pack(send_indices_of_rank[rank], false);
4268 
4269  requests.push_back(MPI_Request());
4270 
4271  const int ierr = MPI_Isend(buffer.data(),
4272  buffer.size(),
4273  MPI_CHAR,
4274  rank,
4275  mpi_tag,
4276  comm,
4277  &requests.back());
4278  AssertThrowMPI(ierr);
4279  }
4280 
4281  // receive indices at the ranks we normally send from
4282  for (const auto rank : send_ranks)
4283  {
4284  if (rank == my_rank)
4285  continue;
4286 
4287  MPI_Status status;
4288  int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, comm, &status);
4289  AssertThrowMPI(ierr);
4290 
4291  int message_length;
4292  ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
4293  AssertThrowMPI(ierr);
4294 
4295  std::vector<char> buffer(message_length);
4296 
4297  ierr = MPI_Recv(buffer.data(),
4298  buffer.size(),
4299  MPI_CHAR,
4300  status.MPI_SOURCE,
4301  mpi_tag,
4302  comm,
4303  MPI_STATUS_IGNORE);
4304  AssertThrowMPI(ierr);
4305 
4306  indices_of_rank[status.MPI_SOURCE] =
4307  Utilities::unpack<std::vector<unsigned int>>(buffer, false);
4308  }
4309 
4310  // make sure all messages have been sent
4311  const int ierr =
4312  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4313  AssertThrowMPI(ierr);
4314 
4315  return indices_of_rank;
4316 #endif
4317  }
4318 
4319 
4320 
4321  template <int structdim, int dim, int spacedim>
4324  const Cache<dim, spacedim> &cache,
4325  const std::vector<std::vector<Point<spacedim>>> &intersection_requests,
4326  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
4327  const std::vector<bool> &marked_vertices,
4328  const double tolerance)
4329  {
4330  using IntersectionRequest = std::vector<Point<spacedim>>;
4331 
4332  using IntersectionAnswer =
4334  structdim,
4335  spacedim>::IntersectionType;
4336 
4337  const auto comm = cache.get_triangulation().get_communicator();
4338 
4340  result;
4341 
4342  auto &send_components = result.send_components;
4343  auto &recv_components = result.recv_components;
4344  auto &recv_ptrs = result.recv_ptrs;
4345 
4346  // search for potential owners
4347  const auto potential_owners = internal::guess_owners_of_entities(
4348  comm, global_bboxes, intersection_requests, tolerance);
4349 
4350  const auto &potential_owners_ranks = std::get<0>(potential_owners);
4351  const auto &potential_owners_ptrs = std::get<1>(potential_owners);
4352  const auto &potential_owners_indices = std::get<2>(potential_owners);
4353 
4354  const auto translate = [&](const unsigned int other_rank) {
4355  const auto ptr = std::find(potential_owners_ranks.begin(),
4356  potential_owners_ranks.end(),
4357  other_rank);
4358 
4359  Assert(ptr != potential_owners_ranks.end(), ExcInternalError());
4360 
4361  const auto other_rank_index =
4362  std::distance(potential_owners_ranks.begin(), ptr);
4363 
4364  return other_rank_index;
4365  };
4366 
4367  Assert(
4368  (marked_vertices.empty()) ||
4369  (marked_vertices.size() == cache.get_triangulation().n_vertices()),
4370  ExcMessage(
4371  "The marked_vertices vector has to be either empty or its size has "
4372  "to equal the number of vertices of the triangulation."));
4373 
4374  // In the case that a marked_vertices vector has been given and none
4375  // of its entries is true, we know that this process does not own
4376  // any of the incoming points (and it will not send any data) so
4377  // that we can take a short cut.
4378  const bool has_relevant_vertices =
4379  (marked_vertices.empty()) ||
4380  (std::find(marked_vertices.begin(), marked_vertices.end(), true) !=
4381  marked_vertices.end());
4382 
4383  // intersection between two cells:
4384  // One rank requests all intersections of owning cell:
4385  // owning cell index, cgal vertices of cell
4386  using RequestType =
4387  std::vector<std::pair<unsigned int, IntersectionRequest>>;
4388  // Other ranks send back all found intersections for requesting cell:
4389  // requesting cell index, cgal vertices of found intersections
4390  using AnswerType =
4391  std::vector<std::pair<unsigned int, IntersectionAnswer>>;
4392 
4393  const auto create_request = [&](const unsigned int other_rank) {
4394  const auto other_rank_index = translate(other_rank);
4395 
4396  RequestType request;
4397  request.reserve(potential_owners_ptrs[other_rank_index + 1] -
4398  potential_owners_ptrs[other_rank_index]);
4399 
4400  for (unsigned int i = potential_owners_ptrs[other_rank_index];
4401  i < potential_owners_ptrs[other_rank_index + 1];
4402  ++i)
4403  request.emplace_back(
4404  potential_owners_indices[i],
4405  intersection_requests[potential_owners_indices[i]]);
4406 
4407  return request;
4408  };
4409 
4410 
4411  // TODO: this is potentially useful in many cases and it would be nice to
4412  // have cache.get_locally_owned_cell_bounding_boxes_rtree(marked_vertices)
4413  const auto construct_locally_owned_cell_bounding_boxes_rtree =
4414  [&cache](const std::vector<bool> &marked_verts) {
4415  const auto cell_marked = [&marked_verts](const auto &cell) {
4416  for (const unsigned int v : cell->vertex_indices())
4417  if (marked_verts[cell->vertex_index(v)])
4418  return true;
4419  return false;
4420  };
4421 
4422  const auto &boxes_and_cells =
4424 
4425  if (marked_verts.empty())
4426  return boxes_and_cells;
4427 
4428  std::vector<std::pair<
4431  potential_boxes_and_cells;
4432 
4433  for (const auto &box_and_cell : boxes_and_cells)
4434  if (cell_marked(box_and_cell.second))
4435  potential_boxes_and_cells.emplace_back(box_and_cell);
4436 
4437  return pack_rtree(potential_boxes_and_cells);
4438  };
4439 
4440 
4441  RTree<
4442  std::pair<BoundingBox<spacedim>,
4444  marked_cell_tree;
4445 
4446  const auto answer_request =
4447  [&](const unsigned int &other_rank,
4448  const RequestType &request) -> AnswerType {
4449  AnswerType answer;
4450 
4451  if (has_relevant_vertices)
4452  {
4453  if (marked_cell_tree.empty())
4454  {
4455  marked_cell_tree =
4456  construct_locally_owned_cell_bounding_boxes_rtree(
4457  marked_vertices);
4458  }
4459 
4460  // process requests
4461  for (unsigned int i = 0; i < request.size(); ++i)
4462  {
4463  // create a bounding box with tolerance
4464  const auto bb = BoundingBox<spacedim>(request[i].second)
4465  .create_extended(tolerance);
4466 
4467  for (const auto &box_cell :
4468  marked_cell_tree |
4469  boost::geometry::index::adaptors::queried(
4470  boost::geometry::index::intersects(bb)))
4471  {
4472 #ifdef DEAL_II_WITH_CGAL
4473  const auto &cell = box_cell.second;
4474  const auto &request_index = request[i].first;
4475  auto requested_intersection = request[i].second;
4476  CGALWrappers::resort_dealii_vertices_to_cgal_order(
4477  structdim, requested_intersection);
4478 
4479  const auto &try_intersection =
4480  CGALWrappers::get_vertices_in_cgal_order(
4481  cell, cache.get_mapping());
4482 
4483  const auto &found_intersections = CGALWrappers::
4484  compute_intersection_of_cells<dim, structdim, spacedim>(
4485  try_intersection, requested_intersection, tolerance);
4486 
4487  if (found_intersections.size() > 0)
4488  {
4489  for (const auto &found_intersection :
4490  found_intersections)
4491  {
4492  answer.emplace_back(request_index,
4493  found_intersection);
4494 
4495  send_components.emplace_back(
4496  std::make_pair(cell->level(), cell->index()),
4497  other_rank,
4498  request_index,
4499  found_intersection);
4500  }
4501  }
4502 #else
4503  (void)other_rank;
4504  (void)box_cell;
4505  Assert(false, ExcNeedsCGAL());
4506 #endif
4507  }
4508  }
4509  }
4510 
4511  return answer;
4512  };
4513 
4514  const auto process_answer = [&](const unsigned int other_rank,
4515  const AnswerType &answer) {
4516  for (unsigned int i = 0; i < answer.size(); ++i)
4517  recv_components.emplace_back(other_rank,
4518  answer[i].first,
4519  answer[i].second);
4520  };
4521 
4522  Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
4523  potential_owners_ranks,
4524  create_request,
4525  answer_request,
4526  process_answer,
4527  comm);
4528 
4529  // sort according to 1) intersection index and 2) rank (keeping the order
4530  // of recv components with same indices and ranks)
4531  std::stable_sort(recv_components.begin(),
4532  recv_components.end(),
4533  [&](const auto &a, const auto &b) {
4534  // intersection index
4535  if (std::get<1>(a) != std::get<1>(b))
4536  return std::get<1>(a) < std::get<1>(b);
4537 
4538  // rank
4539  return std::get<0>(a) < std::get<0>(b);
4540  });
4541 
4542  // sort according to 1) rank and 2) intersection index (keeping the
4543  // order of recv components with same indices and ranks)
4544  std::stable_sort(send_components.begin(),
4545  send_components.end(),
4546  [&](const auto &a, const auto &b) {
4547  // rank
4548  if (std::get<1>(a) != std::get<1>(b))
4549  return std::get<1>(a) < std::get<1>(b);
4550 
4551  // intersection idx
4552  return std::get<2>(a) < std::get<2>(b);
4553  });
4554 
4555  // construct recv_ptrs
4556  recv_ptrs.assign(intersection_requests.size() + 1, 0);
4557  for (const auto &rc : recv_components)
4558  ++recv_ptrs[std::get<1>(rc) + 1];
4559  for (unsigned int i = 0; i < intersection_requests.size(); ++i)
4560  recv_ptrs[i + 1] += recv_ptrs[i];
4561 
4562  return result;
4563  }
4564 
4565  } // namespace internal
4566 
4567 
4568 
4569  template <int spacedim>
4570  unsigned int
4571  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
4572  const Point<spacedim> &p)
4573  {
4574  auto id_and_v = std::min_element(
4575  vertices.begin(),
4576  vertices.end(),
4577  [&](const std::pair<const unsigned int, Point<spacedim>> &p1,
4578  const std::pair<const unsigned int, Point<spacedim>> &p2) -> bool {
4579  return p1.second.distance(p) < p2.second.distance(p);
4580  });
4581  return id_and_v->first;
4582  }
4583 
4584 
4585  template <int dim, int spacedim>
4586  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4587  Point<dim>>
4589  const Cache<dim, spacedim> &cache,
4590  const Point<spacedim> &p,
4592  &cell_hint,
4593  const std::vector<bool> &marked_vertices,
4594  const double tolerance)
4595  {
4596  const auto &mesh = cache.get_triangulation();
4597  const auto &mapping = cache.get_mapping();
4598  const auto &vertex_to_cells = cache.get_vertex_to_cell_map();
4599  const auto &vertex_to_cell_centers =
4601  const auto &used_vertices_rtree = cache.get_used_vertices_rtree();
4602 
4603  return find_active_cell_around_point(mapping,
4604  mesh,
4605  p,
4606  vertex_to_cells,
4607  vertex_to_cell_centers,
4608  cell_hint,
4609  marked_vertices,
4610  used_vertices_rtree,
4611  tolerance);
4612  }
4613 
4614  template <int spacedim>
4615  std::vector<std::vector<BoundingBox<spacedim>>>
4617  const std::vector<BoundingBox<spacedim>> &local_bboxes,
4618  const MPI_Comm mpi_communicator)
4619  {
4620 #ifndef DEAL_II_WITH_MPI
4621  (void)local_bboxes;
4622  (void)mpi_communicator;
4623  Assert(false,
4624  ExcMessage(
4625  "GridTools::exchange_local_bounding_boxes() requires MPI."));
4626  return {};
4627 #else
4628  // Step 1: preparing data to be sent
4629  unsigned int n_bboxes = local_bboxes.size();
4630  // Dimension of the array to be exchanged (number of double)
4631  int n_local_data = 2 * spacedim * n_bboxes;
4632  // data array stores each entry of each point describing the bounding
4633  // boxes
4634  std::vector<double> loc_data_array(n_local_data);
4635  for (unsigned int i = 0; i < n_bboxes; ++i)
4636  for (unsigned int d = 0; d < spacedim; ++d)
4637  {
4638  // Extracting the coordinates of each boundary point
4639  loc_data_array[2 * i * spacedim + d] =
4640  local_bboxes[i].get_boundary_points().first[d];
4641  loc_data_array[2 * i * spacedim + spacedim + d] =
4642  local_bboxes[i].get_boundary_points().second[d];
4643  }
4644 
4645  // Step 2: exchanging the size of local data
4646  unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_communicator);
4647 
4648  // Vector to store the size of loc_data_array for every process
4649  std::vector<int> size_all_data(n_procs);
4650 
4651  // Exchanging the number of bboxes
4652  int ierr = MPI_Allgather(&n_local_data,
4653  1,
4654  MPI_INT,
4655  size_all_data.data(),
4656  1,
4657  MPI_INT,
4658  mpi_communicator);
4659  AssertThrowMPI(ierr);
4660 
4661  // Now computing the displacement, relative to recvbuf,
4662  // at which to store the incoming data
4663  std::vector<int> rdispls(n_procs);
4664  rdispls[0] = 0;
4665  for (unsigned int i = 1; i < n_procs; ++i)
4666  rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
4667 
4668  // Step 3: exchange the data and bounding boxes:
4669  // Allocating a vector to contain all the received data
4670  std::vector<double> data_array(rdispls.back() + size_all_data.back());
4671 
4672  ierr = MPI_Allgatherv(loc_data_array.data(),
4673  n_local_data,
4674  MPI_DOUBLE,
4675  data_array.data(),
4676  size_all_data.data(),
4677  rdispls.data(),
4678  MPI_DOUBLE,
4679  mpi_communicator);
4680  AssertThrowMPI(ierr);
4681 
4682  // Step 4: create the array of bboxes for output
4683  std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
4684  unsigned int begin_idx = 0;
4685  for (unsigned int i = 0; i < n_procs; ++i)
4686  {
4687  // Number of local bounding boxes
4688  unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
4689  global_bboxes[i].resize(n_bbox_i);
4690  for (unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
4691  {
4692  Point<spacedim> p1, p2; // boundary points for bbox
4693  for (unsigned int d = 0; d < spacedim; ++d)
4694  {
4695  p1[d] = data_array[begin_idx + 2 * bbox * spacedim + d];
4696  p2[d] =
4697  data_array[begin_idx + 2 * bbox * spacedim + spacedim + d];
4698  }
4699  BoundingBox<spacedim> loc_bbox(std::make_pair(p1, p2));
4700  global_bboxes[i][bbox] = loc_bbox;
4701  }
4702  // Shifting the first index to the start of the next vector
4703  begin_idx += size_all_data[i];
4704  }
4705  return global_bboxes;
4706 #endif // DEAL_II_WITH_MPI
4707  }
4708 
4709 
4710 
4711  template <int spacedim>
4714  const std::vector<BoundingBox<spacedim>> &local_description,
4715  const MPI_Comm mpi_communicator)
4716  {
4717 #ifndef DEAL_II_WITH_MPI
4718  (void)mpi_communicator;
4719  // Building a tree with the only boxes available without MPI
4720  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> boxes_index(
4721  local_description.size());
4722  // Adding to each box the rank of the process owning it
4723  for (unsigned int i = 0; i < local_description.size(); ++i)
4724  boxes_index[i] = std::make_pair(local_description[i], 0u);
4725  return pack_rtree(boxes_index);
4726 #else
4727  // Exchanging local bounding boxes
4728  const std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes =
4729  Utilities::MPI::all_gather(mpi_communicator, local_description);
4730 
4731  // Preparing to flatten the vector
4732  const unsigned int n_procs =
4733  Utilities::MPI::n_mpi_processes(mpi_communicator);
4734  // The i'th element of the following vector contains the index of the
4735  // first local bounding box from the process of rank i
4736  std::vector<unsigned int> bboxes_position(n_procs);
4737 
4738  unsigned int tot_bboxes = 0;
4739  for (const auto &process_bboxes : global_bboxes)
4740  tot_bboxes += process_bboxes.size();
4741 
4742  // Now flattening the vector
4743  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
4744  flat_global_bboxes;
4745  flat_global_bboxes.reserve(tot_bboxes);
4746  unsigned int process_index = 0;
4747  for (const auto &process_bboxes : global_bboxes)
4748  {
4749  // Initialize a vector containing bounding boxes and rank of a process
4750  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
4751  boxes_and_indices(process_bboxes.size());
4752 
4753  // Adding to each box the rank of the process owning it
4754  for (unsigned int i = 0; i < process_bboxes.size(); ++i)
4755  boxes_and_indices[i] =
4756  std::make_pair(process_bboxes[i], process_index);
4757 
4758  flat_global_bboxes.insert(flat_global_bboxes.end(),
4759  boxes_and_indices.begin(),
4760  boxes_and_indices.end());
4761 
4762  ++process_index;
4763  }
4764 
4765  // Build a tree out of the bounding boxes. We avoid using the
4766  // insert method so that boost uses the packing algorithm
4767  return RTree<std::pair<BoundingBox<spacedim>, unsigned int>>(
4768  flat_global_bboxes.begin(), flat_global_bboxes.end());
4769 #endif // DEAL_II_WITH_MPI
4770  }
4771 
4772 
4773 
4774  template <int dim, int spacedim>
4775  void
4778  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
4779  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group)
4780  {
4781  // 1) determine for each vertex a vertex it coincides with and
4782  // put it into a map
4783  {
4784  static const int lookup_table_2d[2][2] =
4785  // flip:
4786  {
4787  {0, 1}, // false
4788  {1, 0} // true
4789  };
4790 
4791  static const int lookup_table_3d[2][2][2][4] =
4792  // orientation flip rotation
4793  {{{
4794  {0, 2, 1, 3}, // false false false
4795  {2, 3, 0, 1} // false false true
4796  },
4797  {
4798  {3, 1, 2, 0}, // false true false
4799  {1, 0, 3, 2} // false true true
4800  }},
4801  {{
4802  {0, 1, 2, 3}, // true false false
4803  {1, 3, 0, 2} // true false true
4804  },
4805  {
4806  {3, 2, 1, 0}, // true true false
4807  {2, 0, 3, 1} // true true true
4808  }}};
4809 
4810  // loop over all periodic face pairs
4811  for (const auto &pair : tria.get_periodic_face_map())
4812  {
4813  if (pair.first.first->level() != pair.second.first.first->level())
4814  continue;
4815 
4816  const auto face_a = pair.first.first->face(pair.first.second);
4817  const auto face_b =
4818  pair.second.first.first->face(pair.second.first.second);
4819  const auto mask = pair.second.second;
4820 
4821  AssertDimension(face_a->n_vertices(), face_b->n_vertices());
4822 
4823  // loop over all vertices on face
4824  for (unsigned int i = 0; i < face_a->n_vertices(); ++i)
4825  {
4826  const bool face_orientation = mask[0];
4827  const bool face_flip = mask[1];
4828  const bool face_rotation = mask[2];
4829 
4830  // find the right local vertex index for the second face
4831  unsigned int j = 0;
4832  switch (dim)
4833  {
4834  case 1:
4835  j = i;
4836  break;
4837  case 2:
4838  j = lookup_table_2d[face_flip][i];
4839  break;
4840  case 3:
4841  j = lookup_table_3d[face_orientation][face_flip]
4842  [face_rotation][i];
4843  break;
4844  default:
4845  AssertThrow(false, ExcNotImplemented());
4846  }
4847 
4848  // get vertex indices and store in map
4849  const auto vertex_a = face_a->vertex_index(i);
4850  const auto vertex_b = face_b->vertex_index(j);
4851  unsigned int temp = std::min(vertex_a, vertex_b);
4852 
4853  auto it_a = vertex_to_coinciding_vertex_group.find(vertex_a);
4854  if (it_a != vertex_to_coinciding_vertex_group.end())
4855  temp = std::min(temp, it_a->second);
4856 
4857  auto it_b = vertex_to_coinciding_vertex_group.find(vertex_b);
4858  if (it_b != vertex_to_coinciding_vertex_group.end())
4859  temp = std::min(temp, it_b->second);
4860 
4861  if (it_a != vertex_to_coinciding_vertex_group.end())
4862  it_a->second = temp;
4863  else
4864  vertex_to_coinciding_vertex_group[vertex_a] = temp;
4865 
4866  if (it_b != vertex_to_coinciding_vertex_group.end())
4867  it_b->second = temp;
4868  else
4869  vertex_to_coinciding_vertex_group[vertex_b] = temp;
4870  }
4871  }
4872 
4873  // 2) compress map: let vertices point to the coinciding vertex with
4874  // the smallest index
4875  for (auto &p : vertex_to_coinciding_vertex_group)
4876  {
4877  if (p.first == p.second)
4878  continue;
4879  unsigned int temp = p.second;
4880  while (temp != vertex_to_coinciding_vertex_group[temp])
4881  temp = vertex_to_coinciding_vertex_group[temp];
4882  p.second = temp;
4883  }
4884 
4885  // 3) create a map: smallest index of coinciding index -> all
4886  // coinciding indices
4887  for (auto p : vertex_to_coinciding_vertex_group)
4888  coinciding_vertex_groups[p.second] = {};
4889 
4890  for (auto p : vertex_to_coinciding_vertex_group)
4891  coinciding_vertex_groups[p.second].push_back(p.first);
4892  }
4893  }
4894 
4895 
4896 
4897  template <int dim, int spacedim>
4898  std::map<unsigned int, std::set<::types::subdomain_id>>
4901  {
4902  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4903  &tria) == nullptr) // nothing to do for a serial triangulation
4904  return {};
4905 
4906  // 1) collect for each vertex on periodic faces all vertices it coincides
4907  // with
4908  std::map<unsigned int, std::vector<unsigned int>> coinciding_vertex_groups;
4909  std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
4910 
4914 
4915  // 2) collect vertices belonging to local cells
4916  std::vector<bool> vertex_of_own_cell(tria.n_vertices(), false);
4917  for (const auto &cell :
4919  for (const unsigned int v : cell->vertex_indices())
4920  vertex_of_own_cell[cell->vertex_index(v)] = true;
4921 
4922  // 3) for each vertex belonging to a locally owned cell, find all ghost
4923  // neighbors (including the periodic own)
4924  std::map<unsigned int, std::set<types::subdomain_id>> result;
4925 
4926  // loop over all active ghost cells
4927  for (const auto &cell : tria.active_cell_iterators())
4928  if (cell->is_ghost())
4929  {
4930  const types::subdomain_id owner = cell->subdomain_id();
4931 
4932  // loop over all its vertices
4933  for (const unsigned int v : cell->vertex_indices())
4934  {
4935  // set owner if vertex belongs to a local cell
4936  if (vertex_of_own_cell[cell->vertex_index(v)])
4937  result[cell->vertex_index(v)].insert(owner);
4938 
4939  // mark also nodes coinciding due to periodicity
4940  auto coinciding_vertex_group =
4941  vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
4942  if (coinciding_vertex_group !=
4944  for (auto coinciding_vertex :
4945  coinciding_vertex_groups[coinciding_vertex_group->second])
4946  if (vertex_of_own_cell[coinciding_vertex])
4947  result[coinciding_vertex].insert(owner);
4948  }
4949  }
4950 
4951  return result;
4952  }
4953 
4954 
4955 
4956  namespace internal
4957  {
4958  template <int dim,
4959  unsigned int n_vertices,
4960  unsigned int n_sub_vertices,
4961  unsigned int n_configurations,
4962  unsigned int n_lines,
4963  unsigned int n_cols,
4964  typename value_type>
4965  void
4967  const std::array<unsigned int, n_configurations> &cut_line_table,
4969  const ndarray<unsigned int, n_lines, 2> &line_to_vertex_table,
4970  const std::vector<value_type> &ls_values,
4971  const std::vector<Point<dim>> &points,
4972  const std::vector<unsigned int> &mask,
4973  const double iso_level,
4974  const double tolerance,
4975  std::vector<Point<dim>> &vertices,
4976  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
4977  const bool write_back_cell_data)
4978  {
4979  // inspired by https://graphics.stanford.edu/~mdfisher/MarchingCubes.html
4980 
4981  constexpr unsigned int X = static_cast<unsigned int>(-1);
4982 
4983  // determine configuration
4984  unsigned int configuration = 0;
4985  for (unsigned int v = 0; v < n_vertices; ++v)
4986  if (ls_values[mask[v]] < iso_level)
4987  configuration |= (1 << v);
4988 
4989  // cell is not cut (nothing to do)
4990  if (cut_line_table[configuration] == 0)
4991  return;
4992 
4993  // helper function to determine where an edge (between index i and j) is
4994  // cut - see also: http://paulbourke.net/geometry/polygonise/
4995  const auto interpolate = [&](const unsigned int i, const unsigned int j) {
4996  if (std::abs(iso_level - ls_values[mask[i]]) < tolerance)
4997  return points[mask[i]];
4998  if (std::abs(iso_level - ls_values[mask[j]]) < tolerance)
4999  return points[mask[j]];
5000  if (std::abs(ls_values[mask[i]] - ls_values[mask[j]]) < tolerance)
5001  return points[mask[i]];
5002 
5003  const double mu = (iso_level - ls_values[mask[i]]) /
5004  (ls_values[mask[j]] - ls_values[mask[i]]);
5005 
5006  return Point<dim>(points[mask[i]] +
5007  mu * (points[mask[j]] - points[mask[i]]));
5008  };
5009 
5010  // determine the position where edges are cut (if they are cut)
5011  std::array<Point<dim>, n_lines> vertex_list_all;
5012  for (unsigned int l = 0; l < n_lines; ++l)
5013  if (cut_line_table[configuration] & (1 << l))
5014  vertex_list_all[l] =
5015  interpolate(line_to_vertex_table[l][0], line_to_vertex_table[l][1]);
5016 
5017  // merge duplicate vertices if possible
5018  unsigned int local_vertex_count = 0;
5019  std::array<Point<dim>, n_lines> vertex_list_reduced;
5020  std::array<unsigned int, n_lines> local_remap;
5021  std::fill(local_remap.begin(), local_remap.end(), X);
5022  for (int i = 0; new_line_table[configuration][i] != X; ++i)
5023  if (local_remap[new_line_table[configuration][i]] == X)
5024  {
5025  vertex_list_reduced[local_vertex_count] =
5026  vertex_list_all[new_line_table[configuration][i]];
5027  local_remap[new_line_table[configuration][i]] = local_vertex_count;
5028  local_vertex_count++;
5029  }
5030 
5031  // write back vertices
5032  const unsigned int n_vertices_old = vertices.size();
5033  for (unsigned int i = 0; i < local_vertex_count; ++i)
5034  vertices.push_back(vertex_list_reduced[i]);
5035 
5036  // write back cells
5037  if (write_back_cell_data && dim > 1)
5038  {
5039  for (unsigned int i = 0; new_line_table[configuration][i] != X;
5040  i += n_sub_vertices)
5041  {
5042  cells.resize(cells.size() + 1);
5043  cells.back().vertices.resize(n_sub_vertices);
5044 
5045  for (unsigned int v = 0; v < n_sub_vertices; ++v)
5046  cells.back().vertices[v] =
5047  local_remap[new_line_table[configuration][i + v]] +
5048  n_vertices_old;
5049  }
5050  }
5051  }
5052  } // namespace internal
5053 
5054 
5055 
5056  template <int dim, typename VectorType>
5058  const Mapping<dim, dim> &mapping,
5059  const FiniteElement<dim, dim> &fe,
5060  const unsigned int n_subdivisions,
5061  const double tolerance)
5062  : n_subdivisions(n_subdivisions)
5063  , tolerance(tolerance)
5064  , fe_values(mapping,
5065  fe,
5066  create_quadrature_rule(n_subdivisions),
5068  {}
5069 
5070 
5071 
5072  template <int dim, typename VectorType>
5075  const unsigned int n_subdivisions)
5076  {
5077  std::vector<Point<dim>> quadrature_points;
5078 
5079  if (dim == 1)
5080  {
5081  for (unsigned int i = 0; i <= n_subdivisions; ++i)
5082  quadrature_points.emplace_back(1.0 / n_subdivisions * i);
5083  }
5084  else if (dim == 2)
5085  {
5086  for (unsigned int j = 0; j <= n_subdivisions; ++j)
5087  for (unsigned int i = 0; i <= n_subdivisions; ++i)
5088  quadrature_points.emplace_back(1.0 / n_subdivisions * i,
5089  1.0 / n_subdivisions * j);
5090  }
5091  else
5092  {
5093  for (unsigned int k = 0; k <= n_subdivisions; ++k)
5094  for (unsigned int j = 0; j <= n_subdivisions; ++j)
5095  for (unsigned int i = 0; i <= n_subdivisions; ++i)
5096  quadrature_points.emplace_back(1.0 / n_subdivisions * i,
5097  1.0 / n_subdivisions * j,
5098  1.0 / n_subdivisions * k);
5099  }
5100 
5101 
5102  return {quadrature_points};
5103  }
5104 
5105 
5106 
5107  template <int dim, typename VectorType>
5108  void
5110  const DoFHandler<dim> &background_dof_handler,
5111  const VectorType &ls_vector,
5112  const double iso_level,
5113  std::vector<Point<dim>> &vertices,
5114  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const
5115  {
5116  AssertThrow(
5117  dim > 1,
5118  ExcMessage(
5119  "Not implemented for dim==1. Use the alternative process()-function "
5120  "not returning a vector of CellData objects."));
5121 
5122  for (const auto &cell : background_dof_handler.active_cell_iterators() |
5124  process_cell(cell, ls_vector, iso_level, vertices, cells);
5125  }
5126 
5127  template <int dim, typename VectorType>
5128  void
5130  const DoFHandler<dim> &background_dof_handler,
5131  const VectorType &ls_vector,
5132  const double iso_level,
5133  std::vector<Point<dim>> &vertices) const
5134  {
5135  for (const auto &cell : background_dof_handler.active_cell_iterators() |
5137  process_cell(cell, ls_vector, iso_level, vertices);
5138 
5139  delete_duplicated_vertices(vertices, 1e-10 /*tol*/);
5140  }
5141 
5142 
5143  template <int dim, typename VectorType>
5144  void
5146  const typename DoFHandler<dim>::active_cell_iterator &cell,
5147  const VectorType &ls_vector,
5148  const double iso_level,
5149  std::vector<Point<dim>> &vertices,
5150  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const
5151  {
5152  AssertThrow(
5153  dim > 1,
5154  ExcMessage(
5155  "Not implemented for dim==1. Use the alternative process_cell()-function "
5156  "not returning a vector of CellData objects."));
5157 
5158  std::vector<value_type> ls_values;
5159 
5160  fe_values.reinit(cell);
5161  ls_values.resize(fe_values.n_quadrature_points);
5162  fe_values.get_function_values(ls_vector, ls_values);
5163  process_cell(
5164  ls_values, fe_values.get_quadrature_points(), iso_level, vertices, cells);
5165  }
5166 
5167  template <int dim, typename VectorType>
5168  void
5170  const typename DoFHandler<dim>::active_cell_iterator &cell,
5171  const VectorType &ls_vector,
5172  const double iso_level,
5173  std::vector<Point<dim>> &vertices) const
5174  {
5175  // This vector is just a placeholder to reuse the process_cell function.
5176  std::vector<CellData<dim == 1 ? 1 : dim - 1>> dummy_cells;
5177 
5178  std::vector<value_type> ls_values;
5179 
5180  fe_values.reinit(cell);
5181  ls_values.resize(fe_values.n_quadrature_points);
5182  fe_values.get_function_values(ls_vector, ls_values);
5183 
5184  process_cell(ls_values,
5185  fe_values.get_quadrature_points(),
5186  iso_level,
5187  vertices,
5188  dummy_cells,
5189  false /*don't write back cell data*/);
5190  }
5191 
5192 
5193  template <int dim, typename VectorType>
5194  void
5196  std::vector<value_type> &ls_values,
5197  const std::vector<Point<dim>> &points,
5198  const double iso_level,
5199  std::vector<Point<dim>> &vertices,
5200  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
5201  const bool write_back_cell_data) const
5202  {
5203  const unsigned p = n_subdivisions + 1;
5204 
5205  if (dim == 1)
5206  {
5207  for (unsigned int i = 0; i < n_subdivisions; ++i)
5208  {
5209  std::vector<unsigned int> mask{i + 0, i + 1};
5210 
5211  // check if a corner node is cut
5212  if (std::abs(iso_level - ls_values[mask[0]]) < tolerance)
5213  vertices.emplace_back(points[mask[0]]);
5214  else if (std::abs(iso_level - ls_values[mask[1]]) < tolerance)
5215  {
5216  if (i + 1 == n_subdivisions)
5217  vertices.emplace_back(points[mask[1]]);
5218  }
5219  // check if the edge is cut
5220  else if (((ls_values[mask[0]] > iso_level) &&
5221  (ls_values[mask[1]] < iso_level)) ||
5222  ((ls_values[mask[0]] < iso_level) &&
5223  (ls_values[mask[1]] > iso_level)))
5224  {
5225  // determine the interpolation weight (0<mu<1)
5226  const double mu = (iso_level - ls_values[mask[0]]) /
5227  (ls_values[mask[1]] - ls_values[mask[0]]);
5228 
5229  // interpolate
5230  vertices.emplace_back(points[mask[0]] +
5231  mu * (points[mask[1]] - points[mask[0]]));
5232  }
5233  }
5234  }
5235  else if (dim == 2)
5236  {
5237  for (unsigned int j = 0; j < n_subdivisions; ++j)
5238  for (unsigned int i = 0; i < n_subdivisions; ++i)
5239  {
5240  std::vector<unsigned int> mask{p * (j + 0) + (i + 0),
5241  p * (j + 0) + (i + 1),
5242  p * (j + 1) + (i + 1),
5243  p * (j + 1) + (i + 0)};
5244 
5245  process_sub_cell(ls_values,
5246  points,
5247  mask,
5248  iso_level,
5249  vertices,
5250  cells,
5251  write_back_cell_data);
5252  }
5253  }
5254  else if (dim == 3)
5255  {
5256  for (unsigned int k = 0; k < n_subdivisions; ++k)
5257  for (unsigned int j = 0; j < n_subdivisions; ++j)
5258  for (unsigned int i = 0; i < n_subdivisions; ++i)
5259  {
5260  std::vector<unsigned int> mask{
5261  p * p * (k + 0) + p * (j + 0) + (i + 0),
5262  p * p * (k + 0) + p * (j + 0) + (i + 1),
5263  p * p * (k + 0) + p * (j + 1) + (i + 1),
5264  p * p * (k + 0) + p * (j + 1) + (i + 0),
5265  p * p * (k + 1) + p * (j + 0) + (i + 0),
5266  p * p * (k + 1) + p * (j + 0) + (i + 1),
5267  p * p * (k + 1) + p * (j + 1) + (i + 1),
5268  p * p * (k + 1) + p * (j + 1) + (i + 0)};
5269 
5270  process_sub_cell(ls_values,
5271  points,
5272  mask,
5273  iso_level,
5274  vertices,
5275  cells,
5276  write_back_cell_data);
5277  }
5278  }
5279  }
5280 
5281 
5282 
5283  template <int dim, typename VectorType>
5284  void
5286  const std::vector<value_type> &ls_values,
5287  const std::vector<Point<2>> &points,
5288  const std::vector<unsigned int> &mask,
5289  const double iso_level,
5290  std::vector<Point<2>> &vertices,
5291  std::vector<CellData<1>> &cells,
5292  const bool write_back_cell_data) const
5293  {
5294  // set up dimension-dependent sizes and tables
5295  constexpr unsigned int n_vertices = 4;
5296  constexpr unsigned int n_sub_vertices = 2;
5297  constexpr unsigned int n_lines = 4;
5298  constexpr unsigned int n_configurations = Utilities::pow(2, n_vertices);
5299  constexpr unsigned int X = static_cast<unsigned int>(-1);
5300 
5301  // table that indicates if an edge is cut (if the i-th bit is set the i-th
5302  // line is cut)
5303  constexpr std::array<unsigned int, n_configurations> cut_line_table = {
5304  {0b0000,
5305  0b0101,
5306  0b0110,
5307  0b0011,
5308  0b1010,
5309  0b0000,
5310  0b1100,
5311  0b1001,
5312  0b1001,
5313  0b1100,
5314  0b0000,
5315  0b1010,
5316  0b0011,
5317  0b0110,
5318  0b0101,
5319  0b0000}};
5320 
5321  // list of the definition of the newly created lines (each line is defined
5322  // by two edges it cuts)
5323  constexpr ndarray<unsigned int, n_configurations, 5> new_line_table = {
5324  {{{X, X, X, X, X}},
5325  {{0, 2, X, X, X}},
5326  {{1, 2, X, X, X}},
5327  {{0, 1, X, X, X}},
5328  {{1, 3, X, X, X}},
5329  {{X, X, X, X, X}},
5330  {{2, 3, X, X, X}},
5331  {{0, 3, X, X, X}},
5332  {{0, 3, X, X, X}},
5333  {{2, 3, X, X, X}},
5334  {{X, X, X, X, X}},
5335  {{1, 3, X, X, X}},
5336  {{0, 1, X, X, X}},
5337  {{2, 1, X, X, X}},
5338  {{0, 2, X, X, X}},
5339  {{X, X, X, X, X}}}};
5340 
5341  // vertices of each line
5342  constexpr ndarray<unsigned int, n_lines, 2> line_to_vertex_table = {
5343  {{{0, 3}}, {{1, 2}}, {{0, 1}}, {{3, 2}}}};
5344 
5345  // run dimension-independent code
5347  n_vertices,
5348  n_sub_vertices,
5349  n_configurations,
5350  n_lines,
5351  5>(cut_line_table,
5352  new_line_table,
5353  line_to_vertex_table,
5354  ls_values,
5355  points,
5356  mask,
5357  iso_level,
5358  tolerance,
5359  vertices,
5360  cells,
5361  write_back_cell_data);
5362  }
5363 
5364 
5365 
5366  template <int dim, typename VectorType>
5367  void
5369  const std::vector<value_type> &ls_values,
5370  const std::vector<Point<3>> &points,
5371  const std::vector<unsigned int> &mask,
5372  const double iso_level,
5373  std::vector<Point<3>> &vertices,
5374  std::vector<CellData<2>> &cells,
5375  const bool write_back_cell_data) const
5376  {
5377  // set up dimension-dependent sizes and tables
5378  constexpr unsigned int n_vertices = 8;
5379  constexpr unsigned int n_sub_vertices = 3;
5380  constexpr unsigned int n_lines = 12;
5381  constexpr unsigned int n_configurations = Utilities::pow(2, n_vertices);
5382  constexpr unsigned int X = static_cast<unsigned int>(-1);
5383 
5384  // clang-format off
5385  // table that indicates if an edge is cut (if the i-th bit is set the i-th
5386  // line is cut)
5387  constexpr std::array<unsigned int, n_configurations> cut_line_table = {{
5388  0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905,
5389  0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 0x190, 0x99, 0x393, 0x29a,
5390  0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93,
5391  0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, 0x53c,
5392  0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9,
5393  0x1a3, 0xaa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6,
5394  0xfaa, 0xea3, 0xda9, 0xca0, 0x460, 0x569, 0x663, 0x76a, 0x66, 0x16f,
5395  0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
5396  0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5,
5397  0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a,
5398  0x256, 0x35f, 0x55, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53,
5399  0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc,
5400  0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9,
5401  0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6,
5402  0x4ca, 0x5c3, 0x6c9, 0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f,
5403  0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
5404  0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5,
5405  0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a,
5406  0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663,
5407  0x569, 0x460, 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
5408  0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, 0xd30, 0xc39,
5409  0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636,
5410  0x13a, 0x33, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f,
5411  0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190,
5412  0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605,
5413  0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0}};
5414  // clang-format on
5415 
5416  // list of the definition of the newly created triangles (each triangles is
5417  // defined by two edges it cuts)
5418  constexpr ndarray<unsigned int, n_configurations, 16> new_line_table = {
5419  {{{X, X, X, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5420  {{0, 8, 3, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5421  {{0, 1, 9, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5422  {{1, 8, 3, 9, 8, 1, X, X, X, X, X, X, X, X, X, X}},
5423  {{1, 2, 10, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5424  {{0, 8, 3, 1, 2, 10, X, X, X, X, X, X, X, X, X, X}},
5425  {{9, 2, 10, 0, 2, 9, X, X, X, X, X, X, X, X, X, X}},
5426  {{2, 8, 3, 2, 10, 8, 10, 9, 8, X, X, X, X, X, X, X}},
5427  {{3, 11, 2, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5428  {{0, 11, 2, 8, 11, 0, X, X, X, X, X, X, X, X, X, X}},
5429  {{1, 9, 0, 2, 3, 11, X, X, X, X, X, X, X, X, X, X}},
5430  {{1, 11, 2, 1, 9, 11, 9, 8, 11, X, X, X, X, X, X, X}},
5431  {{3, 10, 1, 11, 10, 3, X, X, X, X, X, X, X, X, X, X}},
5432  {{0, 10, 1, 0, 8, 10, 8, 11, 10, X, X, X, X, X, X, X}},
5433  {{3, 9, 0, 3, 11, 9, 11, 10, 9, X, X, X, X, X, X, X}},
5434  {{9, 8, 10, 10, 8, 11, X, X, X, X, X, X, X, X, X, X}},
5435  {{4, 7, 8, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5436  {{4, 3, 0, 7, 3, 4, X, X, X, X, X, X, X, X, X, X}},
5437  {{0, 1, 9, 8, 4, 7, X, X, X, X, X, X, X, X, X, X}},
5438  {{4, 1, 9, 4, 7, 1, 7, 3, 1, X, X, X, X, X, X, X}},
5439  {{1, 2, 10, 8, 4, 7, X, X, X, X, X, X, X, X, X, X}},
5440  {{3, 4, 7, 3, 0, 4, 1, 2, 10, X, X, X, X, X, X, X}},
5441  {{9, 2, 10, 9, 0, 2, 8, 4, 7, X, X, X, X, X, X, X}},
5442  {{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, X, X, X, X}},
5443  {{8, 4, 7, 3, 11, 2, X, X, X, X, X, X, X, X, X, X}},
5444  {{11, 4, 7, 11, 2, 4, 2, 0, 4, X, X, X, X, X, X, X}},
5445  {{9, 0, 1, 8, 4, 7, 2, 3, 11, X, X, X, X, X, X, X}},
5446  {{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, X, X, X, X}},
5447  {{3, 10, 1, 3, 11, 10, 7, 8, 4, X, X, X, X, X, X, X}},
5448  {{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, X, X, X, X}},
5449  {{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, X, X, X, X}},
5450  {{4, 7, 11, 4, 11, 9, 9, 11, 10, X, X, X, X, X, X, X}},
5451  {{9, 5, 4, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5452  {{9, 5, 4, 0, 8, 3, X, X, X, X, X, X, X, X, X, X}},
5453  {{0, 5, 4, 1, 5, 0, X, X, X, X, X, X, X, X, X, X}},
5454  {{8, 5, 4, 8, 3, 5, 3, 1, 5, X, X, X, X, X, X, X}},
5455  {{1, 2, 10, 9, 5, 4, X, X, X, X, X, X, X, X, X, X}},
5456  {{3, 0, 8, 1, 2, 10, 4, 9, 5, X, X, X, X, X, X, X}},
5457  {{5, 2, 10, 5, 4, 2, 4, 0, 2, X, X, X, X, X, X, X}},
5458  {{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, X, X, X, X}},
5459  {{9, 5, 4, 2, 3, 11, X, X, X, X, X, X, X, X, X, X}},
5460  {{0, 11, 2, 0, 8, 11, 4, 9, 5, X, X, X, X, X, X, X}},
5461  {{0, 5, 4, 0, 1, 5, 2, 3, 11, X, X, X, X, X, X, X}},
5462  {{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, X, X, X, X}},
5463  {{10, 3, 11, 10, 1, 3, 9, 5, 4, X, X, X, X, X, X, X}},
5464  {{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, X, X, X, X}},
5465  {{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, X, X, X, X}},
5466  {{5, 4, 8, 5, 8, 10, 10, 8, 11, X, X, X, X, X, X, X}},
5467  {{9, 7, 8, 5, 7, 9, X, X, X, X, X, X, X, X, X, X}},
5468  {{9, 3, 0, 9, 5, 3, 5, 7, 3, X, X, X, X, X, X, X}},
5469  {{0, 7, 8, 0, 1, 7, 1, 5, 7, X, X, X, X, X, X, X}},
5470  {{1, 5, 3, 3, 5, 7, X, X, X, X, X, X, X, X, X, X}},
5471  {{9, 7, 8, 9, 5, 7, 10, 1, 2, X, X, X, X, X, X, X}},
5472  {{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, X, X, X, X}},
5473  {{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, X, X, X, X}},
5474  {{2, 10, 5, 2, 5, 3, 3, 5, 7, X, X, X, X, X, X, X}},
5475  {{7, 9, 5, 7, 8, 9, 3, 11, 2, X, X, X, X, X, X, X}},
5476  {{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, X, X, X, X}},
5477  {{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, X, X, X, X}},
5478  {{11, 2, 1, 11, 1, 7, 7, 1, 5, X, X, X, X, X, X, X}},
5479  {{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, X, X, X, X}},
5480  {{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, X}},
5481  {{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, X}},
5482  {{11, 10, 5, 7, 11, 5, X, X, X, X, X, X, X, X, X, X}},
5483  {{10, 6, 5, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5484  {{0, 8, 3, 5, 10, 6, X, X, X, X, X, X, X, X, X, X}},
5485  {{9, 0, 1, 5, 10, 6, X, X, X, X, X, X, X, X, X, X}},
5486  {{1, 8, 3, 1, 9, 8, 5, 10, 6, X, X, X, X, X, X, X}},
5487  {{1, 6, 5, 2, 6, 1, X, X, X, X, X, X, X, X, X, X}},
5488  {{1, 6, 5, 1, 2, 6, 3, 0, 8, X, X, X, X, X, X, X}},
5489  {{9, 6, 5, 9, 0, 6, 0, 2, 6, X, X, X, X, X, X, X}},
5490  {{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, X, X, X, X}},
5491  {{2, 3, 11, 10, 6, 5, X, X, X, X, X, X, X, X, X, X}},
5492  {{11, 0, 8, 11, 2, 0, 10, 6, 5, X, X, X, X, X, X, X}},
5493  {{0, 1, 9, 2, 3, 11, 5, 10, 6, X, X, X, X, X, X, X}},
5494  {{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, X, X, X, X}},
5495  {{6, 3, 11, 6, 5, 3, 5, 1, 3, X, X, X, X, X, X, X}},
5496  {{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, X, X, X, X}},
5497  {{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, X, X, X, X}},
5498  {{6, 5, 9, 6, 9, 11, 11, 9, 8, X, X, X, X, X, X, X}},
5499  {{5, 10, 6, 4, 7, 8, X, X, X, X, X, X, X, X, X, X}},
5500  {{4, 3, 0, 4, 7, 3, 6, 5, 10, X, X, X, X, X, X, X}},
5501  {{1, 9, 0, 5, 10, 6, 8, 4, 7, X, X, X, X, X, X, X}},
5502  {{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, X, X, X, X}},
5503  {{6, 1, 2, 6, 5, 1, 4, 7, 8, X, X, X, X, X, X, X}},
5504  {{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, X, X, X, X}},
5505  {{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, X, X, X, X}},
5506  {{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, X}},
5507  {{3, 11, 2, 7, 8, 4, 10, 6, 5, X, X, X, X, X, X, X}},
5508  {{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, X, X, X, X}},
5509  {{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, X, X, X, X}},
5510  {{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, X}},
5511  {{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, X, X, X, X}},
5512  {{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, X}},
5513  {{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, X}},
5514  {{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, X, X, X, X}},
5515  {{10, 4, 9, 6, 4, 10, X, X, X, X, X, X, X, X, X, X}},
5516  {{4, 10, 6, 4, 9, 10, 0, 8, 3, X, X, X, X, X, X, X}},
5517  {{10, 0, 1, 10, 6, 0, 6, 4, 0, X, X, X, X, X, X, X}},
5518  {{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, X, X, X, X}},
5519  {{1, 4, 9, 1, 2, 4, 2, 6, 4, X, X, X, X, X, X, X}},
5520  {{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, X, X, X, X}},
5521  {{0, 2, 4, 4, 2, 6, X, X, X, X, X, X, X, X, X, X}},
5522  {{8, 3, 2, 8, 2, 4, 4, 2, 6, X, X, X, X, X, X, X}},
5523  {{10, 4, 9, 10, 6, 4, 11, 2, 3, X, X, X, X, X, X, X}},
5524  {{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, X, X, X, X}},
5525  {{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, X, X, X, X}},
5526  {{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, X}},
5527  {{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, X, X, X, X}},
5528  {{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, X}},
5529  {{3, 11, 6, 3, 6, 0, 0, 6, 4, X, X, X, X, X, X, X}},
5530  {{6, 4, 8, 11, 6, 8, X, X, X, X, X, X, X, X, X, X}},
5531  {{7, 10, 6, 7, 8, 10, 8, 9, 10, X, X, X, X, X, X, X}},
5532  {{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, X, X, X, X}},
5533  {{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, X, X, X, X}},
5534  {{10, 6, 7, 10, 7, 1, 1, 7, 3, X, X, X, X, X, X, X}},
5535  {{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, X, X, X, X}},
5536  {{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, X}},
5537  {{7, 8, 0, 7, 0, 6, 6, 0, 2, X, X, X, X, X, X, X}},
5538  {{7, 3, 2, 6, 7, 2, X, X, X, X, X, X, X, X, X, X}},
5539  {{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, X, X, X, X}},
5540  {{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, X}},
5541  {{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, X}},
5542  {{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, X, X, X, X}},
5543  {{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, X}},
5544  {{0, 9, 1, 11, 6, 7, X, X, X, X, X, X, X, X, X, X}},
5545  {{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, X, X, X, X}},
5546  {{7, 11, 6, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5547  {{7, 6, 11, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5548  {{3, 0, 8, 11, 7, 6, X, X, X, X, X, X, X, X, X, X}},
5549  {{0, 1, 9, 11, 7, 6, X, X, X, X, X, X, X, X, X, X}},
5550  {{8, 1, 9, 8, 3, 1, 11, 7, 6, X, X, X, X, X, X, X}},
5551  {{10, 1, 2, 6, 11, 7, X, X, X, X, X, X, X, X, X, X}},
5552  {{1, 2, 10, 3, 0, 8, 6, 11, 7, X, X, X, X, X, X, X}},
5553  {{2, 9, 0, 2, 10, 9, 6, 11, 7, X, X, X, X, X, X, X}},
5554  {{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, X, X, X, X}},
5555  {{7, 2, 3, 6, 2, 7, X, X, X, X, X, X, X, X, X, X}},
5556  {{7, 0, 8, 7, 6, 0, 6, 2, 0, X, X, X, X, X, X, X}},
5557  {{2, 7, 6, 2, 3, 7, 0, 1, 9, X, X, X, X, X, X, X}},
5558  {{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, X, X, X, X}},
5559  {{10, 7, 6, 10, 1, 7, 1, 3, 7, X, X, X, X, X, X, X}},
5560  {{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, X, X, X, X}},
5561  {{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, X, X, X, X}},
5562  {{7, 6, 10, 7, 10, 8, 8, 10, 9, X, X, X, X, X, X, X}},
5563  {{6, 8, 4, 11, 8, 6, X, X, X, X, X, X, X, X, X, X}},
5564  {{3, 6, 11, 3, 0, 6, 0, 4, 6, X, X, X, X, X, X, X}},
5565  {{8, 6, 11, 8, 4, 6, 9, 0, 1, X, X, X, X, X, X, X}},
5566  {{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, X, X, X, X}},
5567  {{6, 8, 4, 6, 11, 8, 2, 10, 1, X, X, X, X, X, X, X}},
5568  {{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, X, X, X, X}},
5569  {{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, X, X, X, X}},
5570  {{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, X}},
5571  {{8, 2, 3, 8, 4, 2, 4, 6, 2, X, X, X, X, X, X, X}},
5572  {{0, 4, 2, 4, 6, 2, X, X, X, X, X, X, X, X, X, X}},
5573  {{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, X, X, X, X}},
5574  {{1, 9, 4, 1, 4, 2, 2, 4, 6, X, X, X, X, X, X, X}},
5575  {{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, X, X, X, X}},
5576  {{10, 1, 0, 10, 0, 6, 6, 0, 4, X, X, X, X, X, X, X}},
5577  {{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, X}},
5578  {{10, 9, 4, 6, 10, 4, X, X, X, X, X, X, X, X, X, X}},
5579  {{4, 9, 5, 7, 6, 11, X, X, X, X, X, X, X, X, X, X}},
5580  {{0, 8, 3, 4, 9, 5, 11, 7, 6, X, X, X, X, X, X, X}},
5581  {{5, 0, 1, 5, 4, 0, 7, 6, 11, X, X, X, X, X, X, X}},
5582  {{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, X, X, X, X}},
5583  {{9, 5, 4, 10, 1, 2, 7, 6, 11, X, X, X, X, X, X, X}},
5584  {{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, X, X, X, X}},
5585  {{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, X, X, X, X}},
5586  {{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, X}},
5587  {{7, 2, 3, 7, 6, 2, 5, 4, 9, X, X, X, X, X, X, X}},
5588  {{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, X, X, X, X}},
5589  {{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, X, X, X, X}},
5590  {{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, X}},
5591  {{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, X, X, X, X}},
5592  {{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, X}},
5593  {{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, X}},
5594  {{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, X, X, X, X}},
5595  {{6, 9, 5, 6, 11, 9, 11, 8, 9, X, X, X, X, X, X, X}},
5596  {{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, X, X, X, X}},
5597  {{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, X, X, X, X}},
5598  {{6, 11, 3, 6, 3, 5, 5, 3, 1, X, X, X, X, X, X, X}},
5599  {{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, X, X, X, X}},
5600  {{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, X}},
5601  {{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, X}},
5602  {{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, X, X, X, X}},
5603  {{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, X, X, X, X}},
5604  {{9, 5, 6, 9, 6, 0, 0, 6, 2, X, X, X, X, X, X, X}},
5605  {{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, X}},
5606  {{1, 5, 6, 2, 1, 6, X, X, X, X, X, X, X, X, X, X}},
5607  {{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, X}},
5608  {{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, X, X, X, X}},
5609  {{0, 3, 8, 5, 6, 10, X, X, X, X, X, X, X, X, X, X}},
5610  {{10, 5, 6, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5611  {{11, 5, 10, 7, 5, 11, X, X, X, X, X, X, X, X, X, X}},
5612  {{11, 5, 10, 11, 7, 5, 8, 3, 0, X, X, X, X, X, X, X}},
5613  {{5, 11, 7, 5, 10, 11, 1, 9, 0, X, X, X, X, X, X, X}},
5614  {{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, X, X, X, X}},
5615  {{11, 1, 2, 11, 7, 1, 7, 5, 1, X, X, X, X, X, X, X}},
5616  {{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, X, X, X, X}},
5617  {{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, X, X, X, X}},
5618  {{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, X}},
5619  {{2, 5, 10, 2, 3, 5, 3, 7, 5, X, X, X, X, X, X, X}},
5620  {{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, X, X, X, X}},
5621  {{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, X, X, X, X}},
5622  {{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, X}},
5623  {{1, 3, 5, 3, 7, 5, X, X, X, X, X, X, X, X, X, X}},
5624  {{0, 8, 7, 0, 7, 1, 1, 7, 5, X, X, X, X, X, X, X}},
5625  {{9, 0, 3, 9, 3, 5, 5, 3, 7, X, X, X, X, X, X, X}},
5626  {{9, 8, 7, 5, 9, 7, X, X, X, X, X, X, X, X, X, X}},
5627  {{5, 8, 4, 5, 10, 8, 10, 11, 8, X, X, X, X, X, X, X}},
5628  {{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, X, X, X, X}},
5629  {{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, X, X, X, X}},
5630  {{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, X}},
5631  {{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, X, X, X, X}},
5632  {{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, X}},
5633  {{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, X}},
5634  {{9, 4, 5, 2, 11, 3, X, X, X, X, X, X, X, X, X, X}},
5635  {{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, X, X, X, X}},
5636  {{5, 10, 2, 5, 2, 4, 4, 2, 0, X, X, X, X, X, X, X}},
5637  {{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, X}},
5638  {{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, X, X, X, X}},
5639  {{8, 4, 5, 8, 5, 3, 3, 5, 1, X, X, X, X, X, X, X}},
5640  {{0, 4, 5, 1, 0, 5, X, X, X, X, X, X, X, X, X, X}},
5641  {{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, X, X, X, X}},
5642  {{9, 4, 5, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5643  {{4, 11, 7, 4, 9, 11, 9, 10, 11, X, X, X, X, X, X, X}},
5644  {{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, X, X, X, X}},
5645  {{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, X, X, X, X}},
5646  {{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, X}},
5647  {{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, X, X, X, X}},
5648  {{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, X}},
5649  {{11, 7, 4, 11, 4, 2, 2, 4, 0, X, X, X, X, X, X, X}},
5650  {{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, X, X, X, X}},
5651  {{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, X, X, X, X}},
5652  {{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, X}},
5653  {{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, X}},
5654  {{1, 10, 2, 8, 7, 4, X, X, X, X, X, X, X, X, X, X}},
5655  {{4, 9, 1, 4, 1, 7, 7, 1, 3, X, X, X, X, X, X, X}},
5656  {{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, X, X, X, X}},
5657  {{4, 0, 3, 7, 4, 3, X, X, X, X, X, X, X, X, X, X}},
5658  {{4, 8, 7, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5659  {{9, 10, 8, 10, 11, 8, X, X, X, X, X, X, X, X, X, X}},
5660  {{3, 0, 9, 3, 9, 11, 11, 9, 10, X, X, X, X, X, X, X}},
5661  {{0, 1, 10, 0, 10, 8, 8, 10, 11, X, X, X, X, X, X, X}},
5662  {{3, 1, 10, 11, 3, 10, X, X, X, X, X, X, X, X, X, X}},
5663  {{1, 2, 11, 1, 11, 9, 9, 11, 8, X, X, X, X, X, X, X}},
5664  {{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, X, X, X, X}},
5665  {{0, 2, 11, 8, 0, 11, X, X, X, X, X, X, X, X, X, X}},
5666  {{3, 2, 11, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5667  {{2, 3, 8, 2, 8, 10, 10, 8, 9, X, X, X, X, X, X, X}},
5668  {{9, 10, 2, 0, 9, 2, X, X, X, X, X, X, X, X, X, X}},
5669  {{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, X, X, X, X}},
5670  {{1, 10, 2, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5671  {{1, 3, 8, 9, 1, 8, X, X, X, X, X, X, X, X, X, X}},
5672  {{0, 9, 1, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5673  {{0, 3, 8, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5674  {{X, X, X, X, X, X, X, X, X, X, X, X, X, X, X, X}}}};
5675 
5676  // vertices of each line
5677  static constexpr ndarray<unsigned int, n_lines, 2> line_to_vertex_table = {
5678  {{{0, 1}},
5679  {{1, 2}},
5680  {{2, 3}},
5681  {{3, 0}},
5682  {{4, 5}},
5683  {{5, 6}},
5684  {{6, 7}},
5685  {{7, 4}},
5686  {{0, 4}},
5687  {{1, 5}},
5688  {{2, 6}},
5689  {{3, 7}}}};
5690 
5691  // run dimension-independent code
5693  n_vertices,
5694  n_sub_vertices,
5695  n_configurations,
5696  n_lines,
5697  16>(cut_line_table,
5698  new_line_table,
5699  line_to_vertex_table,
5700  ls_values,
5701  points,
5702  mask,
5703  iso_level,
5704  tolerance,
5705  vertices,
5706  cells,
5707  write_back_cell_data);
5708  }
5709 
5710 } /* namespace GridTools */
5711 
5712 
5713 // explicit instantiations
5714 #include "grid_tools.inst"
5715 
void distribute(VectorType &vec) const
std::pair< std::vector< std::pair< int, int > >, std::vector< int > > query(const QueryType &queries)
DistributedTree(const MPI_Comm comm, const std::vector< BoundingBox< dim, Number >> &bounding_boxes)
BoundingBox< spacedim, Number > create_extended(const Number amount) const
BoundingBox< spacedim, Number > create_extended_relative(const Number relative_amount) const
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
types::global_dof_index n_dofs() const
Definition: fe_q.h:551
const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_vertex_to_cell_map() const
const std::vector< std::vector< Tensor< 1, spacedim > > > & get_vertex_to_cell_centers_directions() const
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_locally_owned_cell_bounding_boxes_rtree() const
const Mapping< dim, spacedim > & get_mapping() const
const Triangulation< dim, spacedim > & get_triangulation() const
const RTree< std::pair< Point< spacedim >, unsigned int > > & get_used_vertices_rtree() const
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_cell_bounding_boxes_rtree() const
void process(const DoFHandler< dim > &background_dof_handler, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells) const
Definition: grid_tools.cc:5109
void process_cell(const typename DoFHandler< dim >::active_cell_iterator &cell, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells) const
Definition: grid_tools.cc:5145
void process_sub_cell(const std::vector< value_type > &, const std::vector< Point< 1 >> &, const std::vector< unsigned int > &, const double, std::vector< Point< 1 >> &, std::vector< CellData< 1 >> &, const bool) const
Definition: grid_tools.h:3014
MarchingCubeAlgorithm(const Mapping< dim, dim > &mapping, const FiniteElement< dim, dim > &fe, const unsigned int n_subdivisions=1, const double tolerance=1e-10)
Definition: grid_tools.cc:5057
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition: grid_tools.cc:5074
const Tensor< 2, 2, double > rotation_matrix
Definition: grid_tools.cc:145
Rotate2d(const double angle)
Definition: grid_tools.cc:134
Point< 2 > operator()(const Point< 2 > &p) const
Definition: grid_tools.cc:139
Rotate3d(const Tensor< 1, 3, double > &axis, const double angle)
Definition: grid_tools.cc:153
Point< 3 > operator()(const Point< 3 > &p) const
Definition: grid_tools.cc:160
const Tensor< 2, 3, double > rotation_matrix
Definition: grid_tools.cc:166
Point< spacedim > operator()(const Point< spacedim > p) const
Definition: grid_tools.cc:178
Scale(const double factor)
Definition: grid_tools.cc:174
Point< spacedim > operator()(const Point< spacedim > p) const
Definition: grid_tools.cc:120
Shift(const Tensor< 1, spacedim > &shift)
Definition: grid_tools.cc:116
const Tensor< 1, spacedim > shift
Definition: grid_tools.cc:126
Abstract base class for mapping classes.
Definition: mapping.h:317
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr numbers::NumberTraits< Number >::real_type square() const
constexpr numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
Quadrature< spacedim > compute_affine_transformation(const std::array< Point< spacedim >, dim+1 > &vertices) const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
size_type n() const
size_type n_rows() const
size_type n_cols() const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
numbers::NumberTraits< Number >::real_type norm() const
virtual std::vector< types::boundary_id > get_boundary_ids() const
const std::vector< bool > & get_used_vertices() const
virtual MPI_Comm get_communicator() const
unsigned int n_quads() const
void load_user_indices(const std::vector< unsigned int > &v)
virtual void clear()
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
bool all_reference_cells_are_hyper_cube() const
void clear_user_data()
face_iterator end_face() const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_lines() const
unsigned int n_raw_faces() const
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
cell_iterator end() const
virtual bool has_hanging_nodes() const
bool vertex_used(const unsigned int index) const
virtual void execute_coarsening_and_refinement()
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, std::bitset< 3 > > > & get_periodic_face_map() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int n_vertices() const
void save_user_indices(std::vector< unsigned int > &v) const
active_face_iterator begin_active_face() const
const std::vector< Point< spacedim > > & get_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
Point< 2 > second
Definition: grid_out.cc:4615
Point< 2 > first
Definition: grid_out.cc:4614
unsigned int vertex_indices[2]
const unsigned int v0
const unsigned int v1
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< active_cell_iterator > active_cell_iterators() const
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsMPI()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1947
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcNeedsCGAL()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:441
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
PackagedOperation< Range > constrained_right_hand_side(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop, const Range &right_hand_side)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
LinearOperator< Range, Domain, Payload > constrained_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop)
@ update_values
Shape function values.
@ update_quadrature_points
Transformed quadrature points.
virtual std::vector< types::manifold_id > get_manifold_ids() const
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:2806
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:2898
void map_boundary_to_manifold_ids(const std::vector< types::boundary_id > &src_boundary_ids, const std::vector< types::manifold_id > &dst_manifold_ids, Triangulation< dim, spacedim > &tria, const std::vector< types::boundary_id > &reset_boundary_ids={})
Definition: grid_tools.cc:2831
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
void assign_co_dimensional_manifold_indicators(Triangulation< dim, spacedim > &tria, const std::function< types::manifold_id(const std::set< types::manifold_id > &)> &disambiguation_function=[](const std::set< types::manifold_id > &manifold_ids) { if(manifold_ids.size()==1) return *manifold_ids.begin();else return numbers::flat_manifold_id;}, bool overwrite_only_flat_manifold_ids=true)
Definition: grid_tools.cc:2926
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Task< RT > new_task(const std::function< RT()> &function)
Expression fabs(const Expression &x)
Expression floor(const Expression &x)
Expression acos(const Expression &x)
void get_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::subdomain_id > &subdomain)
Definition: dof_tools.cc:1596
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::tuple< BoundingBox< MeshType::space_dimension >, bool > compute_cell_predicate_bounding_box(const typename MeshType::cell_iterator &parent_cell, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
Definition: grid_tools.cc:1112
bool fix_up_object(const Iterator &object)
Definition: grid_tools.cc:2538
double objective_function(const Iterator &object, const Point< spacedim > &object_mid_point)
Definition: grid_tools.cc:2350
void fix_up_faces(const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, std::integral_constant< int, dim >, std::integral_constant< int, spacedim >)
Definition: grid_tools.cc:2720
double minimal_diameter(const Iterator &object)
Definition: grid_tools.cc:2511
Point< Iterator::AccessorType::space_dimension > get_face_midpoint(const Iterator &object, const unsigned int f, std::integral_constant< int, 1 >)
Definition: grid_tools.cc:2446
std::vector< std::pair< typename Triangulation< dim, spacedim >::active_cell_iterator, Point< dim > > > find_all_locally_owned_active_cells_around_point(const Cache< dim, spacedim > &cache, const Point< spacedim > &point, typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint, const std::vector< bool > &marked_vertices, const double tolerance, const bool enforce_unique_mapping)
Definition: grid_tools.cc:3742
void laplace_solve(const SparseMatrix< double > &S, const AffineConstraints< double > &constraints, Vector< double > &u)
Definition: grid_tools.cc:273
DistributedComputePointLocationsInternal< dim, spacedim > distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< bool > &marked_vertices, const double tolerance, const bool perform_handshake, const bool enforce_unique_mapping=false)
Definition: grid_tools.cc:3942
bool compare_point_association(const unsigned int a, const unsigned int b, const Tensor< 1, spacedim > &point_direction, const std::vector< Tensor< 1, spacedim >> &center_directions)
Definition: grid_tools.cc:800
std::tuple< std::vector< unsigned int >, std::vector< unsigned int >, std::vector< unsigned int > > guess_owners_of_entities(const MPI_Comm comm, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< T > &entities, const double tolerance)
Definition: grid_tools.cc:3603
DistributedComputeIntersectionLocationsInternal< structdim, spacedim > distributed_compute_intersection_locations(const Cache< dim, spacedim > &cache, const std::vector< std::vector< Point< spacedim >>> &intersection_requests, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< bool > &marked_vertices, const double tolerance)
Definition: grid_tools.cc:4323
void process_sub_cell(const std::array< unsigned int, n_configurations > &cut_line_table, const ndarray< unsigned int, n_configurations, n_cols > &new_line_table, const ndarray< unsigned int, n_lines, 2 > &line_to_vertex_table, const std::vector< value_type > &ls_values, const std::vector< Point< dim >> &points, const std::vector< unsigned int > &mask, const double iso_level, const double tolerance, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells, const bool write_back_cell_data)
Definition: grid_tools.cc:4966
void set_subdomain_id_in_zorder_recursively(IT cell, unsigned int &current_proc_idx, unsigned int &current_cell_idx, const unsigned int n_active_cells, const unsigned int n_partitions)
Definition: grid_tools.cc:1978
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
std::vector< BoundingBox< MeshType::space_dimension > > compute_mesh_predicate_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const unsigned int refinement_level=0, const bool allow_merge=false, const unsigned int max_boxes=numbers::invalid_unsigned_int)
Definition: grid_tools.cc:1168
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:257
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:2007
void collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int >> &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
Definition: grid_tools.cc:4776
Triangulation< dim, spacedim >::DistortedCellList fix_up_distorted_child_cells(const typename Triangulation< dim, spacedim >::DistortedCellList &distorted_cells, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2764
void rotate(const double angle, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:201
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim >> &vertices, const Point< spacedim > &p)
Definition: grid_tools.cc:4571
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
return_type distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &local_points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const double tolerance=1e-10, const std::vector< bool > &marked_vertices={}, const bool enforce_unique_mapping=true)
Definition: grid_tools.cc:3529
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:1080
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={}, const double tolerance=1.e-10)
void delete_duplicated_vertices(std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata, std::vector< unsigned int > &considered_vertices, const double tol=1e-12)
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2298
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:3014
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:191
return_type compute_point_locations_try_all(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:3359
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, const MPI_Comm mpi_communicator)
Definition: grid_tools.cc:4616
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2113
std::vector< std::vector< Tensor< 1, spacedim > > > vertex_to_cell_centers_directions(const Triangulation< dim, spacedim > &mesh, const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator >> &vertex_to_cells)
Definition: grid_tools.cc:761
void laplace_transform(const std::map< unsigned int, Point< dim >> &new_points, Triangulation< dim > &tria, const Function< dim, double > *coefficient=nullptr, const bool solve_for_absolute_positions=false)
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:4899
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:1315
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:2282
return_type compute_point_locations(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:3330
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:1768
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
std::vector< types::subdomain_id > get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const std::vector< CellId > &cell_ids)
Definition: grid_tools.cc:2212
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
Definition: grid_tools.cc:415
double diameter(const Triangulation< dim, spacedim > &tria)
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
std::vector< std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > > find_all_active_cells_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const double tolerance, const std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim >> &first_cell, const std::vector< std::set< typename MeshType< dim, spacedim >::active_cell_iterator >> *vertex_to_cells=nullptr)
std::map< unsigned int, types::global_vertex_index > compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:1413
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, const MPI_Comm mpi_communicator)
Definition: grid_tools.cc:4713
@ valid
Iterator points to a valid object.
static const char U
void create_laplace_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, SparseMatrixType &matrix, const Function< spacedim, typename SparseMatrixType::value_type > *const a=nullptr, const AffineConstraints< typename SparseMatrixType::value_type > &constraints=AffineConstraints< typename SparseMatrixType::value_type >())
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:192
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:490
std::string to_string(const T &t)
Definition: patterns.h:2391
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 2, 2, Number > rotation_matrix_2d(const Number &angle)
Tensor< 2, 3, Number > rotation_matrix_3d(const Tensor< 1, 3, Number > &axis, const Number &angle)