Reference documentation for deal.II version GIT 6d02bd1105 2022-05-22 03:35:01+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 - 2022 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 
16 #include <deal.II/base/mpi.h>
17 #include <deal.II/base/mpi.templates.h>
19 #include <deal.II/base/mpi_consensus_algorithms.templates.h>
22 
27 
30 #include <deal.II/dofs/dof_tools.h>
31 
32 #include <deal.II/fe/fe_nothing.h>
33 #include <deal.II/fe/fe_q.h>
34 #include <deal.II/fe/fe_values.h>
35 #include <deal.II/fe/mapping_q.h>
36 #include <deal.II/fe/mapping_q1.h>
37 
42 #include <deal.II/grid/manifold.h>
43 #include <deal.II/grid/tria.h>
46 
50 #include <deal.II/lac/solver_cg.h>
54 #include <deal.II/lac/vector.h>
56 
59 
61 
62 
64 #include <boost/random/mersenne_twister.hpp>
65 #include <boost/random/uniform_real_distribution.hpp>
67 
68 #include <array>
69 #include <cmath>
70 #include <iostream>
71 #include <list>
72 #include <numeric>
73 #include <set>
74 #include <tuple>
75 #include <unordered_map>
76 
78 
79 
80 namespace GridTools
81 {
82  template <int dim, int spacedim>
83  double
85  {
86  // we can't deal with distributed meshes since we don't have all
87  // vertices locally. there is one exception, however: if the mesh has
88  // never been refined. the way to test this is not to ask
89  // tria.n_levels()==1, since this is something that can happen on one
90  // processor without being true on all. however, we can ask for the
91  // global number of active cells and use that
92 #if defined(DEAL_II_WITH_P4EST) && defined(DEBUG)
94  dynamic_cast<
96  Assert(p_tria->n_global_active_cells() == tria.n_cells(0),
98 #endif
99 
100  // the algorithm used simply traverses all cells and picks out the
101  // boundary vertices. it may or may not be faster to simply get all
102  // vectors, don't mark boundary vertices, and compute the distances
103  // thereof, but at least as the mesh is refined, it seems better to
104  // first mark boundary nodes, as marking is O(N) in the number of
105  // cells/vertices, while computing the maximal distance is O(N*N)
106  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
107  std::vector<bool> boundary_vertices(vertices.size(), false);
108 
110  tria.begin_active();
112  tria.end();
113  for (; cell != endc; ++cell)
114  for (const unsigned int face : cell->face_indices())
115  if (cell->face(face)->at_boundary())
116  for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
117  boundary_vertices[cell->face(face)->vertex_index(i)] = true;
118 
119  // now traverse the list of boundary vertices and check distances.
120  // since distances are symmetric, we only have to check one half
121  double max_distance_sqr = 0;
122  std::vector<bool>::const_iterator pi = boundary_vertices.begin();
123  const unsigned int N = boundary_vertices.size();
124  for (unsigned int i = 0; i < N; ++i, ++pi)
125  {
126  std::vector<bool>::const_iterator pj = pi + 1;
127  for (unsigned int j = i + 1; j < N; ++j, ++pj)
128  if ((*pi == true) && (*pj == true) &&
129  ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
130  max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
131  }
132 
133  return std::sqrt(max_distance_sqr);
134  }
135 
136 
137 
138  template <int dim, int spacedim>
139  double
141  const Mapping<dim, spacedim> & mapping)
142  {
143  // get the degree of the mapping if possible. if not, just assume 1
144  unsigned int mapping_degree = 1;
145  if (const auto *p = dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
146  mapping_degree = p->get_degree();
147  else if (const auto *p =
148  dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
149  mapping_degree = p->get_degree();
150 
151  // then initialize an appropriate quadrature formula
152  const QGauss<dim> quadrature_formula(mapping_degree + 1);
153  const unsigned int n_q_points = quadrature_formula.size();
154 
155  // we really want the JxW values from the FEValues object, but it
156  // wants a finite element. create a cheap element as a dummy
157  // element
158  FE_Nothing<dim, spacedim> dummy_fe;
159  FEValues<dim, spacedim> fe_values(mapping,
160  dummy_fe,
161  quadrature_formula,
163 
165  cell = triangulation.begin_active(),
166  endc = triangulation.end();
167 
168  double local_volume = 0;
169 
170  // compute the integral quantities by quadrature
171  for (; cell != endc; ++cell)
172  if (cell->is_locally_owned())
173  {
174  fe_values.reinit(cell);
175  for (unsigned int q = 0; q < n_q_points; ++q)
176  local_volume += fe_values.JxW(q);
177  }
178 
179  double global_volume = 0;
180 
181 #ifdef DEAL_II_WITH_MPI
183  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
184  &triangulation))
185  global_volume =
186  Utilities::MPI::sum(local_volume, p_tria->get_communicator());
187  else
188 #endif
189  global_volume = local_volume;
190 
191  return global_volume;
192  }
193 
194 
195 
196  namespace
197  {
212  template <int dim>
213  struct TransformR2UAffine
214  {
215  static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
217  };
218 
219 
220  /*
221  Octave code:
222  M=[0 1; 1 1];
223  K1 = transpose(M) * inverse (M*transpose(M));
224  printf ("{%f, %f},\n", K1' );
225  */
226  template <>
228  [1] = {{-1.000000}, {1.000000}};
229 
230  template <>
232  {1.000000, 0.000000};
233 
234 
235  /*
236  Octave code:
237  M=[0 1 0 1;0 0 1 1;1 1 1 1];
238  K2 = transpose(M) * inverse (M*transpose(M));
239  printf ("{%f, %f, %f},\n", K2' );
240  */
241  template <>
243  [2] = {{-0.500000, -0.500000},
244  {0.500000, -0.500000},
245  {-0.500000, 0.500000},
246  {0.500000, 0.500000}};
247 
248  /*
249  Octave code:
250  M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
251  K3 = transpose(M) * inverse (M*transpose(M))
252  printf ("{%f, %f, %f, %f},\n", K3' );
253  */
254  template <>
256  {0.750000, 0.250000, 0.250000, -0.250000};
257 
258 
259  template <>
261  [3] = {
262  {-0.250000, -0.250000, -0.250000},
263  {0.250000, -0.250000, -0.250000},
264  {-0.250000, 0.250000, -0.250000},
265  {0.250000, 0.250000, -0.250000},
266  {-0.250000, -0.250000, 0.250000},
267  {0.250000, -0.250000, 0.250000},
268  {-0.250000, 0.250000, 0.250000},
269  {0.250000, 0.250000, 0.250000}
270 
271  };
272 
273 
274  template <>
276  {0.500000,
277  0.250000,
278  0.250000,
279  0.000000,
280  0.250000,
281  0.000000,
282  0.000000,
283  -0.250000};
284  } // namespace
285 
286 
287 
288  template <int dim, int spacedim>
289  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
291  {
293 
294  // A = vertex * KA
296 
297  for (unsigned int d = 0; d < spacedim; ++d)
298  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
299  for (unsigned int e = 0; e < dim; ++e)
300  A[d][e] += vertices[v][d] * TransformR2UAffine<dim>::KA[v][e];
301 
302  // b = vertex * Kb
304  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
306 
307  return std::make_pair(A, b);
308  }
309 
310 
311 
312  template <int dim>
316  const Quadrature<dim> & quadrature)
317  {
318  FE_Nothing<dim> fe;
319  FEValues<dim> fe_values(mapping, fe, quadrature, update_jacobians);
320 
321  Vector<double> aspect_ratio_vector(triangulation.n_active_cells());
322 
323  // loop over cells of processor
324  for (const auto &cell : triangulation.active_cell_iterators())
325  {
326  if (cell->is_locally_owned())
327  {
328  double aspect_ratio_cell = 0.0;
329 
330  fe_values.reinit(cell);
331 
332  // loop over quadrature points
333  for (unsigned int q = 0; q < quadrature.size(); ++q)
334  {
335  const Tensor<2, dim, double> jacobian =
336  Tensor<2, dim, double>(fe_values.jacobian(q));
337 
338  // We intentionally do not want to throw an exception in case of
339  // inverted elements since this is not the task of this
340  // function. Instead, inf is written into the vector in case of
341  // inverted elements.
342  if (determinant(jacobian) <= 0)
343  {
344  aspect_ratio_cell = std::numeric_limits<double>::infinity();
345  }
346  else
347  {
349  for (unsigned int i = 0; i < dim; ++i)
350  for (unsigned int j = 0; j < dim; ++j)
351  J(i, j) = jacobian[i][j];
352 
353  J.compute_svd();
354 
355  double const max_sv = J.singular_value(0);
356  double const min_sv = J.singular_value(dim - 1);
357  double const ar = max_sv / min_sv;
358 
359  // Take the max between the previous and the current
360  // aspect ratio value; if we had previously encountered
361  // an inverted cell, we will have placed an infinity
362  // in the aspect_ratio_cell variable, and that value
363  // will survive this max operation.
364  aspect_ratio_cell = std::max(aspect_ratio_cell, ar);
365  }
366  }
367 
368  // fill vector
369  aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
370  }
371  }
372 
373  return aspect_ratio_vector;
374  }
375 
376 
377 
378  template <int dim>
379  double
382  const Quadrature<dim> & quadrature)
383  {
384  Vector<double> aspect_ratio_vector =
385  compute_aspect_ratio_of_cells(mapping, triangulation, quadrature);
386 
388  aspect_ratio_vector,
390  }
391 
392 
393 
394  template <int dim, int spacedim>
397  {
398  using iterator =
400  const auto predicate = [](const iterator &) { return true; };
401 
402  return compute_bounding_box(
403  tria, std::function<bool(const iterator &)>(predicate));
404  }
405 
406 
407 
408  // Generic functions for appending face data in 2D or 3D. TODO: we can
409  // remove these once we have 'if constexpr'.
410  namespace internal
411  {
412  inline void
413  append_face_data(const CellData<1> &face_data, SubCellData &subcell_data)
414  {
415  subcell_data.boundary_lines.push_back(face_data);
416  }
417 
418 
419 
420  inline void
421  append_face_data(const CellData<2> &face_data, SubCellData &subcell_data)
422  {
423  subcell_data.boundary_quads.push_back(face_data);
424  }
425 
426 
427 
428  // Lexical comparison for sorting CellData objects.
429  template <int structdim>
431  {
432  bool
434  const CellData<structdim> &b) const
435  {
436  // Check vertices:
437  if (std::lexicographical_compare(std::begin(a.vertices),
438  std::end(a.vertices),
439  std::begin(b.vertices),
440  std::end(b.vertices)))
441  return true;
442  // it should never be necessary to check the material or manifold
443  // ids as a 'tiebreaker' (since they must be equal if the vertex
444  // indices are equal). Assert it anyway:
445 #ifdef DEBUG
446  if (std::equal(std::begin(a.vertices),
447  std::end(a.vertices),
448  std::begin(b.vertices)))
449  {
450  Assert(a.material_id == b.material_id &&
451  a.manifold_id == b.manifold_id,
452  ExcMessage(
453  "Two CellData objects with equal vertices must "
454  "have the same material/boundary ids and manifold "
455  "ids."));
456  }
457 #endif
458  return false;
459  }
460  };
461 
462 
472  template <int dim>
474  {
475  public:
479  template <class FaceIteratorType>
480  void
481  insert_face_data(const FaceIteratorType &face)
482  {
483  CellData<dim - 1> face_cell_data(face->n_vertices());
484  for (unsigned int vertex_n = 0; vertex_n < face->n_vertices();
485  ++vertex_n)
486  face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
487  face_cell_data.boundary_id = face->boundary_id();
488  face_cell_data.manifold_id = face->manifold_id();
489 
490  face_data.insert(std::move(face_cell_data));
491  }
492 
497  get()
498  {
499  SubCellData subcell_data;
500 
501  for (const CellData<dim - 1> &face_cell_data : face_data)
502  internal::append_face_data(face_cell_data, subcell_data);
503  return subcell_data;
504  }
505 
506 
507  private:
510  };
511 
512 
513  // Do nothing for dim=1:
514  template <>
515  class FaceDataHelper<1>
516  {
517  public:
518  template <class FaceIteratorType>
519  void
520  insert_face_data(const FaceIteratorType &)
521  {}
522 
524  get()
525  {
526  return SubCellData();
527  }
528  };
529  } // namespace internal
530 
531 
532 
533  template <int dim, int spacedim>
534  std::
535  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
537  {
538  Assert(1 <= tria.n_levels(),
539  ExcMessage("The input triangulation must be non-empty."));
540 
541  std::vector<Point<spacedim>> vertices;
542  std::vector<CellData<dim>> cells;
543 
544  unsigned int max_level_0_vertex_n = 0;
545  for (const auto &cell : tria.cell_iterators_on_level(0))
546  for (const unsigned int cell_vertex_n : cell->vertex_indices())
547  max_level_0_vertex_n =
548  std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
549  vertices.resize(max_level_0_vertex_n + 1);
550 
552  std::set<CellData<1>, internal::CellDataComparator<1>>
553  line_data; // only used in 3D
554 
555  for (const auto &cell : tria.cell_iterators_on_level(0))
556  {
557  // Save cell data
558  CellData<dim> cell_data(cell->n_vertices());
559  for (const unsigned int cell_vertex_n : cell->vertex_indices())
560  {
561  Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
562  ExcInternalError());
563  vertices[cell->vertex_index(cell_vertex_n)] =
564  cell->vertex(cell_vertex_n);
565  cell_data.vertices[cell_vertex_n] =
566  cell->vertex_index(cell_vertex_n);
567  }
568  cell_data.material_id = cell->material_id();
569  cell_data.manifold_id = cell->manifold_id();
570  cells.push_back(cell_data);
571 
572  // Save face data
573  if (dim > 1)
574  {
575  for (const unsigned int face_n : cell->face_indices())
576  // We don't need to insert anything if we have default values
577  {
578  const auto face = cell->face(face_n);
579  if (face->boundary_id() != numbers::internal_face_boundary_id ||
580  face->manifold_id() != numbers::flat_manifold_id)
581  face_data.insert_face_data(face);
582  }
583  }
584  // Save line data
585  if (dim == 3)
586  {
587  for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
588  {
589  const auto line = cell->line(line_n);
590  // We don't need to insert anything if we have default values
591  if (line->boundary_id() != numbers::internal_face_boundary_id ||
592  line->manifold_id() != numbers::flat_manifold_id)
593  {
594  CellData<1> line_cell_data(line->n_vertices());
595  for (unsigned int vertex_n : line->vertex_indices())
596  line_cell_data.vertices[vertex_n] =
597  line->vertex_index(vertex_n);
598  line_cell_data.boundary_id = line->boundary_id();
599  line_cell_data.manifold_id = line->manifold_id();
600  line_data.insert(std::move(line_cell_data));
601  }
602  }
603  }
604  }
605 
606  // Double-check that there are no unused vertices:
607 #ifdef DEBUG
608  {
609  std::vector<bool> used_vertices(vertices.size());
610  for (const CellData<dim> &cell_data : cells)
611  for (const auto v : cell_data.vertices)
612  used_vertices[v] = true;
613  Assert(std::find(used_vertices.begin(), used_vertices.end(), false) ==
614  used_vertices.end(),
615  ExcMessage("The level zero vertices should form a contiguous "
616  "range."));
617  }
618 #endif
619 
620  SubCellData subcell_data = face_data.get();
621 
622  if (dim == 3)
623  for (const CellData<1> &face_line_data : line_data)
624  subcell_data.boundary_lines.push_back(face_line_data);
625 
626  return std::tuple<std::vector<Point<spacedim>>,
627  std::vector<CellData<dim>>,
628  SubCellData>(std::move(vertices),
629  std::move(cells),
630  std::move(subcell_data));
631  }
632 
633 
634 
635  template <int dim, int spacedim>
636  void
638  std::vector<CellData<dim>> & cells,
639  SubCellData & subcelldata)
640  {
641  Assert(
642  subcelldata.check_consistency(dim),
643  ExcMessage(
644  "Invalid SubCellData supplied according to ::check_consistency(). "
645  "This is caused by data containing objects for the wrong dimension."));
646 
647  // first check which vertices are actually used
648  std::vector<bool> vertex_used(vertices.size(), false);
649  for (unsigned int c = 0; c < cells.size(); ++c)
650  for (unsigned int v = 0; v < cells[c].vertices.size(); ++v)
651  {
652  Assert(cells[c].vertices[v] < vertices.size(),
653  ExcMessage("Invalid vertex index encountered! cells[" +
654  Utilities::int_to_string(c) + "].vertices[" +
655  Utilities::int_to_string(v) + "]=" +
656  Utilities::int_to_string(cells[c].vertices[v]) +
657  " is invalid, because only " +
659  " vertices were supplied."));
660  vertex_used[cells[c].vertices[v]] = true;
661  }
662 
663 
664  // then renumber the vertices that are actually used in the same order as
665  // they were beforehand
666  const unsigned int invalid_vertex = numbers::invalid_unsigned_int;
667  std::vector<unsigned int> new_vertex_numbers(vertices.size(),
668  invalid_vertex);
669  unsigned int next_free_number = 0;
670  for (unsigned int i = 0; i < vertices.size(); ++i)
671  if (vertex_used[i] == true)
672  {
673  new_vertex_numbers[i] = next_free_number;
674  ++next_free_number;
675  }
676 
677  // next replace old vertex numbers by the new ones
678  for (unsigned int c = 0; c < cells.size(); ++c)
679  for (auto &v : cells[c].vertices)
680  v = new_vertex_numbers[v];
681 
682  // same for boundary data
683  for (unsigned int c = 0; c < subcelldata.boundary_lines.size(); // NOLINT
684  ++c)
685  for (unsigned int v = 0;
686  v < subcelldata.boundary_lines[c].vertices.size();
687  ++v)
688  {
689  Assert(subcelldata.boundary_lines[c].vertices[v] <
690  new_vertex_numbers.size(),
691  ExcMessage(
692  "Invalid vertex index in subcelldata.boundary_lines. "
693  "subcelldata.boundary_lines[" +
694  Utilities::int_to_string(c) + "].vertices[" +
695  Utilities::int_to_string(v) + "]=" +
697  subcelldata.boundary_lines[c].vertices[v]) +
698  " is invalid, because only " +
700  " vertices were supplied."));
701  subcelldata.boundary_lines[c].vertices[v] =
702  new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]];
703  }
704 
705  for (unsigned int c = 0; c < subcelldata.boundary_quads.size(); // NOLINT
706  ++c)
707  for (unsigned int v = 0;
708  v < subcelldata.boundary_quads[c].vertices.size();
709  ++v)
710  {
711  Assert(subcelldata.boundary_quads[c].vertices[v] <
712  new_vertex_numbers.size(),
713  ExcMessage(
714  "Invalid vertex index in subcelldata.boundary_quads. "
715  "subcelldata.boundary_quads[" +
716  Utilities::int_to_string(c) + "].vertices[" +
717  Utilities::int_to_string(v) + "]=" +
719  subcelldata.boundary_quads[c].vertices[v]) +
720  " is invalid, because only " +
722  " vertices were supplied."));
723 
724  subcelldata.boundary_quads[c].vertices[v] =
725  new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]];
726  }
727 
728  // finally copy over the vertices which we really need to a new array and
729  // replace the old one by the new one
730  std::vector<Point<spacedim>> tmp;
731  tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(), true));
732  for (unsigned int v = 0; v < vertices.size(); ++v)
733  if (vertex_used[v] == true)
734  tmp.push_back(vertices[v]);
735  swap(vertices, tmp);
736  }
737 
738 
739 
740  template <int dim, int spacedim>
741  void
743  std::vector<CellData<dim>> & cells,
744  SubCellData & subcelldata,
745  std::vector<unsigned int> & considered_vertices,
746  const double tol)
747  {
748  AssertIndexRange(2, vertices.size());
749  std::vector<unsigned int> new_vertex_numbers(vertices.size());
750  std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
751 
752  // if the considered_vertices vector is empty, consider all vertices
753  if (considered_vertices.size() == 0)
754  considered_vertices = new_vertex_numbers;
755  Assert(considered_vertices.size() <= vertices.size(), ExcInternalError());
756 
757  // The algorithm below improves upon the naive O(n^2) algorithm by first
758  // sorting vertices by their value in one component and then only
759  // comparing vertices for equality which are nearly equal in that
760  // component. For example, if @p vertices form a cube, then we will only
761  // compare points that have the same x coordinate when we try to find
762  // duplicated vertices.
763 
764  // Start by finding the longest coordinate direction. This minimizes the
765  // number of points that need to be compared against each-other in a
766  // single set for typical geometries.
767  const BoundingBox<spacedim> bbox(vertices);
768 
769  unsigned int longest_coordinate_direction = 0;
770  double longest_coordinate_length = bbox.side_length(0);
771  for (unsigned int d = 1; d < spacedim; ++d)
772  {
773  const double coordinate_length = bbox.side_length(d);
774  if (longest_coordinate_length < coordinate_length)
775  {
776  longest_coordinate_length = coordinate_length;
777  longest_coordinate_direction = d;
778  }
779  }
780 
781  // Sort vertices (while preserving their vertex numbers) along that
782  // coordinate direction:
783  std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
784  sorted_vertices.reserve(vertices.size());
785  for (const unsigned int vertex_n : considered_vertices)
786  {
787  AssertIndexRange(vertex_n, vertices.size());
788  sorted_vertices.emplace_back(vertex_n, vertices[vertex_n]);
789  }
790  std::sort(sorted_vertices.begin(),
791  sorted_vertices.end(),
792  [&](const std::pair<unsigned int, Point<spacedim>> &a,
793  const std::pair<unsigned int, Point<spacedim>> &b) {
794  return a.second[longest_coordinate_direction] <
795  b.second[longest_coordinate_direction];
796  });
797 
798  auto within_tolerance = [=](const Point<spacedim> &a,
799  const Point<spacedim> &b) {
800  for (unsigned int d = 0; d < spacedim; ++d)
801  if (std::abs(a[d] - b[d]) > tol)
802  return false;
803  return true;
804  };
805 
806  // Find a range of numbers that have the same component in the longest
807  // coordinate direction:
808  auto range_start = sorted_vertices.begin();
809  while (range_start != sorted_vertices.end())
810  {
811  auto range_end = range_start + 1;
812  while (range_end != sorted_vertices.end() &&
813  std::abs(range_end->second[longest_coordinate_direction] -
814  range_start->second[longest_coordinate_direction]) <
815  tol)
816  ++range_end;
817 
818  // preserve behavior with older versions of this function by replacing
819  // higher vertex numbers by lower vertex numbers
820  std::sort(range_start,
821  range_end,
822  [](const std::pair<unsigned int, Point<spacedim>> &a,
823  const std::pair<unsigned int, Point<spacedim>> &b) {
824  return a.first < b.first;
825  });
826 
827  // Now de-duplicate [range_start, range_end)
828  //
829  // We have identified all points that are within a strip of width 'tol'
830  // in one coordinate direction. Now we need to figure out which of these
831  // are also close in other coordinate directions. If two are close, we
832  // can mark the second one for deletion.
833  for (auto reference = range_start; reference != range_end; ++reference)
834  {
835  if (reference->first != numbers::invalid_unsigned_int)
836  for (auto it = reference + 1; it != range_end; ++it)
837  {
838  if (within_tolerance(reference->second, it->second))
839  {
840  new_vertex_numbers[it->first] = reference->first;
841  // skip the replaced vertex in the future
842  it->first = numbers::invalid_unsigned_int;
843  }
844  }
845  }
846  range_start = range_end;
847  }
848 
849  // now we got a renumbering list. simply renumber all vertices
850  // (non-duplicate vertices get renumbered to themselves, so nothing bad
851  // happens). after that, the duplicate vertices will be unused, so call
852  // delete_unused_vertices() to do that part of the job.
853  for (auto &cell : cells)
854  for (auto &vertex_index : cell.vertices)
855  vertex_index = new_vertex_numbers[vertex_index];
856  for (auto &quad : subcelldata.boundary_quads)
857  for (auto &vertex_index : quad.vertices)
858  vertex_index = new_vertex_numbers[vertex_index];
859  for (auto &line : subcelldata.boundary_lines)
860  for (auto &vertex_index : line.vertices)
861  vertex_index = new_vertex_numbers[vertex_index];
862 
863  delete_unused_vertices(vertices, cells, subcelldata);
864  }
865 
866 
867 
868  template <int dim, int spacedim>
869  std::size_t
871  const std::vector<Point<spacedim>> &all_vertices,
872  std::vector<CellData<dim>> & cells)
873  {
874  // This function is presently only implemented for hypercube and simplex
875  // volumetric (codimension 0) elements.
876 
877  if (dim == 1)
878  return 0;
879  if (dim == 2 && spacedim == 3)
880  Assert(false, ExcNotImplemented());
881 
882  std::size_t n_negative_cells = 0;
883  for (auto &cell : cells)
884  {
885  const ArrayView<const unsigned int> vertices(cell.vertices);
886  if (GridTools::cell_measure(all_vertices, vertices) < 0)
887  {
888  const unsigned int n_vertices = vertices.size();
889 
890  if (ReferenceCell::n_vertices_to_type(dim, n_vertices)
891  .is_hyper_cube())
892  {
893  ++n_negative_cells;
894 
895  if (dim == 2)
896  {
897  // flip the cell across the y = x line in 2D
898  std::swap(cell.vertices[1], cell.vertices[2]);
899  }
900  else if (dim == 3)
901  {
902  // swap the front and back faces in 3D
903  std::swap(cell.vertices[0], cell.vertices[2]);
904  std::swap(cell.vertices[1], cell.vertices[3]);
905  std::swap(cell.vertices[4], cell.vertices[6]);
906  std::swap(cell.vertices[5], cell.vertices[7]);
907  }
908  }
909 
910  else if (ReferenceCell::n_vertices_to_type(dim, n_vertices)
911  .is_simplex())
912  {
913  ++n_negative_cells;
914  // By basic rules for computing determinants we can just swap
915  // two vertices to fix a negative volume. Arbitrarily pick the
916  // last two.
917  std::swap(cell.vertices[n_vertices - 2],
918  cell.vertices[n_vertices - 1]);
919  }
920  else
921  {
922  AssertThrow(false, ExcNotImplemented());
923  }
924  // Check whether the resulting cell is now ok.
925  // If not, then the grid is seriously broken and
926  // we just give up.
927  AssertThrow(GridTools::cell_measure(all_vertices, vertices) > 0,
928  ExcInternalError());
929  }
930  }
931  return n_negative_cells;
932  }
933 
934 
935  template <int dim, int spacedim>
936  void
938  const std::vector<Point<spacedim>> &all_vertices,
939  std::vector<CellData<dim>> & cells)
940  {
941  const std::size_t n_negative_cells =
942  invert_cells_with_negative_measure(all_vertices, cells);
943 
944  // We assume that all cells of a grid have
945  // either positive or negative volumes but
946  // not both mixed. Although above reordering
947  // might work also on single cells, grids
948  // with both kind of cells are very likely to
949  // be broken. Check for this here.
950  AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
951  ExcMessage(
952  std::string(
953  "This function assumes that either all cells have positive "
954  "volume, or that all cells have been specified in an "
955  "inverted vertex order so that their volume is negative. "
956  "(In the latter case, this class automatically inverts "
957  "every cell.) However, the mesh you have specified "
958  "appears to have both cells with positive and cells with "
959  "negative volume. You need to check your mesh which "
960  "cells these are and how they got there.\n"
961  "As a hint, of the total ") +
962  std::to_string(cells.size()) + " cells in the mesh, " +
963  std::to_string(n_negative_cells) +
964  " appear to have a negative volume."));
965  }
966 
967 
968 
969  // Functions and classes for consistently_order_cells
970  namespace
971  {
977  struct CheapEdge
978  {
982  CheapEdge(const unsigned int v0, const unsigned int v1)
983  : v0(v0)
984  , v1(v1)
985  {}
986 
991  bool
992  operator<(const CheapEdge &e) const
993  {
994  return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1)));
995  }
996 
997  private:
1001  const unsigned int v0, v1;
1002  };
1003 
1004 
1013  template <int dim>
1014  bool
1015  is_consistent(const std::vector<CellData<dim>> &cells)
1016  {
1017  std::set<CheapEdge> edges;
1018 
1019  for (typename std::vector<CellData<dim>>::const_iterator c =
1020  cells.begin();
1021  c != cells.end();
1022  ++c)
1023  {
1024  // construct the edges in reverse order. for each of them,
1025  // ensure that the reverse edge is not yet in the list of
1026  // edges (return false if the reverse edge already *is* in
1027  // the list) and then add the actual edge to it; std::set
1028  // eliminates duplicates automatically
1029  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1030  {
1031  const CheapEdge reverse_edge(
1033  c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 0)]);
1034  if (edges.find(reverse_edge) != edges.end())
1035  return false;
1036 
1037 
1038  // ok, not. insert edge in correct order
1039  const CheapEdge correct_edge(
1041  c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)]);
1042  edges.insert(correct_edge);
1043  }
1044  }
1045 
1046  // no conflicts found, so return true
1047  return true;
1048  }
1049 
1050 
1057  template <int dim>
1058  struct ParallelEdges
1059  {
1065  static const unsigned int starter_edges[dim];
1066 
1071  static const unsigned int n_other_parallel_edges = (1 << (dim - 1)) - 1;
1072  static const unsigned int
1075  };
1076 
1077  template <>
1078  const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2};
1079 
1080  template <>
1081  const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1},
1082  {0},
1083  {3},
1084  {2}};
1085 
1086  template <>
1087  const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8};
1088 
1089  template <>
1090  const unsigned int ParallelEdges<3>::parallel_edges[12][3] = {
1091  {1, 4, 5}, // line 0
1092  {0, 4, 5}, // line 1
1093  {3, 6, 7}, // line 2
1094  {2, 6, 7}, // line 3
1095  {0, 1, 5}, // line 4
1096  {0, 1, 4}, // line 5
1097  {2, 3, 7}, // line 6
1098  {2, 3, 6}, // line 7
1099  {9, 10, 11}, // line 8
1100  {8, 10, 11}, // line 9
1101  {8, 9, 11}, // line 10
1102  {8, 9, 10} // line 11
1103  };
1104 
1105 
1110  struct AdjacentCell
1111  {
1115  AdjacentCell()
1118  {}
1119 
1123  AdjacentCell(const unsigned int cell_index,
1124  const unsigned int edge_within_cell)
1127  {}
1128 
1129 
1130  unsigned int cell_index;
1131  unsigned int edge_within_cell;
1132  };
1133 
1134 
1135 
1136  template <int dim>
1137  class AdjacentCells;
1138 
1144  template <>
1145  class AdjacentCells<2>
1146  {
1147  public:
1152  using const_iterator = const AdjacentCell *;
1153 
1162  void
1163  push_back(const AdjacentCell &adjacent_cell)
1164  {
1166  adjacent_cells[0] = adjacent_cell;
1167  else
1168  {
1171  ExcInternalError());
1172  adjacent_cells[1] = adjacent_cell;
1173  }
1174  }
1175 
1176 
1181  const_iterator
1182  begin() const
1183  {
1184  return adjacent_cells;
1185  }
1186 
1187 
1193  const_iterator
1194  end() const
1195  {
1196  // check whether the current object stores zero, one, or two
1197  // adjacent cells, and use this to point to the element past the
1198  // last valid one
1200  return adjacent_cells;
1202  return adjacent_cells + 1;
1203  else
1204  return adjacent_cells + 2;
1205  }
1206 
1207  private:
1214  AdjacentCell adjacent_cells[2];
1215  };
1216 
1217 
1218 
1226  template <>
1227  class AdjacentCells<3> : public std::vector<AdjacentCell>
1228  {};
1229 
1230 
1240  template <int dim>
1241  class Edge
1242  {
1243  public:
1249  Edge(const CellData<dim> &cell, const unsigned int edge_number)
1250  : orientation_status(not_oriented)
1251  {
1253  ExcInternalError());
1254 
1255  // copy vertices for this particular line
1256  vertex_indices[0] =
1257  cell
1259  vertex_indices[1] =
1260  cell
1262 
1263  // bring them into standard orientation
1264  if (vertex_indices[0] > vertex_indices[1])
1266  }
1267 
1272  bool
1273  operator<(const Edge<dim> &e) const
1274  {
1275  return ((vertex_indices[0] < e.vertex_indices[0]) ||
1276  ((vertex_indices[0] == e.vertex_indices[0]) &&
1277  (vertex_indices[1] < e.vertex_indices[1])));
1278  }
1279 
1283  bool
1284  operator==(const Edge<dim> &e) const
1285  {
1286  return ((vertex_indices[0] == e.vertex_indices[0]) &&
1287  (vertex_indices[1] == e.vertex_indices[1]));
1288  }
1289 
1294  unsigned int vertex_indices[2];
1295 
1300  enum OrientationStatus
1301  {
1302  not_oriented,
1303  forward,
1304  backward
1305  };
1306 
1307  OrientationStatus orientation_status;
1308 
1313  AdjacentCells<dim> adjacent_cells;
1314  };
1315 
1316 
1317 
1322  template <int dim>
1323  struct Cell
1324  {
1330  Cell(const CellData<dim> &c, const std::vector<Edge<dim>> &edge_list)
1331  {
1332  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1333  vertex_indices[i] = c.vertices[i];
1334 
1335  // now for each of the edges of this cell, find the location inside the
1336  // given edge_list array and store than index
1337  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1338  {
1339  const Edge<dim> e(c, l);
1340  edge_indices[l] =
1341  (std::lower_bound(edge_list.begin(), edge_list.end(), e) -
1342  edge_list.begin());
1343  Assert(edge_indices[l] < edge_list.size(), ExcInternalError());
1344  Assert(edge_list[edge_indices[l]] == e, ExcInternalError())
1345  }
1346  }
1347 
1352 
1358  };
1359 
1360 
1361 
1362  template <int dim>
1363  class EdgeDeltaSet;
1364 
1374  template <>
1375  class EdgeDeltaSet<2>
1376  {
1377  public:
1381  using const_iterator = const unsigned int *;
1382 
1387  EdgeDeltaSet()
1388  {
1390  }
1391 
1392 
1396  void
1397  clear()
1398  {
1400  }
1401 
1406  void
1407  insert(const unsigned int edge_index)
1408  {
1410  edge_indices[0] = edge_index;
1411  else
1412  {
1414  ExcInternalError());
1415  edge_indices[1] = edge_index;
1416  }
1417  }
1418 
1419 
1423  const_iterator
1424  begin() const
1425  {
1426  return edge_indices;
1427  }
1428 
1429 
1433  const_iterator
1434  end() const
1435  {
1436  // check whether the current object stores zero, one, or two
1437  // indices, and use this to point to the element past the
1438  // last valid one
1440  return edge_indices;
1442  return edge_indices + 1;
1443  else
1444  return edge_indices + 2;
1445  }
1446 
1447  private:
1451  unsigned int edge_indices[2];
1452  };
1453 
1454 
1455 
1467  template <>
1468  class EdgeDeltaSet<3> : public std::set<unsigned int>
1469  {};
1470 
1471 
1472 
1477  template <int dim>
1478  std::vector<Edge<dim>>
1479  build_edges(const std::vector<CellData<dim>> &cells)
1480  {
1481  // build the edge list for all cells. because each cell has
1482  // GeometryInfo<dim>::lines_per_cell edges, the total number
1483  // of edges is this many times the number of cells. of course
1484  // some of them will be duplicates, and we throw them out below
1485  std::vector<Edge<dim>> edge_list;
1486  edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
1487  for (unsigned int i = 0; i < cells.size(); ++i)
1488  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1489  edge_list.emplace_back(cells[i], l);
1490 
1491  // next sort the edge list and then remove duplicates
1492  std::sort(edge_list.begin(), edge_list.end());
1493  edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
1494  edge_list.end());
1495 
1496  return edge_list;
1497  }
1498 
1499 
1500 
1505  template <int dim>
1506  std::vector<Cell<dim>>
1507  build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
1508  std::vector<Edge<dim>> & edges)
1509  {
1510  std::vector<Cell<dim>> cell_list;
1511  cell_list.reserve(cells.size());
1512  for (unsigned int i = 0; i < cells.size(); ++i)
1513  {
1514  // create our own data structure for the cells and let it
1515  // connect to the edges array
1516  cell_list.emplace_back(cells[i], edges);
1517 
1518  // then also inform the edges that they are adjacent
1519  // to the current cell, and where within this cell
1520  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1521  edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back(
1522  AdjacentCell(i, l));
1523  }
1524  Assert(cell_list.size() == cells.size(), ExcInternalError());
1525 
1526  return cell_list;
1527  }
1528 
1529 
1530 
1535  template <int dim>
1536  unsigned int
1537  get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
1538  const std::vector<Edge<dim>> &edges,
1539  const unsigned int current_cell)
1540  {
1541  for (unsigned int c = current_cell; c < cells.size(); ++c)
1542  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1543  if (edges[cells[c].edge_indices[l]].orientation_status ==
1544  Edge<dim>::not_oriented)
1545  return c;
1546 
1548  }
1549 
1550 
1551 
1557  template <int dim>
1558  void
1559  orient_one_set_of_parallel_edges(const std::vector<Cell<dim>> &cells,
1560  std::vector<Edge<dim>> & edges,
1561  const unsigned int cell,
1562  const unsigned int local_edge)
1563  {
1564  // choose the direction of the first edge. we have free choice
1565  // here and could simply choose "forward" if that's what pleases
1566  // us. however, for backward compatibility with the previous
1567  // implementation used till 2016, let us just choose the
1568  // direction so that it matches what we have in the given cell.
1569  //
1570  // in fact, in what can only be assumed to be a bug in the
1571  // original implementation, after orienting all edges, the code
1572  // that rotates the cells so that they match edge orientations
1573  // (see the rotate_cell() function below) rotated the cell two
1574  // more times by 90 degrees. this is ok -- it simply flips all
1575  // edge orientations, which leaves them valid. rather than do
1576  // the same in the current implementation, we can achieve the
1577  // same effect by modifying the rule above to choose the
1578  // direction of the starting edge of this parallel set
1579  // *opposite* to what it looks like in the current cell
1580  //
1581  // this bug only existed in the 2d implementation since there
1582  // were different implementations for 2d and 3d. consequently,
1583  // only replicate it for the 2d case and be "intuitive" in 3d.
1584  if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
1586  local_edge, 0)])
1587  // orient initial edge *opposite* to the way it is in the cell
1588  // (see above for the reason)
1589  edges[cells[cell].edge_indices[local_edge]].orientation_status =
1590  (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
1591  else
1592  {
1593  Assert(
1594  edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
1595  cells[cell].vertex_indices
1597  ExcInternalError());
1598  Assert(
1599  edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
1600  cells[cell].vertex_indices
1602  ExcInternalError());
1603 
1604  // orient initial edge *opposite* to the way it is in the cell
1605  // (see above for the reason)
1606  edges[cells[cell].edge_indices[local_edge]].orientation_status =
1607  (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
1608  }
1609 
1610  // walk outward from the given edge as described in
1611  // the algorithm in the paper that documents all of
1612  // this
1613  //
1614  // note that in 2d, each of the Deltas can at most
1615  // contain two elements, whereas in 3d it can be arbitrarily many
1616  EdgeDeltaSet<dim> Delta_k;
1617  EdgeDeltaSet<dim> Delta_k_minus_1;
1618  Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]);
1619 
1620  while (Delta_k_minus_1.begin() !=
1621  Delta_k_minus_1.end()) // while set is not empty
1622  {
1623  Delta_k.clear();
1624 
1625  for (typename EdgeDeltaSet<dim>::const_iterator delta =
1626  Delta_k_minus_1.begin();
1627  delta != Delta_k_minus_1.end();
1628  ++delta)
1629  {
1630  Assert(edges[*delta].orientation_status !=
1631  Edge<dim>::not_oriented,
1632  ExcInternalError());
1633 
1634  // now go through the cells adjacent to this edge
1635  for (typename AdjacentCells<dim>::const_iterator adjacent_cell =
1636  edges[*delta].adjacent_cells.begin();
1637  adjacent_cell != edges[*delta].adjacent_cells.end();
1638  ++adjacent_cell)
1639  {
1640  const unsigned int K = adjacent_cell->cell_index;
1641  const unsigned int delta_is_edge_in_K =
1642  adjacent_cell->edge_within_cell;
1643 
1644  // figure out the direction of delta with respect to the cell
1645  // K (in the orientation in which the user has given it to us)
1646  const unsigned int first_edge_vertex =
1647  (edges[*delta].orientation_status == Edge<dim>::forward ?
1648  edges[*delta].vertex_indices[0] :
1649  edges[*delta].vertex_indices[1]);
1650  const unsigned int first_edge_vertex_in_K =
1651  cells[K]
1653  delta_is_edge_in_K, 0)];
1654  Assert(
1655  first_edge_vertex == first_edge_vertex_in_K ||
1656  first_edge_vertex ==
1657  cells[K].vertex_indices[GeometryInfo<
1658  dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
1659  ExcInternalError());
1660 
1661  // now figure out which direction the each of the "opposite"
1662  // edges needs to be oriented into.
1663  for (unsigned int o_e = 0;
1665  ++o_e)
1666  {
1667  // get the index of the opposite edge and select which its
1668  // first vertex needs to be based on how the current edge
1669  // is oriented in the current cell
1670  const unsigned int opposite_edge =
1671  cells[K].edge_indices[ParallelEdges<
1672  dim>::parallel_edges[delta_is_edge_in_K][o_e]];
1673  const unsigned int first_opposite_edge_vertex =
1674  cells[K].vertex_indices
1676  ParallelEdges<
1677  dim>::parallel_edges[delta_is_edge_in_K][o_e],
1678  (first_edge_vertex == first_edge_vertex_in_K ? 0 :
1679  1))];
1680 
1681  // then determine the orientation of the edge based on
1682  // whether the vertex we want to be the edge's first
1683  // vertex is already the first vertex of the edge, or
1684  // whether it points in the opposite direction
1685  const typename Edge<dim>::OrientationStatus
1686  opposite_edge_orientation =
1687  (edges[opposite_edge].vertex_indices[0] ==
1688  first_opposite_edge_vertex ?
1689  Edge<dim>::forward :
1690  Edge<dim>::backward);
1691 
1692  // see if the opposite edge (there is only one in 2d) has
1693  // already been oriented.
1694  if (edges[opposite_edge].orientation_status ==
1695  Edge<dim>::not_oriented)
1696  {
1697  // the opposite edge is not yet oriented. do orient it
1698  // and add it to Delta_k
1699  edges[opposite_edge].orientation_status =
1700  opposite_edge_orientation;
1701  Delta_k.insert(opposite_edge);
1702  }
1703  else
1704  {
1705  // this opposite edge has already been oriented. it
1706  // should be consistent with the current one in 2d,
1707  // while in 3d it may in fact be mis-oriented, and in
1708  // that case the mesh will not be orientable. indicate
1709  // this by throwing an exception that we can catch
1710  // further up; this has the advantage that we can
1711  // propagate through a couple of functions without
1712  // having to do error checking and without modifying
1713  // the 'cells' array that the user gave us
1714  if (dim == 2)
1715  {
1716  Assert(edges[opposite_edge].orientation_status ==
1717  opposite_edge_orientation,
1719  }
1720  else if (dim == 3)
1721  {
1722  if (edges[opposite_edge].orientation_status !=
1723  opposite_edge_orientation)
1724  throw ExcMeshNotOrientable();
1725  }
1726  else
1727  Assert(false, ExcNotImplemented());
1728  }
1729  }
1730  }
1731  }
1732 
1733  // finally copy the new set to the previous one
1734  // (corresponding to increasing 'k' by one in the
1735  // algorithm)
1736  Delta_k_minus_1 = Delta_k;
1737  }
1738  }
1739 
1740 
1748  template <int dim>
1749  void
1750  rotate_cell(const std::vector<Cell<dim>> &cell_list,
1751  const std::vector<Edge<dim>> &edge_list,
1752  const unsigned int cell_index,
1753  std::vector<CellData<dim>> & raw_cells)
1754  {
1755  // find the first vertex of the cell. this is the vertex where dim edges
1756  // originate, so for each of the edges record which the starting vertex is
1757  unsigned int starting_vertex_of_edge[GeometryInfo<dim>::lines_per_cell];
1758  for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
1759  {
1760  Assert(edge_list[cell_list[cell_index].edge_indices[e]]
1761  .orientation_status != Edge<dim>::not_oriented,
1762  ExcInternalError());
1763  if (edge_list[cell_list[cell_index].edge_indices[e]]
1764  .orientation_status == Edge<dim>::forward)
1765  starting_vertex_of_edge[e] =
1766  edge_list[cell_list[cell_index].edge_indices[e]]
1767  .vertex_indices[0];
1768  else
1769  starting_vertex_of_edge[e] =
1770  edge_list[cell_list[cell_index].edge_indices[e]]
1771  .vertex_indices[1];
1772  }
1773 
1774  // find the vertex number that appears dim times. this will then be
1775  // the vertex at which we want to locate the origin of the cell's
1776  // coordinate system (i.e., vertex 0)
1777  unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int;
1778  switch (dim)
1779  {
1780  case 2:
1781  {
1782  // in 2d, we can simply enumerate the possibilities where the
1783  // origin may be located because edges zero and one don't share
1784  // any vertices, and the same for edges two and three
1785  if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
1786  (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
1787  origin_vertex_of_cell = starting_vertex_of_edge[0];
1788  else if ((starting_vertex_of_edge[1] ==
1789  starting_vertex_of_edge[2]) ||
1790  (starting_vertex_of_edge[1] ==
1791  starting_vertex_of_edge[3]))
1792  origin_vertex_of_cell = starting_vertex_of_edge[1];
1793  else
1794  Assert(false, ExcInternalError());
1795 
1796  break;
1797  }
1798 
1799  case 3:
1800  {
1801  // one could probably do something similar in 3d, but that seems
1802  // more complicated than one wants to write down. just go
1803  // through the list of possible starting vertices and check
1804  for (origin_vertex_of_cell = 0;
1805  origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
1806  ++origin_vertex_of_cell)
1807  if (std::count(starting_vertex_of_edge,
1808  starting_vertex_of_edge +
1810  cell_list[cell_index]
1811  .vertex_indices[origin_vertex_of_cell]) == dim)
1812  break;
1813  Assert(origin_vertex_of_cell <
1815  ExcInternalError());
1816 
1817  break;
1818  }
1819 
1820  default:
1821  Assert(false, ExcNotImplemented());
1822  }
1823 
1824  // now rotate raw_cells[cell_index] in such a way that its orientation
1825  // matches that of cell_list[cell_index]
1826  switch (dim)
1827  {
1828  case 2:
1829  {
1830  // in 2d, we can literally rotate the cell until its origin
1831  // matches the one that we have determined above should be
1832  // the origin vertex
1833  //
1834  // when doing a rotation, take into account the ordering of
1835  // vertices (not in clockwise or counter-clockwise sense)
1836  while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell)
1837  {
1838  const unsigned int tmp = raw_cells[cell_index].vertices[0];
1839  raw_cells[cell_index].vertices[0] =
1840  raw_cells[cell_index].vertices[1];
1841  raw_cells[cell_index].vertices[1] =
1842  raw_cells[cell_index].vertices[3];
1843  raw_cells[cell_index].vertices[3] =
1844  raw_cells[cell_index].vertices[2];
1845  raw_cells[cell_index].vertices[2] = tmp;
1846  }
1847  break;
1848  }
1849 
1850  case 3:
1851  {
1852  // in 3d, the situation is a bit more complicated. from above, we
1853  // now know which vertex is at the origin (because 3 edges
1854  // originate from it), but that still leaves 3 possible rotations
1855  // of the cube. the important realization is that we can choose
1856  // any of them: in all 3 rotations, all edges originate from the
1857  // one vertex, and that fixes the directions of all 12 edges in
1858  // the cube because these 3 cover all 3 equivalence classes!
1859  // consequently, we can select an arbitrary one among the
1860  // permutations -- for example the following ones:
1861  static const unsigned int cube_permutations[8][8] = {
1862  {0, 1, 2, 3, 4, 5, 6, 7},
1863  {1, 5, 3, 7, 0, 4, 2, 6},
1864  {2, 6, 0, 4, 3, 7, 1, 5},
1865  {3, 2, 1, 0, 7, 6, 5, 4},
1866  {4, 0, 6, 2, 5, 1, 7, 3},
1867  {5, 4, 7, 6, 1, 0, 3, 2},
1868  {6, 7, 4, 5, 2, 3, 0, 1},
1869  {7, 3, 5, 1, 6, 2, 4, 0}};
1870 
1871  unsigned int
1872  temp_vertex_indices[GeometryInfo<dim>::vertices_per_cell];
1873  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1874  temp_vertex_indices[v] =
1875  raw_cells[cell_index]
1876  .vertices[cube_permutations[origin_vertex_of_cell][v]];
1877  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1878  raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
1879 
1880  break;
1881  }
1882 
1883  default:
1884  {
1885  Assert(false, ExcNotImplemented());
1886  }
1887  }
1888  }
1889 
1890 
1896  template <int dim>
1897  void
1898  reorient(std::vector<CellData<dim>> &cells)
1899  {
1900  // first build the arrays that connect cells to edges and the other
1901  // way around
1902  std::vector<Edge<dim>> edge_list = build_edges(cells);
1903  std::vector<Cell<dim>> cell_list =
1904  build_cells_and_connect_edges(cells, edge_list);
1905 
1906  // then loop over all cells and start orienting parallel edge sets
1907  // of cells that still have non-oriented edges
1908  unsigned int next_cell_with_unoriented_edge = 0;
1909  while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
1910  cell_list, edge_list, next_cell_with_unoriented_edge)) !=
1912  {
1913  // see which edge sets are still not oriented
1914  //
1915  // we do not need to look at each edge because if we orient edge
1916  // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there
1917  // will be 3 other edges that are also oriented). there are only
1918  // dim independent sets of edges, so loop over these.
1919  //
1920  // we need to check whether each one of these starter edges may
1921  // already be oriented because the line (sheet) that connects
1922  // globally parallel edges may be self-intersecting in the
1923  // current cell
1924  for (unsigned int l = 0; l < dim; ++l)
1925  if (edge_list[cell_list[next_cell_with_unoriented_edge]
1927  .orientation_status == Edge<dim>::not_oriented)
1928  orient_one_set_of_parallel_edges(
1929  cell_list,
1930  edge_list,
1931  next_cell_with_unoriented_edge,
1933 
1934  // ensure that we have really oriented all edges now, not just
1935  // the starter edges
1936  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1937  Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
1938  .edge_indices[l]]
1939  .orientation_status != Edge<dim>::not_oriented,
1940  ExcInternalError());
1941  }
1942 
1943  // now that we have oriented all edges, we need to rotate cells
1944  // so that the edges point in the right direction with the now
1945  // rotated coordinate system
1946  for (unsigned int c = 0; c < cells.size(); ++c)
1947  rotate_cell(cell_list, edge_list, c, cells);
1948  }
1949 
1950 
1951  // overload of the function above for 1d -- there is nothing
1952  // to orient in that case
1953  void
1954  reorient(std::vector<CellData<1>> &)
1955  {}
1956  } // namespace
1957 
1958  template <int dim>
1959  void
1961  {
1962  Assert(cells.size() != 0,
1963  ExcMessage(
1964  "List of elements to orient must have at least one cell"));
1965 
1966  // there is nothing for us to do in 1d
1967  if (dim == 1)
1968  return;
1969 
1970  // check if grids are already consistent. if so, do
1971  // nothing. if not, then do the reordering
1972  if (!is_consistent(cells))
1973  try
1974  {
1975  reorient(cells);
1976  }
1977  catch (const ExcMeshNotOrientable &)
1978  {
1979  // the mesh is not orientable. this is acceptable if we are in 3d,
1980  // as class Triangulation knows how to handle this, but it is
1981  // not in 2d; in that case, re-throw the exception
1982  if (dim < 3)
1983  throw;
1984  }
1985  }
1986 
1987 
1988  // define some transformations
1989  namespace internal
1990  {
1991  template <int spacedim>
1992  class Shift
1993  {
1994  public:
1996  : shift(shift)
1997  {}
2000  {
2001  return p + shift;
2002  }
2003 
2004  private:
2006  };
2007 
2008 
2009  // Transformation to rotate around one of the cartesian axes.
2010  class Rotate3d
2011  {
2012  public:
2013  Rotate3d(const Tensor<1, 3, double> &axis, const double angle)
2014  : rotation_matrix(
2015  Physics::Transformations::Rotations::rotation_matrix_3d(axis,
2016  angle))
2017  {}
2018 
2019  Point<3>
2020  operator()(const Point<3> &p) const
2021  {
2022  return static_cast<Point<3>>(rotation_matrix * p);
2023  }
2024 
2025  private:
2027  };
2028 
2029 
2030  template <int spacedim>
2031  class Scale
2032  {
2033  public:
2034  explicit Scale(const double factor)
2035  : factor(factor)
2036  {}
2039  {
2040  return p * factor;
2041  }
2042 
2043  private:
2044  const double factor;
2045  };
2046  } // namespace internal
2047 
2048 
2049  template <int dim, int spacedim>
2050  void
2051  shift(const Tensor<1, spacedim> & shift_vector,
2053  {
2055  }
2056 
2057 
2058  template <int dim>
2059  void
2061  const double angle,
2063  {
2065  }
2066 
2067 
2068  template <int dim>
2069  void
2070  rotate(const double angle,
2071  const unsigned int axis,
2073  {
2074  Assert(axis < 3, ExcMessage("Invalid axis given!"));
2075 
2076  Tensor<1, 3, double> vector;
2077  vector[axis] = 1.;
2078 
2080  }
2081 
2082 
2083  template <int dim, int spacedim>
2084  void
2085  scale(const double scaling_factor,
2087  {
2088  Assert(scaling_factor > 0, ExcScalingFactorNotPositive(scaling_factor));
2090  }
2091 
2092 
2093  namespace internal
2094  {
2100  inline void
2102  const AffineConstraints<double> &constraints,
2103  Vector<double> & u)
2104  {
2105  const unsigned int n_dofs = S.n();
2106  const auto op = linear_operator(S);
2107  const auto SF = constrained_linear_operator(constraints, op);
2109  prec.initialize(S, 1.2);
2110 
2111  SolverControl control(n_dofs, 1.e-10, false, false);
2113  SolverCG<Vector<double>> solver(control, mem);
2114 
2115  Vector<double> f(n_dofs);
2116 
2117  const auto constrained_rhs =
2118  constrained_right_hand_side(constraints, op, f);
2119  solver.solve(SF, u, constrained_rhs, prec);
2120 
2121  constraints.distribute(u);
2122  }
2123  } // namespace internal
2124 
2125 
2126  // Implementation for dimensions except 1
2127  template <int dim>
2128  void
2129  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
2131  const Function<dim> * coefficient,
2132  const bool solve_for_absolute_positions)
2133  {
2134  if (dim == 1)
2135  Assert(false, ExcNotImplemented());
2136 
2137  // first provide everything that is needed for solving a Laplace
2138  // equation.
2139  FE_Q<dim> q1(1);
2140 
2141  DoFHandler<dim> dof_handler(triangulation);
2142  dof_handler.distribute_dofs(q1);
2143 
2144  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
2145  DoFTools::make_sparsity_pattern(dof_handler, dsp);
2146  dsp.compress();
2147 
2148  SparsityPattern sparsity_pattern;
2149  sparsity_pattern.copy_from(dsp);
2150  sparsity_pattern.compress();
2151 
2152  SparseMatrix<double> S(sparsity_pattern);
2153 
2154  QGauss<dim> quadrature(4);
2155 
2157  StaticMappingQ1<dim>::mapping, dof_handler, quadrature, S, coefficient);
2158 
2159  // set up the boundary values for the laplace problem
2160  std::array<AffineConstraints<double>, dim> constraints;
2161  typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
2162  new_points.end();
2163 
2164  // fill these maps using the data given by new_points
2165  for (const auto &cell : dof_handler.active_cell_iterators())
2166  {
2167  // loop over all vertices of the cell and see if it is listed in the map
2168  // given as first argument of the function
2169  for (const unsigned int vertex_no : cell->vertex_indices())
2170  {
2171  const unsigned int vertex_index = cell->vertex_index(vertex_no);
2172  const Point<dim> & vertex_point = cell->vertex(vertex_no);
2173 
2174  const typename std::map<unsigned int, Point<dim>>::const_iterator
2175  map_iter = new_points.find(vertex_index);
2176 
2177  if (map_iter != map_end)
2178  for (unsigned int i = 0; i < dim; ++i)
2179  {
2180  constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
2181  constraints[i].set_inhomogeneity(
2182  cell->vertex_dof_index(vertex_no, 0),
2183  (solve_for_absolute_positions ?
2184  map_iter->second(i) :
2185  map_iter->second(i) - vertex_point[i]));
2186  }
2187  }
2188  }
2189 
2190  for (unsigned int i = 0; i < dim; ++i)
2191  constraints[i].close();
2192 
2193  // solve the dim problems with different right hand sides.
2194  Vector<double> us[dim];
2195  for (unsigned int i = 0; i < dim; ++i)
2196  us[i].reinit(dof_handler.n_dofs());
2197 
2198  // solve linear systems in parallel
2199  Threads::TaskGroup<> tasks;
2200  for (unsigned int i = 0; i < dim; ++i)
2201  tasks +=
2202  Threads::new_task(&internal::laplace_solve, S, constraints[i], us[i]);
2203  tasks.join_all();
2204 
2205  // change the coordinates of the points of the triangulation
2206  // according to the computed values
2207  std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2208  for (const auto &cell : dof_handler.active_cell_iterators())
2209  for (const unsigned int vertex_no : cell->vertex_indices())
2210  if (vertex_touched[cell->vertex_index(vertex_no)] == false)
2211  {
2212  Point<dim> &v = cell->vertex(vertex_no);
2213 
2214  const types::global_dof_index dof_index =
2215  cell->vertex_dof_index(vertex_no, 0);
2216  for (unsigned int i = 0; i < dim; ++i)
2217  if (solve_for_absolute_positions)
2218  v(i) = us[i](dof_index);
2219  else
2220  v(i) += us[i](dof_index);
2221 
2222  vertex_touched[cell->vertex_index(vertex_no)] = true;
2223  }
2224  }
2225 
2226  template <int dim, int spacedim>
2227  std::map<unsigned int, Point<spacedim>>
2229  {
2230  std::map<unsigned int, Point<spacedim>> vertex_map;
2232  cell = tria.begin_active(),
2233  endc = tria.end();
2234  for (; cell != endc; ++cell)
2235  {
2236  for (unsigned int i : cell->face_indices())
2237  {
2238  const typename Triangulation<dim, spacedim>::face_iterator &face =
2239  cell->face(i);
2240  if (face->at_boundary())
2241  {
2242  for (unsigned j = 0; j < face->n_vertices(); ++j)
2243  {
2244  const Point<spacedim> &vertex = face->vertex(j);
2245  const unsigned int vertex_index = face->vertex_index(j);
2246  vertex_map[vertex_index] = vertex;
2247  }
2248  }
2249  }
2250  }
2251  return vertex_map;
2252  }
2253 
2258  template <int dim, int spacedim>
2259  void
2260  distort_random(const double factor,
2262  const bool keep_boundary,
2263  const unsigned int seed)
2264  {
2265  // if spacedim>dim we need to make sure that we perturb
2266  // points but keep them on
2267  // the manifold. however, this isn't implemented right now
2268  Assert(spacedim == dim, ExcNotImplemented());
2269 
2270 
2271  // find the smallest length of the
2272  // lines adjacent to the
2273  // vertex. take the initial value
2274  // to be larger than anything that
2275  // might be found: the diameter of
2276  // the triangulation, here
2277  // estimated by adding up the
2278  // diameters of the coarse grid
2279  // cells.
2280  double almost_infinite_length = 0;
2281  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
2282  triangulation.begin(0);
2283  cell != triangulation.end(0);
2284  ++cell)
2285  almost_infinite_length += cell->diameter();
2286 
2287  std::vector<double> minimal_length(triangulation.n_vertices(),
2288  almost_infinite_length);
2289 
2290  // also note if a vertex is at the boundary
2291  std::vector<bool> at_boundary(keep_boundary ? triangulation.n_vertices() :
2292  0,
2293  false);
2294  // for parallel::shared::Triangulation we need to work on all vertices,
2295  // not just the ones related to locally owned cells;
2296  const bool is_parallel_shared =
2298  &triangulation) != nullptr);
2299  for (const auto &cell : triangulation.active_cell_iterators())
2300  if (is_parallel_shared || cell->is_locally_owned())
2301  {
2302  if (dim > 1)
2303  {
2304  for (unsigned int i = 0; i < cell->n_lines(); ++i)
2305  {
2307  line = cell->line(i);
2308 
2309  if (keep_boundary && line->at_boundary())
2310  {
2311  at_boundary[line->vertex_index(0)] = true;
2312  at_boundary[line->vertex_index(1)] = true;
2313  }
2314 
2315  minimal_length[line->vertex_index(0)] =
2316  std::min(line->diameter(),
2317  minimal_length[line->vertex_index(0)]);
2318  minimal_length[line->vertex_index(1)] =
2319  std::min(line->diameter(),
2320  minimal_length[line->vertex_index(1)]);
2321  }
2322  }
2323  else // dim==1
2324  {
2325  if (keep_boundary)
2326  for (unsigned int vertex = 0; vertex < 2; ++vertex)
2327  if (cell->at_boundary(vertex) == true)
2328  at_boundary[cell->vertex_index(vertex)] = true;
2329 
2330  minimal_length[cell->vertex_index(0)] =
2331  std::min(cell->diameter(),
2332  minimal_length[cell->vertex_index(0)]);
2333  minimal_length[cell->vertex_index(1)] =
2334  std::min(cell->diameter(),
2335  minimal_length[cell->vertex_index(1)]);
2336  }
2337  }
2338 
2339  // create a random number generator for the interval [-1,1]
2340  boost::random::mt19937 rng(seed);
2341  boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
2342 
2343  // If the triangulation is distributed, we need to
2344  // exchange the moved vertices across mpi processes
2345  if (auto distributed_triangulation =
2347  &triangulation))
2348  {
2349  const std::vector<bool> locally_owned_vertices =
2351  std::vector<bool> vertex_moved(triangulation.n_vertices(), false);
2352 
2353  // Next move vertices on locally owned cells
2354  for (const auto &cell : triangulation.active_cell_iterators())
2355  if (cell->is_locally_owned())
2356  {
2357  for (const unsigned int vertex_no : cell->vertex_indices())
2358  {
2359  const unsigned global_vertex_no =
2360  cell->vertex_index(vertex_no);
2361 
2362  // ignore this vertex if we shall keep the boundary and
2363  // this vertex *is* at the boundary, if it is already moved
2364  // or if another process moves this vertex
2365  if ((keep_boundary && at_boundary[global_vertex_no]) ||
2366  vertex_moved[global_vertex_no] ||
2367  !locally_owned_vertices[global_vertex_no])
2368  continue;
2369 
2370  // first compute a random shift vector
2371  Point<spacedim> shift_vector;
2372  for (unsigned int d = 0; d < spacedim; ++d)
2373  shift_vector(d) = uniform_distribution(rng);
2374 
2375  shift_vector *= factor * minimal_length[global_vertex_no] /
2376  std::sqrt(shift_vector.square());
2377 
2378  // finally move the vertex
2379  cell->vertex(vertex_no) += shift_vector;
2380  vertex_moved[global_vertex_no] = true;
2381  }
2382  }
2383 
2384  distributed_triangulation->communicate_locally_moved_vertices(
2385  locally_owned_vertices);
2386  }
2387  else
2388  // if this is a sequential triangulation, we could in principle
2389  // use the algorithm above, but we'll use an algorithm that we used
2390  // before the parallel::distributed::Triangulation was introduced
2391  // in order to preserve backward compatibility
2392  {
2393  // loop over all vertices and compute their new locations
2394  const unsigned int n_vertices = triangulation.n_vertices();
2395  std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
2396  const std::vector<Point<spacedim>> &old_vertex_locations =
2397  triangulation.get_vertices();
2398 
2399  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2400  {
2401  // ignore this vertex if we will keep the boundary and
2402  // this vertex *is* at the boundary
2403  if (keep_boundary && at_boundary[vertex])
2404  new_vertex_locations[vertex] = old_vertex_locations[vertex];
2405  else
2406  {
2407  // compute a random shift vector
2408  Point<spacedim> shift_vector;
2409  for (unsigned int d = 0; d < spacedim; ++d)
2410  shift_vector(d) = uniform_distribution(rng);
2411 
2412  shift_vector *= factor * minimal_length[vertex] /
2413  std::sqrt(shift_vector.square());
2414 
2415  // record new vertex location
2416  new_vertex_locations[vertex] =
2417  old_vertex_locations[vertex] + shift_vector;
2418  }
2419  }
2420 
2421  // now do the actual move of the vertices
2422  for (const auto &cell : triangulation.active_cell_iterators())
2423  for (const unsigned int vertex_no : cell->vertex_indices())
2424  cell->vertex(vertex_no) =
2425  new_vertex_locations[cell->vertex_index(vertex_no)];
2426  }
2427 
2428  // Correct hanging nodes if necessary
2429  if (dim >= 2)
2430  {
2431  // We do the same as in GridTools::transform
2432  //
2433  // exclude hanging nodes at the boundaries of artificial cells:
2434  // these may belong to ghost cells for which we know the exact
2435  // location of vertices, whereas the artificial cell may or may
2436  // not be further refined, and so we cannot know whether
2437  // the location of the hanging node is correct or not
2439  cell = triangulation.begin_active(),
2440  endc = triangulation.end();
2441  for (; cell != endc; ++cell)
2442  if (!cell->is_artificial())
2443  for (const unsigned int face : cell->face_indices())
2444  if (cell->face(face)->has_children() &&
2445  !cell->face(face)->at_boundary())
2446  {
2447  // this face has hanging nodes
2448  if (dim == 2)
2449  cell->face(face)->child(0)->vertex(1) =
2450  (cell->face(face)->vertex(0) +
2451  cell->face(face)->vertex(1)) /
2452  2;
2453  else if (dim == 3)
2454  {
2455  cell->face(face)->child(0)->vertex(1) =
2456  .5 * (cell->face(face)->vertex(0) +
2457  cell->face(face)->vertex(1));
2458  cell->face(face)->child(0)->vertex(2) =
2459  .5 * (cell->face(face)->vertex(0) +
2460  cell->face(face)->vertex(2));
2461  cell->face(face)->child(1)->vertex(3) =
2462  .5 * (cell->face(face)->vertex(1) +
2463  cell->face(face)->vertex(3));
2464  cell->face(face)->child(2)->vertex(3) =
2465  .5 * (cell->face(face)->vertex(2) +
2466  cell->face(face)->vertex(3));
2467 
2468  // center of the face
2469  cell->face(face)->child(0)->vertex(3) =
2470  .25 * (cell->face(face)->vertex(0) +
2471  cell->face(face)->vertex(1) +
2472  cell->face(face)->vertex(2) +
2473  cell->face(face)->vertex(3));
2474  }
2475  }
2476  }
2477  }
2478 
2479 
2480 
2481  template <int dim, template <int, int> class MeshType, int spacedim>
2482  unsigned int
2483  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
2484  const Point<spacedim> & p,
2485  const std::vector<bool> & marked_vertices)
2486  {
2487  // first get the underlying triangulation from the mesh and determine
2488  // vertices and used vertices
2490 
2491  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
2492 
2493  Assert(tria.get_vertices().size() == marked_vertices.size() ||
2494  marked_vertices.size() == 0,
2496  marked_vertices.size()));
2497 
2498  // marked_vertices is expected to be a subset of used_vertices. Thus,
2499  // comparing the range marked_vertices.begin() to marked_vertices.end() with
2500  // the range used_vertices.begin() to used_vertices.end() the element in the
2501  // second range must be valid if the element in the first range is valid.
2502  Assert(
2503  marked_vertices.size() == 0 ||
2504  std::equal(marked_vertices.begin(),
2505  marked_vertices.end(),
2506  tria.get_used_vertices().begin(),
2507  [](bool p, bool q) { return !p || q; }),
2508  ExcMessage(
2509  "marked_vertices should be a subset of used vertices in the triangulation "
2510  "but marked_vertices contains one or more vertices that are not used vertices!"));
2511 
2512  // If marked_indices is empty, consider all used_vertices for finding the
2513  // closest vertex to the point. Otherwise, marked_indices is used.
2514  const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
2516  marked_vertices;
2517 
2518  // At the beginning, the first used vertex is considered to be the closest
2519  // one.
2520  std::vector<bool>::const_iterator first =
2521  std::find(vertices_to_use.begin(), vertices_to_use.end(), true);
2522 
2523  // Assert that at least one vertex is actually used
2524  Assert(first != vertices_to_use.end(), ExcInternalError());
2525 
2526  unsigned int best_vertex = std::distance(vertices_to_use.begin(), first);
2527  double best_dist = (p - vertices[best_vertex]).norm_square();
2528 
2529  // For all remaining vertices, test
2530  // whether they are any closer
2531  for (unsigned int j = best_vertex + 1; j < vertices.size(); ++j)
2532  if (vertices_to_use[j])
2533  {
2534  const double dist = (p - vertices[j]).norm_square();
2535  if (dist < best_dist)
2536  {
2537  best_vertex = j;
2538  best_dist = dist;
2539  }
2540  }
2541 
2542  return best_vertex;
2543  }
2544 
2545 
2546 
2547  template <int dim, template <int, int> class MeshType, int spacedim>
2548  unsigned int
2550  const MeshType<dim, spacedim> &mesh,
2551  const Point<spacedim> & p,
2552  const std::vector<bool> & marked_vertices)
2553  {
2554  // Take a shortcut in the simple case.
2555  if (mapping.preserves_vertex_locations() == true)
2556  return find_closest_vertex(mesh, p, marked_vertices);
2557 
2558  // first get the underlying triangulation from the mesh and determine
2559  // vertices and used vertices
2561 
2562  auto vertices = extract_used_vertices(tria, mapping);
2563 
2564  Assert(tria.get_vertices().size() == marked_vertices.size() ||
2565  marked_vertices.size() == 0,
2567  marked_vertices.size()));
2568 
2569  // marked_vertices is expected to be a subset of used_vertices. Thus,
2570  // comparing the range marked_vertices.begin() to marked_vertices.end()
2571  // with the range used_vertices.begin() to used_vertices.end() the element
2572  // in the second range must be valid if the element in the first range is
2573  // valid.
2574  Assert(
2575  marked_vertices.size() == 0 ||
2576  std::equal(marked_vertices.begin(),
2577  marked_vertices.end(),
2578  tria.get_used_vertices().begin(),
2579  [](bool p, bool q) { return !p || q; }),
2580  ExcMessage(
2581  "marked_vertices should be a subset of used vertices in the triangulation "
2582  "but marked_vertices contains one or more vertices that are not used vertices!"));
2583 
2584  // Remove from the map unwanted elements.
2585  if (marked_vertices.size() != 0)
2586  for (auto it = vertices.begin(); it != vertices.end();)
2587  {
2588  if (marked_vertices[it->first] == false)
2589  {
2590  it = vertices.erase(it);
2591  }
2592  else
2593  {
2594  ++it;
2595  }
2596  }
2597 
2598  return find_closest_vertex(vertices, p);
2599  }
2600 
2601 
2602 
2603  template <int dim, template <int, int> class MeshType, int spacedim>
2604 #ifndef _MSC_VER
2605  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
2606 #else
2607  std::vector<
2608  typename ::internal::
2609  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2610 #endif
2611  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &mesh,
2612  const unsigned int vertex)
2613  {
2614  // make sure that the given vertex is
2615  // an active vertex of the underlying
2616  // triangulation
2617  AssertIndexRange(vertex, mesh.get_triangulation().n_vertices());
2618  Assert(mesh.get_triangulation().get_used_vertices()[vertex],
2619  ExcVertexNotUsed(vertex));
2620 
2621  // use a set instead of a vector
2622  // to ensure that cells are inserted only
2623  // once
2624  std::set<typename ::internal::
2625  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2627 
2628  // go through all active cells and look if the vertex is part of that cell
2629  //
2630  // in 1d, this is all we need to care about. in 2d/3d we also need to worry
2631  // that the vertex might be a hanging node on a face or edge of a cell; in
2632  // this case, we would want to add those cells as well on whose faces the
2633  // vertex is located but for which it is not a vertex itself.
2634  //
2635  // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging
2636  // node can only be in the middle of a face and we can query the neighboring
2637  // cell from the current cell. on the other hand, in 3d a hanging node
2638  // vertex can also be on an edge but there can be many other cells on
2639  // this edge and we can not access them from the cell we are currently
2640  // on.
2641  //
2642  // so, in the 3d case, if we run the algorithm as in 2d, we catch all
2643  // those cells for which the vertex we seek is on a *subface*, but we
2644  // miss the case of cells for which the vertex we seek is on a
2645  // sub-edge for which there is no corresponding sub-face (because the
2646  // immediate neighbor behind this face is not refined), see for example
2647  // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we
2648  // haven't yet found the vertex for the current cell we also need to
2649  // look at the mid-points of edges
2650  //
2651  // as a final note, deciding whether a neighbor is actually coarser is
2652  // simple in the case of isotropic refinement (we just need to look at
2653  // the level of the current and the neighboring cell). however, this
2654  // isn't so simple if we have used anisotropic refinement since then
2655  // the level of a cell is not indicative of whether it is coarser or
2656  // not than the current cell. ultimately, we want to add all cells on
2657  // which the vertex is, independent of whether they are coarser or
2658  // finer and so in the 2d case below we simply add *any* *active* neighbor.
2659  // in the worst case, we add cells multiple times to the adjacent_cells
2660  // list, but std::set throws out those cells already entered
2661  for (const auto &cell : mesh.active_cell_iterators())
2662  {
2663  for (const unsigned int v : cell->vertex_indices())
2664  if (cell->vertex_index(v) == vertex)
2665  {
2666  // OK, we found a cell that contains
2667  // the given vertex. We add it
2668  // to the list.
2669  adjacent_cells.insert(cell);
2670 
2671  // as explained above, in 2+d we need to check whether
2672  // this vertex is on a face behind which there is a
2673  // (possibly) coarser neighbor. if this is the case,
2674  // then we need to also add this neighbor
2675  if (dim >= 2)
2676  for (const auto face :
2677  cell->reference_cell().faces_for_given_vertex(v))
2678  if (!cell->at_boundary(face) &&
2679  cell->neighbor(face)->is_active())
2680  {
2681  // there is a (possibly) coarser cell behind a
2682  // face to which the vertex belongs. the
2683  // vertex we are looking at is then either a
2684  // vertex of that coarser neighbor, or it is a
2685  // hanging node on one of the faces of that
2686  // cell. in either case, it is adjacent to the
2687  // vertex, so add it to the list as well (if
2688  // the cell was already in the list then the
2689  // std::set makes sure that we get it only
2690  // once)
2691  adjacent_cells.insert(cell->neighbor(face));
2692  }
2693 
2694  // in any case, we have found a cell, so go to the next cell
2695  goto next_cell;
2696  }
2697 
2698  // in 3d also loop over the edges
2699  if (dim >= 3)
2700  {
2701  for (unsigned int e = 0; e < cell->n_lines(); ++e)
2702  if (cell->line(e)->has_children())
2703  // the only place where this vertex could have been
2704  // hiding is on the mid-edge point of the edge we
2705  // are looking at
2706  if (cell->line(e)->child(0)->vertex_index(1) == vertex)
2707  {
2708  adjacent_cells.insert(cell);
2709 
2710  // jump out of this tangle of nested loops
2711  goto next_cell;
2712  }
2713  }
2714 
2715  // in more than 3d we would probably have to do the same as
2716  // above also for even lower-dimensional objects
2717  Assert(dim <= 3, ExcNotImplemented());
2718 
2719  // move on to the next cell if we have found the
2720  // vertex on the current one
2721  next_cell:;
2722  }
2723 
2724  // if this was an active vertex then there needs to have been
2725  // at least one cell to which it is adjacent!
2726  Assert(adjacent_cells.size() > 0, ExcInternalError());
2727 
2728  // return the result as a vector, rather than the set we built above
2729  return std::vector<
2730  typename ::internal::
2731  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
2732  adjacent_cells.begin(), adjacent_cells.end());
2733  }
2734 
2735 
2736 
2737  template <int dim, int spacedim>
2738  std::vector<std::vector<Tensor<1, spacedim>>>
2740  const Triangulation<dim, spacedim> &mesh,
2741  const std::vector<
2743  &vertex_to_cells)
2744  {
2745  const std::vector<Point<spacedim>> &vertices = mesh.get_vertices();
2746  const unsigned int n_vertices = vertex_to_cells.size();
2747 
2748  AssertDimension(vertices.size(), n_vertices);
2749 
2750 
2751  std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
2752  n_vertices);
2753  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2754  if (mesh.vertex_used(vertex))
2755  {
2756  const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
2757  vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
2758 
2759  typename std::set<typename Triangulation<dim, spacedim>::
2760  active_cell_iterator>::iterator it =
2761  vertex_to_cells[vertex].begin();
2762  for (unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
2763  {
2764  vertex_to_cell_centers[vertex][cell] =
2765  (*it)->center() - vertices[vertex];
2766  vertex_to_cell_centers[vertex][cell] /=
2767  vertex_to_cell_centers[vertex][cell].norm();
2768  }
2769  }
2770  return vertex_to_cell_centers;
2771  }
2772 
2773 
2774  namespace internal
2775  {
2776  template <int spacedim>
2777  bool
2779  const unsigned int a,
2780  const unsigned int b,
2781  const Tensor<1, spacedim> & point_direction,
2782  const std::vector<Tensor<1, spacedim>> &center_directions)
2783  {
2784  const double scalar_product_a = center_directions[a] * point_direction;
2785  const double scalar_product_b = center_directions[b] * point_direction;
2786 
2787  // The function is supposed to return if a is before b. We are looking
2788  // for the alignment of point direction and center direction, therefore
2789  // return if the scalar product of a is larger.
2790  return (scalar_product_a > scalar_product_b);
2791  }
2792  } // namespace internal
2793 
2794  template <int dim, template <int, int> class MeshType, int spacedim>
2795 #ifndef _MSC_VER
2796  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
2797 #else
2798  std::pair<typename ::internal::
2799  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
2800  Point<dim>>
2801 #endif
2803  const Mapping<dim, spacedim> & mapping,
2804  const MeshType<dim, spacedim> &mesh,
2805  const Point<spacedim> & p,
2806  const std::vector<
2807  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
2808  & vertex_to_cells,
2809  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
2810  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
2811  const std::vector<bool> & marked_vertices,
2812  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree,
2813  const double tolerance,
2814  const RTree<
2815  std::pair<BoundingBox<spacedim>,
2817  *relevant_cell_bounding_boxes_rtree)
2818  {
2819  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2820  Point<dim>>
2821  cell_and_position;
2822  cell_and_position.first = mesh.end();
2823 
2824  // To handle points at the border we keep track of points which are close to
2825  // the unit cell:
2826  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2827  Point<dim>>
2828  cell_and_position_approx;
2829 
2830  if (relevant_cell_bounding_boxes_rtree != nullptr &&
2831  !relevant_cell_bounding_boxes_rtree->empty())
2832  {
2833  // create a bounding box around point p with 2*tolerance as side length.
2834  auto p1 = p;
2835  auto p2 = p;
2836 
2837  for (unsigned int d = 0; d < spacedim; ++d)
2838  {
2839  p1[d] = p1[d] - tolerance;
2840  p2[d] = p2[d] + tolerance;
2841  }
2842 
2843  BoundingBox<spacedim> bb({p1, p2});
2844 
2845  if (relevant_cell_bounding_boxes_rtree->qbegin(
2846  boost::geometry::index::intersects(bb)) ==
2847  relevant_cell_bounding_boxes_rtree->qend())
2848  return cell_and_position;
2849  }
2850 
2851  bool found_cell = false;
2852  bool approx_cell = false;
2853 
2854  unsigned int closest_vertex_index = 0;
2855  Tensor<1, spacedim> vertex_to_point;
2856  auto current_cell = cell_hint;
2857 
2858  while (found_cell == false)
2859  {
2860  // First look at the vertices of the cell cell_hint. If it's an
2861  // invalid cell, then query for the closest global vertex
2862  if (current_cell.state() == IteratorState::valid)
2863  {
2864  const auto cell_vertices = mapping.get_vertices(current_cell);
2865  const unsigned int closest_vertex =
2866  find_closest_vertex_of_cell<dim, spacedim>(current_cell,
2867  p,
2868  mapping);
2869  vertex_to_point = p - cell_vertices[closest_vertex];
2870  closest_vertex_index = current_cell->vertex_index(closest_vertex);
2871  }
2872  else
2873  {
2874  if (!used_vertices_rtree.empty())
2875  {
2876  // If we have an rtree at our disposal, use it.
2877  using ValueType = std::pair<Point<spacedim>, unsigned int>;
2878  std::function<bool(const ValueType &)> marked;
2879  if (marked_vertices.size() == mesh.n_vertices())
2880  marked = [&marked_vertices](const ValueType &value) -> bool {
2881  return marked_vertices[value.second];
2882  };
2883  else
2884  marked = [](const ValueType &) -> bool { return true; };
2885 
2886  std::vector<std::pair<Point<spacedim>, unsigned int>> res;
2887  used_vertices_rtree.query(
2888  boost::geometry::index::nearest(p, 1) &&
2889  boost::geometry::index::satisfies(marked),
2890  std::back_inserter(res));
2891 
2892  // We should have one and only one result
2893  AssertDimension(res.size(), 1);
2894  closest_vertex_index = res[0].second;
2895  }
2896  else
2897  {
2898  closest_vertex_index = GridTools::find_closest_vertex(
2899  mapping, mesh, p, marked_vertices);
2900  }
2901  vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
2902  }
2903 
2904  const double vertex_point_norm = vertex_to_point.norm();
2905  if (vertex_point_norm > 0)
2906  vertex_to_point /= vertex_point_norm;
2907 
2908  const unsigned int n_neighbor_cells =
2909  vertex_to_cells[closest_vertex_index].size();
2910 
2911  // Create a corresponding map of vectors from vertex to cell center
2912  std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
2913 
2914  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
2915  neighbor_permutation[i] = i;
2916 
2917  auto comp = [&](const unsigned int a, const unsigned int b) -> bool {
2918  return internal::compare_point_association<spacedim>(
2919  a,
2920  b,
2921  vertex_to_point,
2922  vertex_to_cell_centers[closest_vertex_index]);
2923  };
2924 
2925  std::sort(neighbor_permutation.begin(),
2926  neighbor_permutation.end(),
2927  comp);
2928  // It is possible the vertex is close
2929  // to an edge, thus we add a tolerance
2930  // to keep also the "best" cell
2931  double best_distance = tolerance;
2932 
2933  // Search all of the cells adjacent to the closest vertex of the cell
2934  // hint. Most likely we will find the point in them.
2935  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
2936  {
2937  try
2938  {
2939  auto cell = vertex_to_cells[closest_vertex_index].begin();
2940  std::advance(cell, neighbor_permutation[i]);
2941 
2942  if (!(*cell)->is_artificial())
2943  {
2944  const Point<dim> p_unit =
2945  mapping.transform_real_to_unit_cell(*cell, p);
2947  tolerance))
2948  {
2949  cell_and_position.first = *cell;
2950  cell_and_position.second = p_unit;
2951  found_cell = true;
2952  approx_cell = false;
2953  break;
2954  }
2955  // The point is not inside this cell: checking how far
2956  // outside it is and whether we want to use this cell as a
2957  // backup if we can't find a cell within which the point
2958  // lies.
2959  const double dist =
2961  if (dist < best_distance)
2962  {
2963  best_distance = dist;
2964  cell_and_position_approx.first = *cell;
2965  cell_and_position_approx.second = p_unit;
2966  approx_cell = true;
2967  }
2968  }
2969  }
2970  catch (typename Mapping<dim>::ExcTransformationFailed &)
2971  {}
2972  }
2973 
2974  if (found_cell == true)
2975  return cell_and_position;
2976  else if (approx_cell == true)
2977  return cell_and_position_approx;
2978 
2979  // The first time around, we check for vertices in the hint_cell. If
2980  // that does not work, we set the cell iterator to an invalid one, and
2981  // look for a global vertex close to the point. If that does not work,
2982  // we are in trouble, and just throw an exception.
2983  //
2984  // If we got here, then we did not find the point. If the
2985  // current_cell.state() here is not IteratorState::valid, it means that
2986  // the user did not provide a hint_cell, and at the beginning of the
2987  // while loop we performed an actual global search on the mesh
2988  // vertices. Not finding the point then means the point is outside the
2989  // domain, or that we've had problems with the algorithm above. Try as a
2990  // last resort the other (simpler) algorithm.
2991  if (current_cell.state() != IteratorState::valid)
2993  mapping, mesh, p, marked_vertices, tolerance);
2994 
2995  current_cell = typename MeshType<dim, spacedim>::active_cell_iterator();
2996  }
2997  return cell_and_position;
2998  }
2999 
3000 
3001 
3002  template <int dim, int spacedim>
3003  unsigned int
3006  const Point<spacedim> & position,
3007  const Mapping<dim, spacedim> & mapping)
3008  {
3009  const auto vertices = mapping.get_vertices(cell);
3010  double minimum_distance = position.distance_square(vertices[0]);
3011  unsigned int closest_vertex = 0;
3012 
3013  for (unsigned int v = 1; v < cell->n_vertices(); ++v)
3014  {
3015  const double vertex_distance = position.distance_square(vertices[v]);
3016  if (vertex_distance < minimum_distance)
3017  {
3018  closest_vertex = v;
3019  minimum_distance = vertex_distance;
3020  }
3021  }
3022  return closest_vertex;
3023  }
3024 
3025 
3026 
3027  namespace internal
3028  {
3029  namespace BoundingBoxPredicate
3030  {
3031  template <class MeshType>
3032  std::tuple<BoundingBox<MeshType::space_dimension>, bool>
3034  const typename MeshType::cell_iterator &parent_cell,
3035  const std::function<
3036  bool(const typename MeshType::active_cell_iterator &)> &predicate)
3037  {
3038  bool has_predicate =
3039  false; // Start assuming there's no cells with predicate inside
3040  std::vector<typename MeshType::active_cell_iterator> active_cells;
3041  if (parent_cell->is_active())
3042  active_cells = {parent_cell};
3043  else
3044  // Finding all active cells descendants of the current one (or the
3045  // current one if it is active)
3046  active_cells = get_active_child_cells<MeshType>(parent_cell);
3047 
3048  const unsigned int spacedim = MeshType::space_dimension;
3049 
3050  // Looking for the first active cell which has the property predicate
3051  unsigned int i = 0;
3052  while (i < active_cells.size() && !predicate(active_cells[i]))
3053  ++i;
3054 
3055  // No active cells or no active cells with property
3056  if (active_cells.size() == 0 || i == active_cells.size())
3057  {
3058  BoundingBox<spacedim> bbox;
3059  return std::make_tuple(bbox, has_predicate);
3060  }
3061 
3062  // The two boundary points defining the boundary box
3063  Point<spacedim> maxp = active_cells[i]->vertex(0);
3064  Point<spacedim> minp = active_cells[i]->vertex(0);
3065 
3066  for (; i < active_cells.size(); ++i)
3067  if (predicate(active_cells[i]))
3068  for (const unsigned int v : active_cells[i]->vertex_indices())
3069  for (unsigned int d = 0; d < spacedim; ++d)
3070  {
3071  minp[d] = std::min(minp[d], active_cells[i]->vertex(v)[d]);
3072  maxp[d] = std::max(maxp[d], active_cells[i]->vertex(v)[d]);
3073  }
3074 
3075  has_predicate = true;
3076  BoundingBox<spacedim> bbox(std::make_pair(minp, maxp));
3077  return std::make_tuple(bbox, has_predicate);
3078  }
3079  } // namespace BoundingBoxPredicate
3080  } // namespace internal
3081 
3082 
3083 
3084  template <class MeshType>
3085  std::vector<BoundingBox<MeshType::space_dimension>>
3087  const MeshType &mesh,
3088  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3089  & predicate,
3090  const unsigned int refinement_level,
3091  const bool allow_merge,
3092  const unsigned int max_boxes)
3093  {
3094  // Algorithm brief description: begin with creating bounding boxes of all
3095  // cells at refinement_level (and coarser levels if there are active cells)
3096  // which have the predicate property. These are then merged
3097 
3098  Assert(
3099  refinement_level <= mesh.n_levels(),
3100  ExcMessage(
3101  "Error: refinement level is higher then total levels in the triangulation!"));
3102 
3103  const unsigned int spacedim = MeshType::space_dimension;
3104  std::vector<BoundingBox<spacedim>> bounding_boxes;
3105 
3106  // Creating a bounding box for all active cell on coarser level
3107 
3108  for (unsigned int i = 0; i < refinement_level; ++i)
3109  for (const typename MeshType::cell_iterator &cell :
3110  mesh.active_cell_iterators_on_level(i))
3111  {
3112  bool has_predicate = false;
3113  BoundingBox<spacedim> bbox;
3114  std::tie(bbox, has_predicate) =
3116  MeshType>(cell, predicate);
3117  if (has_predicate)
3118  bounding_boxes.push_back(bbox);
3119  }
3120 
3121  // Creating a Bounding Box for all cells on the chosen refinement_level
3122  for (const typename MeshType::cell_iterator &cell :
3123  mesh.cell_iterators_on_level(refinement_level))
3124  {
3125  bool has_predicate = false;
3126  BoundingBox<spacedim> bbox;
3127  std::tie(bbox, has_predicate) =
3129  MeshType>(cell, predicate);
3130  if (has_predicate)
3131  bounding_boxes.push_back(bbox);
3132  }
3133 
3134  if (!allow_merge)
3135  // If merging is not requested return the created bounding_boxes
3136  return bounding_boxes;
3137  else
3138  {
3139  // Merging part of the algorithm
3140  // Part 1: merging neighbors
3141  // This array stores the indices of arrays we have already merged
3142  std::vector<unsigned int> merged_boxes_idx;
3143  bool found_neighbors = true;
3144 
3145  // We merge only neighbors which can be expressed by a single bounding
3146  // box e.g. in 1d [0,1] and [1,2] can be described with [0,2] without
3147  // losing anything
3148  while (found_neighbors)
3149  {
3150  found_neighbors = false;
3151  for (unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
3152  {
3153  if (std::find(merged_boxes_idx.begin(),
3154  merged_boxes_idx.end(),
3155  i) == merged_boxes_idx.end())
3156  for (unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
3157  if (std::find(merged_boxes_idx.begin(),
3158  merged_boxes_idx.end(),
3159  j) == merged_boxes_idx.end() &&
3160  bounding_boxes[i].get_neighbor_type(
3161  bounding_boxes[j]) ==
3163  {
3164  bounding_boxes[i].merge_with(bounding_boxes[j]);
3165  merged_boxes_idx.push_back(j);
3166  found_neighbors = true;
3167  }
3168  }
3169  }
3170 
3171  // Copying the merged boxes into merged_b_boxes
3172  std::vector<BoundingBox<spacedim>> merged_b_boxes;
3173  for (unsigned int i = 0; i < bounding_boxes.size(); ++i)
3174  if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
3175  merged_boxes_idx.end())
3176  merged_b_boxes.push_back(bounding_boxes[i]);
3177 
3178  // Part 2: if there are too many bounding boxes, merging smaller boxes
3179  // This has sense only in dimension 2 or greater, since in dimension 1,
3180  // neighboring intervals can always be merged without problems
3181  if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
3182  {
3183  std::vector<double> volumes;
3184  for (unsigned int i = 0; i < merged_b_boxes.size(); ++i)
3185  volumes.push_back(merged_b_boxes[i].volume());
3186 
3187  while (merged_b_boxes.size() > max_boxes)
3188  {
3189  unsigned int min_idx =
3190  std::min_element(volumes.begin(), volumes.end()) -
3191  volumes.begin();
3192  volumes.erase(volumes.begin() + min_idx);
3193  // Finding a neighbor
3194  bool not_removed = true;
3195  for (unsigned int i = 0;
3196  i < merged_b_boxes.size() && not_removed;
3197  ++i)
3198  // We merge boxes if we have "attached" or "mergeable"
3199  // neighbors, even though mergeable should be dealt with in
3200  // Part 1
3201  if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
3202  merged_b_boxes[min_idx]) ==
3204  merged_b_boxes[i].get_neighbor_type(
3205  merged_b_boxes[min_idx]) ==
3207  {
3208  merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
3209  merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
3210  not_removed = false;
3211  }
3212  Assert(!not_removed,
3213  ExcMessage("Error: couldn't merge bounding boxes!"));
3214  }
3215  }
3216  Assert(merged_b_boxes.size() <= max_boxes,
3217  ExcMessage(
3218  "Error: couldn't reach target number of bounding boxes!"));
3219  return merged_b_boxes;
3220  }
3221  }
3222 
3223 
3224 
3225  template <int spacedim>
3226 #ifndef DOXYGEN
3227  std::tuple<std::vector<std::vector<unsigned int>>,
3228  std::map<unsigned int, unsigned int>,
3229  std::map<unsigned int, std::vector<unsigned int>>>
3230 #else
3231  return_type
3232 #endif
3234  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3235  const std::vector<Point<spacedim>> & points)
3236  {
3237  unsigned int n_procs = global_bboxes.size();
3238  std::vector<std::vector<unsigned int>> point_owners(n_procs);
3239  std::map<unsigned int, unsigned int> map_owners_found;
3240  std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3241 
3242  unsigned int n_points = points.size();
3243  for (unsigned int pt = 0; pt < n_points; ++pt)
3244  {
3245  // Keep track of how many processes we guess to own the point
3246  std::vector<unsigned int> owners_found;
3247  // Check in which other processes the point might be
3248  for (unsigned int rk = 0; rk < n_procs; ++rk)
3249  {
3250  for (const BoundingBox<spacedim> &bbox : global_bboxes[rk])
3251  if (bbox.point_inside(points[pt]))
3252  {
3253  point_owners[rk].emplace_back(pt);
3254  owners_found.emplace_back(rk);
3255  break; // We can check now the next process
3256  }
3257  }
3258  Assert(owners_found.size() > 0,
3259  ExcMessage("No owners found for the point " +
3260  std::to_string(pt)));
3261  if (owners_found.size() == 1)
3262  map_owners_found[pt] = owners_found[0];
3263  else
3264  // Multiple owners
3265  map_owners_guessed[pt] = owners_found;
3266  }
3267 
3268  return std::make_tuple(std::move(point_owners),
3269  std::move(map_owners_found),
3270  std::move(map_owners_guessed));
3271  }
3272 
3273  template <int spacedim>
3274 #ifndef DOXYGEN
3275  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
3276  std::map<unsigned int, unsigned int>,
3277  std::map<unsigned int, std::vector<unsigned int>>>
3278 #else
3279  return_type
3280 #endif
3282  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
3283  const std::vector<Point<spacedim>> & points)
3284  {
3285  std::map<unsigned int, std::vector<unsigned int>> point_owners;
3286  std::map<unsigned int, unsigned int> map_owners_found;
3287  std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3288  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> search_result;
3289 
3290  unsigned int n_points = points.size();
3291  for (unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
3292  {
3293  search_result.clear(); // clearing last output
3294 
3295  // Running tree search
3296  covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
3297  std::back_inserter(search_result));
3298 
3299  // Keep track of how many processes we guess to own the point
3300  std::set<unsigned int> owners_found;
3301  // Check in which other processes the point might be
3302  for (const auto &rank_bbox : search_result)
3303  {
3304  // Try to add the owner to the owners found,
3305  // and check if it was already present
3306  const bool pt_inserted = owners_found.insert(pt_n).second;
3307  if (pt_inserted)
3308  point_owners[rank_bbox.second].emplace_back(pt_n);
3309  }
3310  Assert(owners_found.size() > 0,
3311  ExcMessage("No owners found for the point " +
3312  std::to_string(pt_n)));
3313  if (owners_found.size() == 1)
3314  map_owners_found[pt_n] = *owners_found.begin();
3315  else
3316  // Multiple owners
3317  std::copy(owners_found.begin(),
3318  owners_found.end(),
3319  std::back_inserter(map_owners_guessed[pt_n]));
3320  }
3321 
3322  return std::make_tuple(std::move(point_owners),
3323  std::move(map_owners_found),
3324  std::move(map_owners_guessed));
3325  }
3326 
3327 
3328  template <int dim, int spacedim>
3329  std::vector<
3330  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3332  {
3333  std::vector<
3334  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3335  vertex_to_cell_map(triangulation.n_vertices());
3337  cell = triangulation.begin_active(),
3338  endc = triangulation.end();
3339  for (; cell != endc; ++cell)
3340  for (const unsigned int i : cell->vertex_indices())
3341  vertex_to_cell_map[cell->vertex_index(i)].insert(cell);
3342 
3343  // Take care of hanging nodes
3344  cell = triangulation.begin_active();
3345  for (; cell != endc; ++cell)
3346  {
3347  for (unsigned int i : cell->face_indices())
3348  {
3349  if ((cell->at_boundary(i) == false) &&
3350  (cell->neighbor(i)->is_active()))
3351  {
3353  adjacent_cell = cell->neighbor(i);
3354  for (unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
3355  vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert(
3356  adjacent_cell);
3357  }
3358  }
3359 
3360  // in 3d also loop over the edges
3361  if (dim == 3)
3362  {
3363  for (unsigned int i = 0; i < cell->n_lines(); ++i)
3364  if (cell->line(i)->has_children())
3365  // the only place where this vertex could have been
3366  // hiding is on the mid-edge point of the edge we
3367  // are looking at
3368  vertex_to_cell_map[cell->line(i)->child(0)->vertex_index(1)]
3369  .insert(cell);
3370  }
3371  }
3372 
3373  return vertex_to_cell_map;
3374  }
3375 
3376 
3377 
3378  template <int dim, int spacedim>
3379  std::map<unsigned int, types::global_vertex_index>
3382  {
3383  std::map<unsigned int, types::global_vertex_index>
3384  local_to_global_vertex_index;
3385 
3386 #ifndef DEAL_II_WITH_MPI
3387 
3388  // without MPI, this function doesn't make sense because on cannot
3389  // use parallel::distributed::Triangulation in any meaningful
3390  // way
3391  (void)triangulation;
3392  Assert(false,
3393  ExcMessage("This function does not make any sense "
3394  "for parallel::distributed::Triangulation "
3395  "objects if you do not have MPI enabled."));
3396 
3397 #else
3398 
3399  using active_cell_iterator =
3401  const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
3403 
3404  // Create a local index for the locally "owned" vertices
3405  types::global_vertex_index next_index = 0;
3406  unsigned int max_cellid_size = 0;
3407  std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
3408  vertices_added;
3409  std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
3410  std::map<types::subdomain_id,
3411  std::vector<std::tuple<types::global_vertex_index,
3413  std::string>>>
3414  vertices_to_send;
3415  active_cell_iterator cell = triangulation.begin_active(),
3416  endc = triangulation.end();
3417  std::set<active_cell_iterator> missing_vert_cells;
3418  std::set<unsigned int> used_vertex_index;
3419  for (; cell != endc; ++cell)
3420  {
3421  if (cell->is_locally_owned())
3422  {
3423  for (const unsigned int i : cell->vertex_indices())
3424  {
3425  types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
3426  typename std::set<active_cell_iterator>::iterator
3427  adjacent_cell = vertex_to_cell[cell->vertex_index(i)].begin(),
3428  end_adj_cell = vertex_to_cell[cell->vertex_index(i)].end();
3429  for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
3430  lowest_subdomain_id =
3431  std::min(lowest_subdomain_id,
3432  (*adjacent_cell)->subdomain_id());
3433 
3434  // See if I "own" this vertex
3435  if (lowest_subdomain_id == cell->subdomain_id())
3436  {
3437  // Check that the vertex we are working on a vertex that has
3438  // not be dealt with yet
3439  if (used_vertex_index.find(cell->vertex_index(i)) ==
3440  used_vertex_index.end())
3441  {
3442  // Set the local index
3443  local_to_global_vertex_index[cell->vertex_index(i)] =
3444  next_index++;
3445 
3446  // Store the information that will be sent to the
3447  // adjacent cells on other subdomains
3448  adjacent_cell =
3449  vertex_to_cell[cell->vertex_index(i)].begin();
3450  for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
3451  if ((*adjacent_cell)->subdomain_id() !=
3452  cell->subdomain_id())
3453  {
3454  std::pair<types::subdomain_id,
3456  tmp((*adjacent_cell)->subdomain_id(),
3457  cell->vertex_index(i));
3458  if (vertices_added.find(tmp) ==
3459  vertices_added.end())
3460  {
3461  vertices_to_send[(*adjacent_cell)
3462  ->subdomain_id()]
3463  .emplace_back(i,
3464  cell->vertex_index(i),
3465  cell->id().to_string());
3466  if (cell->id().to_string().size() >
3467  max_cellid_size)
3468  max_cellid_size =
3469  cell->id().to_string().size();
3470  vertices_added.insert(tmp);
3471  }
3472  }
3473  used_vertex_index.insert(cell->vertex_index(i));
3474  }
3475  }
3476  else
3477  {
3478  // We don't own the vertex so we will receive its global
3479  // index
3480  vertices_to_recv[lowest_subdomain_id].insert(
3481  cell->vertex_index(i));
3482  missing_vert_cells.insert(cell);
3483  }
3484  }
3485  }
3486 
3487  // Some hanging nodes are vertices of ghost cells. They need to be
3488  // received.
3489  if (cell->is_ghost())
3490  {
3491  for (unsigned int i : cell->face_indices())
3492  {
3493  if (cell->at_boundary(i) == false)
3494  {
3495  if (cell->neighbor(i)->is_active())
3496  {
3497  typename Triangulation<dim,
3498  spacedim>::active_cell_iterator
3499  adjacent_cell = cell->neighbor(i);
3500  if ((adjacent_cell->is_locally_owned()))
3501  {
3502  types::subdomain_id adj_subdomain_id =
3503  adjacent_cell->subdomain_id();
3504  if (cell->subdomain_id() < adj_subdomain_id)
3505  for (unsigned int j = 0;
3506  j < cell->face(i)->n_vertices();
3507  ++j)
3508  {
3509  vertices_to_recv[cell->subdomain_id()].insert(
3510  cell->face(i)->vertex_index(j));
3511  missing_vert_cells.insert(cell);
3512  }
3513  }
3514  }
3515  }
3516  }
3517  }
3518  }
3519 
3520  // Get the size of the largest CellID string
3521  max_cellid_size =
3522  Utilities::MPI::max(max_cellid_size, triangulation.get_communicator());
3523 
3524  // Make indices global by getting the number of vertices owned by each
3525  // processors and shifting the indices accordingly
3527  int ierr = MPI_Exscan(&next_index,
3528  &shift,
3529  1,
3531  MPI_SUM,
3532  triangulation.get_communicator());
3533  AssertThrowMPI(ierr);
3534 
3535  std::map<unsigned int, types::global_vertex_index>::iterator
3536  global_index_it = local_to_global_vertex_index.begin(),
3537  global_index_end = local_to_global_vertex_index.end();
3538  for (; global_index_it != global_index_end; ++global_index_it)
3539  global_index_it->second += shift;
3540 
3541 
3542  const int mpi_tag = Utilities::MPI::internal::Tags::
3544  const int mpi_tag2 = Utilities::MPI::internal::Tags::
3546 
3547 
3548  // In a first message, send the global ID of the vertices and the local
3549  // positions in the cells. In a second messages, send the cell ID as a
3550  // resize string. This is done in two messages so that types are not mixed
3551 
3552  // Send the first message
3553  std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
3554  vertices_to_send.size());
3555  std::vector<MPI_Request> first_requests(vertices_to_send.size());
3556  typename std::map<types::subdomain_id,
3557  std::vector<std::tuple<types::global_vertex_index,
3559  std::string>>>::iterator
3560  vert_to_send_it = vertices_to_send.begin(),
3561  vert_to_send_end = vertices_to_send.end();
3562  for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3563  ++vert_to_send_it, ++i)
3564  {
3565  int destination = vert_to_send_it->first;
3566  const unsigned int n_vertices = vert_to_send_it->second.size();
3567  const int buffer_size = 2 * n_vertices;
3568  vertices_send_buffers[i].resize(buffer_size);
3569 
3570  // fill the buffer
3571  for (unsigned int j = 0; j < n_vertices; ++j)
3572  {
3573  vertices_send_buffers[i][2 * j] =
3574  std::get<0>(vert_to_send_it->second[j]);
3575  vertices_send_buffers[i][2 * j + 1] =
3576  local_to_global_vertex_index[std::get<1>(
3577  vert_to_send_it->second[j])];
3578  }
3579 
3580  // Send the message
3581  ierr = MPI_Isend(vertices_send_buffers[i].data(),
3582  buffer_size,
3584  destination,
3585  mpi_tag,
3586  triangulation.get_communicator(),
3587  &first_requests[i]);
3588  AssertThrowMPI(ierr);
3589  }
3590 
3591  // Receive the first message
3592  std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
3593  vertices_to_recv.size());
3594  typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
3595  vert_to_recv_it = vertices_to_recv.begin(),
3596  vert_to_recv_end = vertices_to_recv.end();
3597  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3598  ++vert_to_recv_it, ++i)
3599  {
3600  int source = vert_to_recv_it->first;
3601  const unsigned int n_vertices = vert_to_recv_it->second.size();
3602  const int buffer_size = 2 * n_vertices;
3603  vertices_recv_buffers[i].resize(buffer_size);
3604 
3605  // Receive the message
3606  ierr = MPI_Recv(vertices_recv_buffers[i].data(),
3607  buffer_size,
3609  source,
3610  mpi_tag,
3611  triangulation.get_communicator(),
3612  MPI_STATUS_IGNORE);
3613  AssertThrowMPI(ierr);
3614  }
3615 
3616 
3617  // Send second message
3618  std::vector<std::vector<char>> cellids_send_buffers(
3619  vertices_to_send.size());
3620  std::vector<MPI_Request> second_requests(vertices_to_send.size());
3621  vert_to_send_it = vertices_to_send.begin();
3622  for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3623  ++vert_to_send_it, ++i)
3624  {
3625  int destination = vert_to_send_it->first;
3626  const unsigned int n_vertices = vert_to_send_it->second.size();
3627  const int buffer_size = max_cellid_size * n_vertices;
3628  cellids_send_buffers[i].resize(buffer_size);
3629 
3630  // fill the buffer
3631  unsigned int pos = 0;
3632  for (unsigned int j = 0; j < n_vertices; ++j)
3633  {
3634  std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
3635  for (unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
3636  {
3637  if (k < cell_id.size())
3638  cellids_send_buffers[i][pos] = cell_id[k];
3639  // if necessary fill up the reserved part of the buffer with an
3640  // invalid value
3641  else
3642  cellids_send_buffers[i][pos] = '-';
3643  }
3644  }
3645 
3646  // Send the message
3647  ierr = MPI_Isend(cellids_send_buffers[i].data(),
3648  buffer_size,
3649  MPI_CHAR,
3650  destination,
3651  mpi_tag2,
3652  triangulation.get_communicator(),
3653  &second_requests[i]);
3654  AssertThrowMPI(ierr);
3655  }
3656 
3657  // Receive the second message
3658  std::vector<std::vector<char>> cellids_recv_buffers(
3659  vertices_to_recv.size());
3660  vert_to_recv_it = vertices_to_recv.begin();
3661  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3662  ++vert_to_recv_it, ++i)
3663  {
3664  int source = vert_to_recv_it->first;
3665  const unsigned int n_vertices = vert_to_recv_it->second.size();
3666  const int buffer_size = max_cellid_size * n_vertices;
3667  cellids_recv_buffers[i].resize(buffer_size);
3668 
3669  // Receive the message
3670  ierr = MPI_Recv(cellids_recv_buffers[i].data(),
3671  buffer_size,
3672  MPI_CHAR,
3673  source,
3674  mpi_tag2,
3675  triangulation.get_communicator(),
3676  MPI_STATUS_IGNORE);
3677  AssertThrowMPI(ierr);
3678  }
3679 
3680 
3681  // Match the data received with the required vertices
3682  vert_to_recv_it = vertices_to_recv.begin();
3683  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3684  ++i, ++vert_to_recv_it)
3685  {
3686  for (unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
3687  {
3688  const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
3689  const types::global_vertex_index global_id_recv =
3690  vertices_recv_buffers[i][2 * j + 1];
3691  const std::string cellid_recv(
3692  &cellids_recv_buffers[i][max_cellid_size * j],
3693  &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
3694  bool found = false;
3695  typename std::set<active_cell_iterator>::iterator
3696  cell_set_it = missing_vert_cells.begin(),
3697  end_cell_set = missing_vert_cells.end();
3698  for (; (found == false) && (cell_set_it != end_cell_set);
3699  ++cell_set_it)
3700  {
3701  typename std::set<active_cell_iterator>::iterator
3702  candidate_cell =
3703  vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
3704  end_cell =
3705  vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
3706  for (; candidate_cell != end_cell; ++candidate_cell)
3707  {
3708  std::string current_cellid =
3709  (*candidate_cell)->id().to_string();
3710  current_cellid.resize(max_cellid_size, '-');
3711  if (current_cellid.compare(cellid_recv) == 0)
3712  {
3713  local_to_global_vertex_index
3714  [(*candidate_cell)->vertex_index(local_pos_recv)] =
3715  global_id_recv;
3716  found = true;
3717 
3718  break;
3719  }
3720  }
3721  }
3722  }
3723  }
3724 #endif
3725 
3726  return local_to_global_vertex_index;
3727  }
3728 
3729 
3730 
3731  template <int dim, int spacedim>
3732  void
3735  DynamicSparsityPattern & cell_connectivity)
3736  {
3737  cell_connectivity.reinit(triangulation.n_active_cells(),
3738  triangulation.n_active_cells());
3739 
3740  // loop over all cells and their neighbors to build the sparsity
3741  // pattern. note that it's a bit hard to enter all the connections when a
3742  // neighbor has children since we would need to find out which of its
3743  // children is adjacent to the current cell. this problem can be omitted
3744  // if we only do something if the neighbor has no children -- in that case
3745  // it is either on the same or a coarser level than we are. in return, we
3746  // have to add entries in both directions for both cells
3747  for (const auto &cell : triangulation.active_cell_iterators())
3748  {
3749  const unsigned int index = cell->active_cell_index();
3750  cell_connectivity.add(index, index);
3751  for (auto f : cell->face_indices())
3752  if ((cell->at_boundary(f) == false) &&
3753  (cell->neighbor(f)->has_children() == false))
3754  {
3755  const unsigned int other_index =
3756  cell->neighbor(f)->active_cell_index();
3757  cell_connectivity.add(index, other_index);
3758  cell_connectivity.add(other_index, index);
3759  }
3760  }
3761  }
3762 
3763 
3764 
3765  template <int dim, int spacedim>
3766  void
3769  DynamicSparsityPattern & cell_connectivity)
3770  {
3771  std::vector<std::vector<unsigned int>> vertex_to_cell(
3772  triangulation.n_vertices());
3773  for (const auto &cell : triangulation.active_cell_iterators())
3774  {
3775  for (const unsigned int v : cell->vertex_indices())
3776  vertex_to_cell[cell->vertex_index(v)].push_back(
3777  cell->active_cell_index());
3778  }
3779 
3780  cell_connectivity.reinit(triangulation.n_active_cells(),
3781  triangulation.n_active_cells());
3782  for (const auto &cell : triangulation.active_cell_iterators())
3783  {
3784  for (const unsigned int v : cell->vertex_indices())
3785  for (unsigned int n = 0;
3786  n < vertex_to_cell[cell->vertex_index(v)].size();
3787  ++n)
3788  cell_connectivity.add(cell->active_cell_index(),
3789  vertex_to_cell[cell->vertex_index(v)][n]);
3790  }
3791  }
3792 
3793 
3794  template <int dim, int spacedim>
3795  void
3798  const unsigned int level,
3799  DynamicSparsityPattern & cell_connectivity)
3800  {
3801  std::vector<std::vector<unsigned int>> vertex_to_cell(
3802  triangulation.n_vertices());
3803  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3804  triangulation.begin(level);
3805  cell != triangulation.end(level);
3806  ++cell)
3807  {
3808  for (const unsigned int v : cell->vertex_indices())
3809  vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
3810  }
3811 
3812  cell_connectivity.reinit(triangulation.n_cells(level),
3813  triangulation.n_cells(level));
3814  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3815  triangulation.begin(level);
3816  cell != triangulation.end(level);
3817  ++cell)
3818  {
3819  for (const unsigned int v : cell->vertex_indices())
3820  for (unsigned int n = 0;
3821  n < vertex_to_cell[cell->vertex_index(v)].size();
3822  ++n)
3823  cell_connectivity.add(cell->index(),
3824  vertex_to_cell[cell->vertex_index(v)][n]);
3825  }
3826  }
3827 
3828 
3829 
3830  template <int dim, int spacedim>
3831  void
3832  partition_triangulation(const unsigned int n_partitions,
3834  const SparsityTools::Partitioner partitioner)
3835  {
3837  &triangulation) == nullptr),
3838  ExcMessage("Objects of type parallel::distributed::Triangulation "
3839  "are already partitioned implicitly and can not be "
3840  "partitioned again explicitly."));
3841 
3842  std::vector<unsigned int> cell_weights;
3843 
3844  // Get cell weighting if a signal has been attached to the triangulation
3845  if (!triangulation.signals.weight.empty())
3846  {
3847  cell_weights.resize(triangulation.n_active_cells(), 0U);
3848 
3849  // In a first step, obtain the weights of the locally owned
3850  // cells. For all others, the weight remains at the zero the
3851  // vector was initialized with above.
3852  for (const auto &cell : triangulation.active_cell_iterators())
3853  if (cell->is_locally_owned())
3854  cell_weights[cell->active_cell_index()] =
3855  triangulation.signals.weight(
3857 
3858  // If this is a parallel triangulation, we then need to also
3859  // get the weights for all other cells. We have asserted above
3860  // that this function can't be used for
3861  // parallel::distributed::Triangulation objects, so the only
3862  // ones we have to worry about here are
3863  // parallel::shared::Triangulation
3864  if (const auto shared_tria =
3866  &triangulation))
3867  Utilities::MPI::sum(cell_weights,
3868  shared_tria->get_communicator(),
3869  cell_weights);
3870 
3871  // verify that the global sum of weights is larger than 0
3872  Assert(std::accumulate(cell_weights.begin(),
3873  cell_weights.end(),
3874  std::uint64_t(0)) > 0,
3875  ExcMessage("The global sum of weights over all active cells "
3876  "is zero. Please verify how you generate weights."));
3877  }
3878 
3879  // Call the other more general function
3880  partition_triangulation(n_partitions,
3881  cell_weights,
3882  triangulation,
3883  partitioner);
3884  }
3885 
3886 
3887 
3888  template <int dim, int spacedim>
3889  void
3890  partition_triangulation(const unsigned int n_partitions,
3891  const std::vector<unsigned int> &cell_weights,
3893  const SparsityTools::Partitioner partitioner)
3894  {
3896  &triangulation) == nullptr),
3897  ExcMessage("Objects of type parallel::distributed::Triangulation "
3898  "are already partitioned implicitly and can not be "
3899  "partitioned again explicitly."));
3900  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
3901 
3902  // check for an easy return
3903  if (n_partitions == 1)
3904  {
3905  for (const auto &cell : triangulation.active_cell_iterators())
3906  cell->set_subdomain_id(0);
3907  return;
3908  }
3909 
3910  // we decompose the domain by first
3911  // generating the connection graph of all
3912  // cells with their neighbors, and then
3913  // passing this graph off to METIS.
3914  // finally defer to the other function for
3915  // partitioning and assigning subdomain ids
3916  DynamicSparsityPattern cell_connectivity;
3917  get_face_connectivity_of_cells(triangulation, cell_connectivity);
3918 
3919  SparsityPattern sp_cell_connectivity;
3920  sp_cell_connectivity.copy_from(cell_connectivity);
3921  partition_triangulation(n_partitions,
3922  cell_weights,
3923  sp_cell_connectivity,
3924  triangulation,
3925  partitioner);
3926  }
3927 
3928 
3929 
3930  template <int dim, int spacedim>
3931  void
3932  partition_triangulation(const unsigned int n_partitions,
3933  const SparsityPattern & cell_connection_graph,
3935  const SparsityTools::Partitioner partitioner)
3936  {
3938  &triangulation) == nullptr),
3939  ExcMessage("Objects of type parallel::distributed::Triangulation "
3940  "are already partitioned implicitly and can not be "
3941  "partitioned again explicitly."));
3942 
3943  std::vector<unsigned int> cell_weights;
3944 
3945  // Get cell weighting if a signal has been attached to the triangulation
3946  if (!triangulation.signals.weight.empty())
3947  {
3948  cell_weights.resize(triangulation.n_active_cells(), 0U);
3949 
3950  // In a first step, obtain the weights of the locally owned
3951  // cells. For all others, the weight remains at the zero the
3952  // vector was initialized with above.
3953  for (const auto &cell : triangulation.active_cell_iterators() |
3955  cell_weights[cell->active_cell_index()] =
3956  triangulation.signals.weight(
3958 
3959  // If this is a parallel triangulation, we then need to also
3960  // get the weights for all other cells. We have asserted above
3961  // that this function can't be used for
3962  // parallel::distribute::Triangulation objects, so the only
3963  // ones we have to worry about here are
3964  // parallel::shared::Triangulation
3965  if (const auto shared_tria =
3967  &triangulation))
3968  Utilities::MPI::sum(cell_weights,
3969  shared_tria->get_communicator(),
3970  cell_weights);
3971 
3972  // verify that the global sum of weights is larger than 0
3973  Assert(std::accumulate(cell_weights.begin(),
3974  cell_weights.end(),
3975  std::uint64_t(0)) > 0,
3976  ExcMessage("The global sum of weights over all active cells "
3977  "is zero. Please verify how you generate weights."));
3978  }
3979 
3980  // Call the other more general function
3981  partition_triangulation(n_partitions,
3982  cell_weights,
3983  cell_connection_graph,
3984  triangulation,
3985  partitioner);
3986  }
3987 
3988 
3989 
3990  template <int dim, int spacedim>
3991  void
3992  partition_triangulation(const unsigned int n_partitions,
3993  const std::vector<unsigned int> &cell_weights,
3994  const SparsityPattern & cell_connection_graph,
3996  const SparsityTools::Partitioner partitioner)
3997  {
3999  &triangulation) == nullptr),
4000  ExcMessage("Objects of type parallel::distributed::Triangulation "
4001  "are already partitioned implicitly and can not be "
4002  "partitioned again explicitly."));
4003  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
4004  Assert(cell_connection_graph.n_rows() == triangulation.n_active_cells(),
4005  ExcMessage("Connectivity graph has wrong size"));
4006  Assert(cell_connection_graph.n_cols() == triangulation.n_active_cells(),
4007  ExcMessage("Connectivity graph has wrong size"));
4008 
4009  // signal that partitioning is going to happen
4010  triangulation.signals.pre_partition();
4011 
4012  // check for an easy return
4013  if (n_partitions == 1)
4014  {
4015  for (const auto &cell : triangulation.active_cell_iterators())
4016  cell->set_subdomain_id(0);
4017  return;
4018  }
4019 
4020  // partition this connection graph and get
4021  // back a vector of indices, one per degree
4022  // of freedom (which is associated with a
4023  // cell)
4024  std::vector<unsigned int> partition_indices(triangulation.n_active_cells());
4025  SparsityTools::partition(cell_connection_graph,
4026  cell_weights,
4027  n_partitions,
4028  partition_indices,
4029  partitioner);
4030 
4031  // finally loop over all cells and set the subdomain ids
4032  for (const auto &cell : triangulation.active_cell_iterators())
4033  cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
4034  }
4035 
4036 
4037  namespace internal
4038  {
4042  template <class IT>
4043  void
4045  unsigned int & current_proc_idx,
4046  unsigned int & current_cell_idx,
4047  const unsigned int n_active_cells,
4048  const unsigned int n_partitions)
4049  {
4050  if (cell->is_active())
4051  {
4052  while (current_cell_idx >=
4053  std::floor(static_cast<uint_least64_t>(n_active_cells) *
4054  (current_proc_idx + 1) / n_partitions))
4055  ++current_proc_idx;
4056  cell->set_subdomain_id(current_proc_idx);
4057  ++current_cell_idx;
4058  }
4059  else
4060  {
4061  for (unsigned int n = 0; n < cell->n_children(); ++n)
4063  current_proc_idx,
4064  current_cell_idx,
4066  n_partitions);
4067  }
4068  }
4069  } // namespace internal
4070 
4071  template <int dim, int spacedim>
4072  void
4073  partition_triangulation_zorder(const unsigned int n_partitions,
4075  const bool group_siblings)
4076  {
4078  &triangulation) == nullptr),
4079  ExcMessage("Objects of type parallel::distributed::Triangulation "
4080  "are already partitioned implicitly and can not be "
4081  "partitioned again explicitly."));
4082  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
4083  Assert(triangulation.signals.weight.empty(), ExcNotImplemented());
4084 
4085  // signal that partitioning is going to happen
4086  triangulation.signals.pre_partition();
4087 
4088  // check for an easy return
4089  if (n_partitions == 1)
4090  {
4091  for (const auto &cell : triangulation.active_cell_iterators())
4092  cell->set_subdomain_id(0);
4093  return;
4094  }
4095 
4096  // Duplicate the coarse cell reordoring
4097  // as done in p4est
4098  std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
4099  std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
4100 
4101  DynamicSparsityPattern cell_connectivity;
4103  0,
4104  cell_connectivity);
4105  coarse_cell_to_p4est_tree_permutation.resize(triangulation.n_cells(0));
4106  SparsityTools::reorder_hierarchical(cell_connectivity,
4107  coarse_cell_to_p4est_tree_permutation);
4108 
4109  p4est_tree_to_coarse_cell_permutation =
4110  Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
4111 
4112  unsigned int current_proc_idx = 0;
4113  unsigned int current_cell_idx = 0;
4114  const unsigned int n_active_cells = triangulation.n_active_cells();
4115 
4116  // set subdomain id for active cell descendants
4117  // of each coarse cell in permuted order
4118  for (unsigned int idx = 0; idx < triangulation.n_cells(0); ++idx)
4119  {
4120  const unsigned int coarse_cell_idx =
4121  p4est_tree_to_coarse_cell_permutation[idx];
4122  typename Triangulation<dim, spacedim>::cell_iterator coarse_cell(
4123  &triangulation, 0, coarse_cell_idx);
4124 
4126  current_proc_idx,
4127  current_cell_idx,
4129  n_partitions);
4130  }
4131 
4132  // if all children of a cell are active (e.g. we
4133  // have a cell that is refined once and no part
4134  // is refined further), p4est places all of them
4135  // on the same processor. The new owner will be
4136  // the processor with the largest number of children
4137  // (ties are broken by picking the lower rank).
4138  // Duplicate this logic here.
4139  if (group_siblings)
4140  {
4142  cell = triangulation.begin(),
4143  endc = triangulation.end();
4144  for (; cell != endc; ++cell)
4145  {
4146  if (cell->is_active())
4147  continue;
4148  bool all_children_active = true;
4149  std::map<unsigned int, unsigned int> map_cpu_n_cells;
4150  for (unsigned int n = 0; n < cell->n_children(); ++n)
4151  if (!cell->child(n)->is_active())
4152  {
4153  all_children_active = false;
4154  break;
4155  }
4156  else
4157  ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
4158 
4159  if (!all_children_active)
4160  continue;
4161 
4162  unsigned int new_owner = cell->child(0)->subdomain_id();
4163  for (std::map<unsigned int, unsigned int>::iterator it =
4164  map_cpu_n_cells.begin();
4165  it != map_cpu_n_cells.end();
4166  ++it)
4167  if (it->second > map_cpu_n_cells[new_owner])
4168  new_owner = it->first;
4169 
4170  for (unsigned int n = 0; n < cell->n_children(); ++n)
4171  cell->child(n)->set_subdomain_id(new_owner);
4172  }
4173  }
4174  }
4175 
4176 
4177  template <int dim, int spacedim>
4178  void
4180  {
4181  unsigned int n_levels = triangulation.n_levels();
4182  for (int lvl = n_levels - 1; lvl >= 0; --lvl)
4183  {
4184  for (const auto &cell : triangulation.cell_iterators_on_level(lvl))
4185  {
4186  if (cell->is_active())
4187  cell->set_level_subdomain_id(cell->subdomain_id());
4188  else
4189  {
4190  Assert(cell->child(0)->level_subdomain_id() !=
4192  ExcInternalError());
4193  cell->set_level_subdomain_id(
4194  cell->child(0)->level_subdomain_id());
4195  }
4196  }
4197  }
4198  }
4199 
4200  namespace internal
4201  {
4202  namespace
4203  {
4204  // Split get_subdomain_association() for p::d::T since we want to compile
4205  // it in 1D but none of the p4est stuff is available in 1D.
4206  template <int dim, int spacedim>
4207  void
4210  & triangulation,
4211  const std::vector<CellId> & cell_ids,
4212  std::vector<types::subdomain_id> &subdomain_ids)
4213  {
4214 #ifndef DEAL_II_WITH_P4EST
4215  (void)triangulation;
4216  (void)cell_ids;
4217  (void)subdomain_ids;
4218  Assert(
4219  false,
4220  ExcMessage(
4221  "You are attempting to use a functionality that is only available "
4222  "if deal.II was configured to use p4est, but cmake did not find a "
4223  "valid p4est library."));
4224 #else
4225  // for parallel distributed triangulations, we will ask the p4est oracle
4226  // about the global partitioning of active cells since this information
4227  // is stored on every process
4228  for (const auto &cell_id : cell_ids)
4229  {
4230  // find descendent from coarse quadrant
4231  typename ::internal::p4est::types<dim>::quadrant p4est_cell,
4233 
4234  ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
4235  for (const auto &child_index : cell_id.get_child_indices())
4236  {
4237  ::internal::p4est::init_quadrant_children<dim>(
4238  p4est_cell, p4est_children);
4239  p4est_cell =
4240  p4est_children[static_cast<unsigned int>(child_index)];
4241  }
4242 
4243  // find owning process, i.e., the subdomain id
4244  const int owner =
4246  const_cast<typename ::internal::p4est::types<dim>::forest
4247  *>(triangulation.get_p4est()),
4248  cell_id.get_coarse_cell_id(),
4249  &p4est_cell,
4251  triangulation.get_communicator()));
4252 
4253  Assert(owner >= 0, ExcMessage("p4est should know the owner."));
4254 
4255  subdomain_ids.push_back(owner);
4256  }
4257 #endif
4258  }
4259 
4260 
4261 
4262  template <int spacedim>
4263  void
4266  const std::vector<CellId> &,
4267  std::vector<types::subdomain_id> &)
4268  {
4269  Assert(false, ExcNotImplemented());
4270  }
4271  } // anonymous namespace
4272  } // namespace internal
4273 
4274 
4275 
4276  template <int dim, int spacedim>
4277  std::vector<types::subdomain_id>
4279  const std::vector<CellId> & cell_ids)
4280  {
4281  std::vector<types::subdomain_id> subdomain_ids;
4282  subdomain_ids.reserve(cell_ids.size());
4283 
4284  if (dynamic_cast<
4286  &triangulation) != nullptr)
4287  {
4288  Assert(false, ExcNotImplemented());
4289  }
4291  *parallel_tria = dynamic_cast<
4293  &triangulation))
4294  {
4295  internal::get_subdomain_association(*parallel_tria,
4296  cell_ids,
4297  subdomain_ids);
4298  }
4299  else if (const parallel::shared::Triangulation<dim, spacedim> *shared_tria =
4301  *>(&triangulation))
4302  {
4303  // for parallel shared triangulations, we need to access true subdomain
4304  // ids which are also valid for artificial cells
4305  const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
4306  shared_tria->get_true_subdomain_ids_of_cells();
4307 
4308  for (const auto &cell_id : cell_ids)
4309  {
4310  const unsigned int active_cell_index =
4311  shared_tria->create_cell_iterator(cell_id)->active_cell_index();
4312  subdomain_ids.push_back(
4313  true_subdomain_ids_of_cells[active_cell_index]);
4314  }
4315  }
4316  else
4317  {
4318  // the most general type of triangulation is the serial one. here, all
4319  // subdomain information is directly available
4320  for (const auto &cell_id : cell_ids)
4321  {
4322  subdomain_ids.push_back(
4323  triangulation.create_cell_iterator(cell_id)->subdomain_id());
4324  }
4325  }
4326 
4327  return subdomain_ids;
4328  }
4329 
4330 
4331 
4332  template <int dim, int spacedim>
4333  void
4335  std::vector<types::subdomain_id> & subdomain)
4336  {
4337  Assert(subdomain.size() == triangulation.n_active_cells(),
4338  ExcDimensionMismatch(subdomain.size(),
4339  triangulation.n_active_cells()));
4340  for (const auto &cell : triangulation.active_cell_iterators())
4341  subdomain[cell->active_cell_index()] = cell->subdomain_id();
4342  }
4343 
4344 
4345 
4346  template <int dim, int spacedim>
4347  unsigned int
4350  const types::subdomain_id subdomain)
4351  {
4352  unsigned int count = 0;
4353  for (const auto &cell : triangulation.active_cell_iterators())
4354  if (cell->subdomain_id() == subdomain)
4355  ++count;
4356 
4357  return count;
4358  }
4359 
4360 
4361 
4362  template <int dim, int spacedim>
4363  std::vector<bool>
4365  {
4366  // start with all vertices
4367  std::vector<bool> locally_owned_vertices =
4368  triangulation.get_used_vertices();
4369 
4370  // if the triangulation is distributed, eliminate those that
4371  // are owned by other processors -- either because the vertex is
4372  // on an artificial cell, or because it is on a ghost cell with
4373  // a smaller subdomain
4374  if (const auto *tr = dynamic_cast<
4376  &triangulation))
4377  for (const auto &cell : triangulation.active_cell_iterators())
4378  if (cell->is_artificial() ||
4379  (cell->is_ghost() &&
4380  (cell->subdomain_id() < tr->locally_owned_subdomain())))
4381  for (const unsigned int v : cell->vertex_indices())
4382  locally_owned_vertices[cell->vertex_index(v)] = false;
4383 
4384  return locally_owned_vertices;
4385  }
4386 
4387 
4388 
4389  template <int dim, int spacedim>
4390  double
4392  const Mapping<dim, spacedim> & mapping)
4393  {
4394  double min_diameter = std::numeric_limits<double>::max();
4395  for (const auto &cell : triangulation.active_cell_iterators())
4396  if (!cell->is_artificial())
4397  min_diameter = std::min(min_diameter, cell->diameter(mapping));
4398 
4399  double global_min_diameter = 0;
4400 
4401 #ifdef DEAL_II_WITH_MPI
4402  if (const parallel::TriangulationBase<dim, spacedim> *p_tria =
4403  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4404  &triangulation))
4405  global_min_diameter =
4406  Utilities::MPI::min(min_diameter, p_tria->get_communicator());
4407  else
4408 #endif
4409  global_min_diameter = min_diameter;
4410 
4411  return global_min_diameter;
4412  }
4413 
4414 
4415 
4416  template <int dim, int spacedim>
4417  double
4419  const Mapping<dim, spacedim> & mapping)
4420  {
4421  double max_diameter = 0.;
4422  for (const auto &cell : triangulation.active_cell_iterators())
4423  if (!cell->is_artificial())
4424  max_diameter = std::max(max_diameter, cell->diameter(mapping));
4425 
4426  double global_max_diameter = 0;
4427 
4428 #ifdef DEAL_II_WITH_MPI
4429  if (const parallel::TriangulationBase<dim, spacedim> *p_tria =
4430  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4431  &triangulation))
4432  global_max_diameter =
4433  Utilities::MPI::max(max_diameter, p_tria->get_communicator());
4434  else
4435 #endif
4436  global_max_diameter = max_diameter;
4437 
4438  return global_max_diameter;
4439  }
4440 
4441 
4442 
4443  namespace internal
4444  {
4445  namespace FixUpDistortedChildCells
4446  {
4447  // compute the mean square
4448  // deviation of the alternating
4449  // forms of the children of the
4450  // given object from that of
4451  // the object itself. for
4452  // objects with
4453  // structdim==spacedim, the
4454  // alternating form is the
4455  // determinant of the jacobian,
4456  // whereas for faces with
4457  // structdim==spacedim-1, the
4458  // alternating form is the
4459  // (signed and scaled) normal
4460  // vector
4461  //
4462  // this average square
4463  // deviation is computed for an
4464  // object where the center node
4465  // has been replaced by the
4466  // second argument to this
4467  // function
4468  template <typename Iterator, int spacedim>
4469  double
4470  objective_function(const Iterator & object,
4471  const Point<spacedim> &object_mid_point)
4472  {
4473  const unsigned int structdim =
4474  Iterator::AccessorType::structure_dimension;
4475  Assert(spacedim == Iterator::AccessorType::dimension,
4476  ExcInternalError());
4477 
4478  // everything below is wrong
4479  // if not for the following
4480  // condition
4481  Assert(object->refinement_case() ==
4483  ExcNotImplemented());
4484  // first calculate the
4485  // average alternating form
4486  // for the parent cell/face
4489  Tensor<spacedim - structdim, spacedim>
4490  parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
4491 
4492  for (const unsigned int i : object->vertex_indices())
4493  parent_vertices[i] = object->vertex(i);
4494 
4496  parent_vertices, parent_alternating_forms);
4497 
4498  const Tensor<spacedim - structdim, spacedim>
4499  average_parent_alternating_form =
4500  std::accumulate(parent_alternating_forms,
4501  parent_alternating_forms +
4504 
4505  // now do the same
4506  // computation for the
4507  // children where we use the
4508  // given location for the
4509  // object mid point instead of
4510  // the one the triangulation
4511  // currently reports
4515  Tensor<spacedim - structdim, spacedim> child_alternating_forms
4518 
4519  for (unsigned int c = 0; c < object->n_children(); ++c)
4520  for (const unsigned int i : object->child(c)->vertex_indices())
4521  child_vertices[c][i] = object->child(c)->vertex(i);
4522 
4523  // replace mid-object
4524  // vertex. note that for
4525  // child i, the mid-object
4526  // vertex happens to have the
4527  // number
4528  // max_children_per_cell-i
4529  for (unsigned int c = 0; c < object->n_children(); ++c)
4530  child_vertices[c][GeometryInfo<structdim>::max_children_per_cell - c -
4531  1] = object_mid_point;
4532 
4533  for (unsigned int c = 0; c < object->n_children(); ++c)
4535  child_vertices[c], child_alternating_forms[c]);
4536 
4537  // on a uniformly refined
4538  // hypercube object, the child
4539  // alternating forms should
4540  // all be smaller by a factor
4541  // of 2^structdim than the
4542  // ones of the parent. as a
4543  // consequence, we'll use the
4544  // squared deviation from
4545  // this ideal value as an
4546  // objective function
4547  double objective = 0;
4548  for (unsigned int c = 0; c < object->n_children(); ++c)
4549  for (const unsigned int i : object->child(c)->vertex_indices())
4550  objective +=
4551  (child_alternating_forms[c][i] -
4552  average_parent_alternating_form / std::pow(2., 1. * structdim))
4553  .norm_square();
4554 
4555  return objective;
4556  }
4557 
4558 
4564  template <typename Iterator>
4566  get_face_midpoint(const Iterator & object,
4567  const unsigned int f,
4568  std::integral_constant<int, 1>)
4569  {
4570  return object->vertex(f);
4571  }
4572 
4573 
4574 
4580  template <typename Iterator>
4582  get_face_midpoint(const Iterator & object,
4583  const unsigned int f,
4584  std::integral_constant<int, 2>)
4585  {
4586  return object->line(f)->center();
4587  }
4588 
4589 
4590 
4596  template <typename Iterator>
4598  get_face_midpoint(const Iterator & object,
4599  const unsigned int f,
4600  std::integral_constant<int, 3>)
4601  {
4602  return object->face(f)->center();
4603  }
4604 
4605 
4606 
4629  template <typename Iterator>
4630  double
4631  minimal_diameter(const Iterator &object)
4632  {
4633  const unsigned int structdim =
4634  Iterator::AccessorType::structure_dimension;
4635 
4636  double diameter = object->diameter();
4637  for (const unsigned int f : object->face_indices())
4638  for (unsigned int e = f + 1; e < object->n_faces(); ++e)
4639  diameter = std::min(
4640  diameter,
4641  get_face_midpoint(object,
4642  f,
4643  std::integral_constant<int, structdim>())
4644  .distance(get_face_midpoint(
4645  object, e, std::integral_constant<int, structdim>())));
4646 
4647  return diameter;
4648  }
4649 
4650 
4651 
4656  template <typename Iterator>
4657  bool
4658  fix_up_object(const Iterator &object)
4659  {
4660  const unsigned int structdim =
4661  Iterator::AccessorType::structure_dimension;
4662  const unsigned int spacedim = Iterator::AccessorType::space_dimension;
4663 
4664  // right now we can only deal with cells that have been refined
4665  // isotropically because that is the only case where we have a cell
4666  // mid-point that can be moved around without having to consider
4667  // boundary information
4668  Assert(object->has_children(), ExcInternalError());
4669  Assert(object->refinement_case() ==
4671  ExcNotImplemented());
4672 
4673  // get the current location of the object mid-vertex:
4674  Point<spacedim> object_mid_point = object->child(0)->vertex(
4676 
4677  // now do a few steepest descent steps to reduce the objective
4678  // function. compute the diameter in the helper function above
4679  unsigned int iteration = 0;
4680  const double diameter = minimal_diameter(object);
4681 
4682  // current value of objective function and initial delta
4683  double current_value = objective_function(object, object_mid_point);
4684  double initial_delta = 0;
4685 
4686  do
4687  {
4688  // choose a step length that is initially 1/4 of the child
4689  // objects' diameter, and a sequence whose sum does not converge
4690  // (to avoid premature termination of the iteration)
4691  const double step_length = diameter / 4 / (iteration + 1);
4692 
4693  // compute the objective function's derivative using a two-sided
4694  // difference formula with eps=step_length/10
4695  Tensor<1, spacedim> gradient;
4696  for (unsigned int d = 0; d < spacedim; ++d)
4697  {
4698  const double eps = step_length / 10;
4699 
4701  h[d] = eps / 2;
4702 
4703  gradient[d] =
4705  object, project_to_object(object, object_mid_point + h)) -
4707  object, project_to_object(object, object_mid_point - h))) /
4708  eps;
4709  }
4710 
4711  // there is nowhere to go
4712  if (gradient.norm() == 0)
4713  break;
4714 
4715  // We need to go in direction -gradient. the optimal value of the
4716  // objective function is zero, so assuming that the model is
4717  // quadratic we would have to go -2*val/||gradient|| in this
4718  // direction, make sure we go at most step_length into this
4719  // direction
4720  object_mid_point -=
4721  std::min(2 * current_value / (gradient * gradient),
4722  step_length / gradient.norm()) *
4723  gradient;
4724  object_mid_point = project_to_object(object, object_mid_point);
4725 
4726  // compute current value of the objective function
4727  const double previous_value = current_value;
4728  current_value = objective_function(object, object_mid_point);
4729 
4730  if (iteration == 0)
4731  initial_delta = (previous_value - current_value);
4732 
4733  // stop if we aren't moving much any more
4734  if ((iteration >= 1) &&
4735  ((previous_value - current_value < 0) ||
4736  (std::fabs(previous_value - current_value) <
4737  0.001 * initial_delta)))
4738  break;
4739 
4740  ++iteration;
4741  }
4742  while (iteration < 20);
4743 
4744  // verify that the new
4745  // location is indeed better
4746  // than the one before. check
4747  // this by comparing whether
4748  // the minimum value of the
4749  // products of parent and
4750  // child alternating forms is
4751  // positive. for cells this
4752  // means that the
4753  // determinants have the same
4754  // sign, for faces that the
4755  // face normals of parent and
4756  // children point in the same
4757  // general direction
4758  double old_min_product, new_min_product;
4759 
4762  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
4763  parent_vertices[i] = object->vertex(i);
4764 
4765  Tensor<spacedim - structdim, spacedim>
4766  parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
4768  parent_vertices, parent_alternating_forms);
4769 
4773 
4774  for (unsigned int c = 0; c < object->n_children(); ++c)
4775  for (const unsigned int i : object->child(c)->vertex_indices())
4776  child_vertices[c][i] = object->child(c)->vertex(i);
4777 
4778  Tensor<spacedim - structdim, spacedim> child_alternating_forms
4781 
4782  for (unsigned int c = 0; c < object->n_children(); ++c)
4784  child_vertices[c], child_alternating_forms[c]);
4785 
4786  old_min_product =
4787  child_alternating_forms[0][0] * parent_alternating_forms[0];
4788  for (unsigned int c = 0; c < object->n_children(); ++c)
4789  for (const unsigned int i : object->child(c)->vertex_indices())
4790  for (const unsigned int j : object->vertex_indices())
4791  old_min_product = std::min<double>(old_min_product,
4792  child_alternating_forms[c][i] *
4793  parent_alternating_forms[j]);
4794 
4795  // for the new minimum value,
4796  // replace mid-object
4797  // vertex. note that for child
4798  // i, the mid-object vertex
4799  // happens to have the number
4800  // max_children_per_cell-i
4801  for (unsigned int c = 0; c < object->n_children(); ++c)
4802  child_vertices[c][GeometryInfo<structdim>::max_children_per_cell - c -
4803  1] = object_mid_point;
4804 
4805  for (unsigned int c = 0; c < object->n_children(); ++c)
4807  child_vertices[c], child_alternating_forms[c]);
4808 
4809  new_min_product =
4810  child_alternating_forms[0][0] * parent_alternating_forms[0];
4811  for (unsigned int c = 0; c < object->n_children(); ++c)
4812  for (const unsigned int i : object->child(c)->vertex_indices())
4813  for (const unsigned int j : object->vertex_indices())
4814  new_min_product = std::min<double>(new_min_product,
4815  child_alternating_forms[c][i] *
4816  parent_alternating_forms[j]);
4817 
4818  // if new minimum value is
4819  // better than before, then set the
4820  // new mid point. otherwise
4821  // return this object as one of
4822  // those that can't apparently
4823  // be fixed
4824  if (new_min_product >= old_min_product)
4825  object->child(0)->vertex(
4827  object_mid_point;
4828 
4829  // return whether after this
4830  // operation we have an object that
4831  // is well oriented
4832  return (std::max(new_min_product, old_min_product) > 0);
4833  }
4834 
4835 
4836 
4837  // possibly fix up the faces of a cell by moving around its mid-points
4838  template <int dim, int spacedim>
4839  void
4841  const typename ::Triangulation<dim, spacedim>::cell_iterator
4842  &cell,
4843  std::integral_constant<int, dim>,
4844  std::integral_constant<int, spacedim>)
4845  {
4846  // see if we first can fix up some of the faces of this object. We can
4847  // mess with faces if and only if the neighboring cell is not even
4848  // more refined than we are (since in that case the sub-faces have
4849  // themselves children that we can't move around any more). however,
4850  // the latter case shouldn't happen anyway: if the current face is
4851  // distorted but the neighbor is even more refined, then the face had
4852  // been deformed before already, and had been ignored at the time; we
4853  // should then also be able to ignore it this time as well
4854  for (auto f : cell->face_indices())
4855  {
4856  Assert(cell->face(f)->has_children(), ExcInternalError());
4857  Assert(cell->face(f)->refinement_case() ==
4859  ExcInternalError());
4860 
4861  bool subface_is_more_refined = false;
4862  for (unsigned int g = 0;
4863  g < GeometryInfo<dim>::max_children_per_face;
4864  ++g)
4865  if (cell->face(f)->child(g)->has_children())
4866  {
4867  subface_is_more_refined = true;
4868  break;
4869  }
4870 
4871  if (subface_is_more_refined == true)
4872  continue;
4873 
4874  // we finally know that we can do something about this face
4875  fix_up_object(cell->face(f));
4876  }
4877  }
4878  } /* namespace FixUpDistortedChildCells */
4879  } /* namespace internal */
4880 
4881 
4882  template <int dim, int spacedim>
4886  &distorted_cells,
4887  Triangulation<dim, spacedim> & /*triangulation*/)
4888  {
4889  static_assert(
4890  dim != 1 && spacedim != 1,
4891  "This function is only valid when dim != 1 or spacedim != 1.");
4892  typename Triangulation<dim, spacedim>::DistortedCellList unfixable_subset;
4893 
4894  // loop over all cells that we have to fix up
4895  for (typename std::list<
4896  typename Triangulation<dim, spacedim>::cell_iterator>::const_iterator
4897  cell_ptr = distorted_cells.distorted_cells.begin();
4898  cell_ptr != distorted_cells.distorted_cells.end();
4899  ++cell_ptr)
4900  {
4901  const typename Triangulation<dim, spacedim>::cell_iterator cell =
4902  *cell_ptr;
4903 
4904  Assert(!cell->is_active(),
4905  ExcMessage(
4906  "This function is only valid for a list of cells that "
4907  "have children (i.e., no cell in the list may be active)."));
4908 
4910  cell,
4911  std::integral_constant<int, dim>(),
4912  std::integral_constant<int, spacedim>());
4913 
4914  // If possible, fix up the object.
4916  unfixable_subset.distorted_cells.push_back(cell);
4917  }
4918 
4919  return unfixable_subset;
4920  }
4921 
4922 
4923 
4924  template <int dim, int spacedim>
4925  void
4927  const bool reset_boundary_ids)
4928  {
4929  const auto src_boundary_ids = tria.get_boundary_ids();
4930  std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
4931  auto m_it = dst_manifold_ids.begin();
4932  for (const auto b : src_boundary_ids)
4933  {
4934  *m_it = static_cast<types::manifold_id>(b);
4935  ++m_it;
4936  }
4937  const std::vector<types::boundary_id> reset_boundary_id =
4938  reset_boundary_ids ?
4939  std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
4940  src_boundary_ids;
4941  map_boundary_to_manifold_ids(src_boundary_ids,
4942  dst_manifold_ids,
4943  tria,
4944  reset_boundary_id);
4945  }
4946 
4947 
4948 
4949  template <int dim, int spacedim>
4950  void
4952  const std::vector<types::boundary_id> &src_boundary_ids,
4953  const std::vector<types::manifold_id> &dst_manifold_ids,
4955  const std::vector<types::boundary_id> &reset_boundary_ids_)
4956  {
4957  AssertDimension(src_boundary_ids.size(), dst_manifold_ids.size());
4958  const auto reset_boundary_ids =
4959  reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
4960  AssertDimension(reset_boundary_ids.size(), src_boundary_ids.size());
4961 
4962  // in 3d, we not only have to copy boundary ids of faces, but also of edges
4963  // because we see them twice (once from each adjacent boundary face),
4964  // we cannot immediately reset their boundary ids. thus, copy first
4965  // and reset later
4966  if (dim >= 3)
4967  for (const auto &cell : tria.active_cell_iterators())
4968  for (auto f : cell->face_indices())
4969  if (cell->face(f)->at_boundary())
4970  for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
4971  {
4972  const auto bid = cell->face(f)->line(e)->boundary_id();
4973  const unsigned int ind = std::find(src_boundary_ids.begin(),
4974  src_boundary_ids.end(),
4975  bid) -
4976  src_boundary_ids.begin();
4977  if (ind < src_boundary_ids.size())
4978  cell->face(f)->line(e)->set_manifold_id(
4979  dst_manifold_ids[ind]);
4980  }
4981 
4982  // now do cells
4983  for (const auto &cell : tria.active_cell_iterators())
4984  for (auto f : cell->face_indices())
4985  if (cell->face(f)->at_boundary())
4986  {
4987  const auto bid = cell->face(f)->boundary_id();
4988  const unsigned int ind =
4989  std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
4990  src_boundary_ids.begin();
4991 
4992  if (ind < src_boundary_ids.size())
4993  {
4994  // assign the manifold id
4995  cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
4996  // then reset boundary id
4997  cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
4998  }
4999 
5000  if (dim >= 3)
5001  for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
5002  {
5003  const auto bid = cell->face(f)->line(e)->boundary_id();
5004  const unsigned int ind = std::find(src_boundary_ids.begin(),
5005  src_boundary_ids.end(),
5006  bid) -
5007  src_boundary_ids.begin();
5008  if (ind < src_boundary_ids.size())
5009  cell->face(f)->line(e)->set_boundary_id(
5010  reset_boundary_ids[ind]);
5011  }
5012  }
5013  }
5014 
5015 
5016  template <int dim, int spacedim>
5017  void
5019  const bool compute_face_ids)
5020  {
5022  cell = tria.begin_active(),
5023  endc = tria.end();
5024 
5025  for (; cell != endc; ++cell)
5026  {
5027  cell->set_manifold_id(cell->material_id());
5028  if (compute_face_ids == true)
5029  {
5030  for (auto f : cell->face_indices())
5031  {
5032  if (cell->at_boundary(f) == false)
5033  cell->face(f)->set_manifold_id(
5034  std::min(cell->material_id(),
5035  cell->neighbor(f)->material_id()));
5036  else
5037  cell->face(f)->set_manifold_id(cell->material_id());
5038  }
5039  }
5040  }
5041  }
5042 
5043 
5044  template <int dim, int spacedim>
5045  void
5048  const std::function<types::manifold_id(
5049  const std::set<types::manifold_id> &)> &disambiguation_function,
5050  bool overwrite_only_flat_manifold_ids)
5051  {
5052  // Easy case first:
5053  if (dim == 1)
5054  return;
5055  const unsigned int n_subobjects =
5056  dim == 2 ? tria.n_lines() : tria.n_lines() + tria.n_quads();
5057 
5058  // If user index is zero, then it has not been set.
5059  std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
5060  std::vector<unsigned int> backup;
5061  tria.save_user_indices(backup);
5063 
5064  unsigned next_index = 1;
5065  for (auto &cell : tria.active_cell_iterators())
5066  {
5067  if (dim > 1)
5068  for (unsigned int l = 0; l < cell->n_lines(); ++l)
5069  {
5070  if (cell->line(l)->user_index() == 0)
5071  {
5072  AssertIndexRange(next_index, n_subobjects + 1);
5073  manifold_ids[next_index].insert(cell->manifold_id());
5074  cell->line(l)->set_user_index(next_index++);
5075  }
5076  else
5077  manifold_ids[cell->line(l)->user_index()].insert(
5078  cell->manifold_id());
5079  }
5080  if (dim > 2)
5081  for (unsigned int l = 0; l < cell->n_faces(); ++l)
5082  {
5083  if (cell->quad(l)->user_index() == 0)
5084  {
5085  AssertIndexRange(next_index, n_subobjects + 1);
5086  manifold_ids[next_index].insert(cell->manifold_id());
5087  cell->quad(l)->set_user_index(next_index++);
5088  }
5089  else
5090  manifold_ids[cell->quad(l)->user_index()].insert(
5091  cell->manifold_id());
5092  }
5093  }
5094  for (auto &cell : tria.active_cell_iterators())
5095  {
5096  if (dim > 1)
5097  for (unsigned int l = 0; l < cell->n_lines(); ++l)
5098  {
5099  const auto id = cell->line(l)->user_index();
5100  // Make sure we change the manifold indicator only once
5101  if (id != 0)
5102  {
5103  if (cell->line(l)->manifold_id() ==
5105  overwrite_only_flat_manifold_ids == false)
5106  cell->line(l)->set_manifold_id(
5107  disambiguation_function(manifold_ids[id]));
5108  cell->line(l)->set_user_index(0);
5109  }
5110  }
5111  if (dim > 2)
5112  for (unsigned int l = 0; l < cell->n_faces(); ++l)
5113  {
5114  const auto id = cell->quad(l)->user_index();
5115  // Make sure we change the manifold indicator only once
5116  if (id != 0)
5117  {
5118  if (cell->quad(l)->manifold_id() ==
5120  overwrite_only_flat_manifold_ids == false)
5121  cell->quad(l)->set_manifold_id(
5122  disambiguation_function(manifold_ids[id]));
5123  cell->quad(l)->set_user_index(0);
5124  }
5125  }
5126  }
5127  tria.load_user_indices(backup);
5128  }
5129 
5130 
5131 
5132  template <int dim, int spacedim>
5133  std::pair<unsigned int, double>
5136  {
5137  double max_ratio = 1;
5138  unsigned int index = 0;
5139 
5140  for (unsigned int i = 0; i < dim; ++i)
5141  for (unsigned int j = i + 1; j < dim; ++j)
5142  {
5143  unsigned int ax = i % dim;
5144  unsigned int next_ax = j % dim;
5145 
5146  double ratio =
5147  cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
5148 
5149  if (ratio > max_ratio)
5150  {
5151  max_ratio = ratio;
5152  index = ax;
5153  }
5154  else if (1.0 / ratio > max_ratio)
5155  {
5156  max_ratio = 1.0 / ratio;
5157  index = next_ax;
5158  }
5159  }
5160  return std::make_pair(index, max_ratio);
5161  }
5162 
5163 
5164  template <int dim, int spacedim>
5165  void
5167  const bool isotropic,
5168  const unsigned int max_iterations)
5169  {
5170  unsigned int iter = 0;
5171  bool continue_refinement = true;
5172 
5173  while (continue_refinement && (iter < max_iterations))
5174  {
5175  if (max_iterations != numbers::invalid_unsigned_int)
5176  iter++;
5177  continue_refinement = false;
5178 
5179  for (const auto &cell : tria.active_cell_iterators())
5180  for (const unsigned int j : cell->face_indices())
5181  if (cell->at_boundary(j) == false &&
5182  cell->neighbor(j)->has_children())
5183  {
5184  if (isotropic)
5185  {
5186  cell->set_refine_flag();
5187  continue_refinement = true;
5188  }
5189  else
5190  continue_refinement |= cell->flag_for_face_refinement(j);
5191  }
5192 
5194  }
5195  }
5196 
5197  template <int dim, int spacedim>
5198  void
5200  const double max_ratio,
5201  const unsigned int max_iterations)
5202  {
5203  unsigned int iter = 0;
5204  bool continue_refinement = true;
5205 
5206  while (continue_refinement && (iter < max_iterations))
5207  {
5208  iter++;
5209  continue_refinement = false;
5210  for (const auto &cell : tria.active_cell_iterators())
5211  {
5212  std::pair<unsigned int, double> info =
5213  GridTools::get_longest_direction<dim, spacedim>(cell);
5214  if (info.second > max_ratio)
5215  {
5216  cell->set_refine_flag(
5217  RefinementCase<dim>::cut_axis(info.first));
5218  continue_refinement = true;
5219  }
5220  }
5222  }
5223  }
5224 
5225 
5226  template <int dim, int spacedim>
5227  void
5229  const double limit_angle_fraction)
5230  {
5231  if (dim == 1)
5232  return; // Nothing to do
5233 
5234  // Check that we don't have hanging nodes
5236  ExcMessage("The input Triangulation cannot "
5237  "have hanging nodes."));
5238 
5240 
5241  bool has_cells_with_more_than_dim_faces_on_boundary = true;
5242  bool has_cells_with_dim_faces_on_boundary = false;
5243 
5244  unsigned int refinement_cycles = 0;
5245 
5246  while (has_cells_with_more_than_dim_faces_on_boundary)
5247  {
5248  has_cells_with_more_than_dim_faces_on_boundary = false;
5249 
5250  for (const auto &cell : tria.active_cell_iterators())
5251  {
5252  unsigned int boundary_face_counter = 0;
5253  for (auto f : cell->face_indices())
5254  if (cell->face(f)->at_boundary())
5255  boundary_face_counter++;
5256  if (boundary_face_counter > dim)
5257  {
5258  has_cells_with_more_than_dim_faces_on_boundary = true;
5259  break;
5260  }
5261  else if (boundary_face_counter == dim)
5262  has_cells_with_dim_faces_on_boundary = true;
5263  }
5264  if (has_cells_with_more_than_dim_faces_on_boundary)
5265  {
5266  tria.refine_global(1);
5267  refinement_cycles++;
5268  }
5269  }
5270 
5271  if (has_cells_with_dim_faces_on_boundary)
5272  {
5273  tria.refine_global(1);
5274  refinement_cycles++;
5275  }
5276  else
5277  {
5278  while (refinement_cycles > 0)
5279  {
5280  for (const auto &cell : tria.active_cell_iterators())
5281  cell->set_coarsen_flag();
5283  refinement_cycles--;
5284  }
5285  return;
5286  }
5287 
5288  std::vector<bool> cells_to_remove(tria.n_active_cells(), false);
5289  std::vector<Point<spacedim>> vertices = tria.get_vertices();
5290 
5291  std::vector<bool> faces_to_remove(tria.n_raw_faces(), false);
5292 
5293  std::vector<CellData<dim>> cells_to_add;
5294  SubCellData subcelldata_to_add;
5295 
5296  // Trick compiler for dimension independent things
5297  const unsigned int v0 = 0, v1 = 1, v2 = (dim > 1 ? 2 : 0),
5298  v3 = (dim > 1 ? 3 : 0);
5299 
5300  for (const auto &cell : tria.active_cell_iterators())
5301  {
5302  double angle_fraction = 0;
5303  unsigned int vertex_at_corner = numbers::invalid_unsigned_int;
5304 
5305  if (dim == 2)
5306  {
5308  p0[spacedim > 1 ? 1 : 0] = 1;
5310  p1[0] = 1;
5311 
5312  if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
5313  {
5314  p0 = cell->vertex(v0) - cell->vertex(v2);
5315  p1 = cell->vertex(v3) - cell->vertex(v2);
5316  vertex_at_corner = v2;
5317  }
5318  else if (cell->face(v3)->at_boundary() &&
5319  cell->face(v1)->at_boundary())
5320  {
5321  p0 = cell->vertex(v2) - cell->vertex(v3);
5322  p1 = cell->vertex(v1) - cell->vertex(v3);
5323  vertex_at_corner = v3;
5324  }
5325  else if (cell->face(1)->at_boundary() &&
5326  cell->face(2)->at_boundary())
5327  {
5328  p0 = cell->vertex(v0) - cell->vertex(v1);
5329  p1 = cell->vertex(v3) - cell->vertex(v1);
5330  vertex_at_corner = v1;
5331  }
5332  else if (cell->face(2)->at_boundary() &&
5333  cell->face(0)->at_boundary())
5334  {
5335  p0 = cell->vertex(v2) - cell->vertex(v0);
5336  p1 = cell->vertex(v1) - cell->vertex(v0);
5337  vertex_at_corner = v0;
5338  }
5339  p0 /= p0.norm();
5340  p1 /= p1.norm();
5341  angle_fraction = std::acos(p0 * p1) / numbers::PI;
5342  }
5343  else
5344  {
5345  Assert(false, ExcNotImplemented());
5346  }
5347 
5348  if (angle_fraction > limit_angle_fraction)
5349  {
5350  auto flags_removal = [&](unsigned int f1,
5351  unsigned int f2,
5352  unsigned int n1,
5353  unsigned int n2) -> void {
5354  cells_to_remove[cell->active_cell_index()] = true;
5355  cells_to_remove[cell->neighbor(n1)->active_cell_index()] = true;
5356  cells_to_remove[cell->neighbor(n2)->active_cell_index()] = true;
5357 
5358  faces_to_remove[cell->face(f1)->index()] = true;
5359  faces_to_remove[cell->face(f2)->index()] = true;
5360 
5361  faces_to_remove[cell->neighbor(n1)->face(f1)->index()] = true;
5362  faces_to_remove[cell->neighbor(n2)->face(f2)->index()] = true;
5363  };
5364 
5365  auto cell_creation = [&](const unsigned int vv0,
5366  const unsigned int vv1,
5367  const unsigned int f0,
5368  const unsigned int f1,
5369 
5370  const unsigned int n0,
5371  const unsigned int v0n0,
5372  const unsigned int v1n0,
5373 
5374  const unsigned int n1,
5375  const unsigned int v0n1,
5376  const unsigned int v1n1) {
5377  CellData<dim> c1, c2;
5378  CellData<1> l1, l2;
5379 
5380  c1.vertices[v0] = cell->vertex_index(vv0);
5381  c1.vertices[v1] = cell->vertex_index(vv1);
5382  c1.vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
5383  c1.vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
5384 
5385  c1.manifold_id = cell->manifold_id();
5386  c1.material_id = cell->material_id();
5387 
5388  c2.vertices[v0] = cell->vertex_index(vv0);
5389  c2.vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
5390  c2.vertices[v2] = cell->vertex_index(vv1);
5391  c2.vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
5392 
5393  c2.manifold_id = cell->manifold_id();
5394  c2.material_id = cell->material_id();
5395 
5396  l1.vertices[0] = cell->vertex_index(vv0);
5397  l1.vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
5398 
5399  l1.boundary_id = cell->line(f0)->boundary_id();
5400  l1.manifold_id = cell->line(f0)->manifold_id();
5401  subcelldata_to_add.boundary_lines.push_back(l1);
5402 
5403  l2.vertices[0] = cell->vertex_index(vv0);
5404  l2.vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
5405 
5406  l2.boundary_id = cell->line(f1)->boundary_id();
5407  l2.manifold_id = cell->line(f1)->manifold_id();
5408  subcelldata_to_add.boundary_lines.push_back(l2);
5409 
5410  cells_to_add.push_back(c1);
5411  cells_to_add.push_back(c2);
5412  };
5413 
5414  if (dim == 2)
5415  {
5416  switch (vertex_at_corner)
5417  {
5418  case 0:
5419  flags_removal(0, 2, 3, 1);
5420  cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
5421  break;
5422  case 1:
5423  flags_removal(1, 2, 3, 0);
5424  cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
5425  break;
5426  case 2:
5427  flags_removal(3, 0, 1, 2);
5428  cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
5429  break;
5430  case 3:
5431  flags_removal(3, 1, 0, 2);
5432  cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
5433  break;
5434  }
5435  }
5436  else
5437  {
5438  Assert(false, ExcNotImplemented());
5439  }
5440  }
5441  }
5442 
5443  // if no cells need to be added, then no regularization is necessary.
5444  // Restore things as they were before this function was called.
5445  if (cells_to_add.size() == 0)
5446  {
5447  while (refinement_cycles > 0)
5448  {
5449  for (const auto &cell : tria.active_cell_iterators())
5450  cell->set_coarsen_flag();
5452  refinement_cycles--;
5453  }
5454  return;
5455  }
5456 
5457  // add the cells that were not marked as skipped
5458  for (const auto &cell : tria.active_cell_iterators())
5459  {
5460  if (cells_to_remove[cell->active_cell_index()] == false)
5461  {
5462  CellData<dim> c(cell->n_vertices());
5463  for (const unsigned int v : cell->vertex_indices())
5464  c.vertices[v] = cell->vertex_index(v);
5465  c.manifold_id = cell->manifold_id();
5466  c.material_id = cell->material_id();
5467  cells_to_add.push_back(c);
5468  }
5469  }
5470 
5471  // Face counter for both dim == 2 and dim == 3
5473  face = tria.begin_active_face(),
5474  endf = tria.end_face();
5475  for (; face != endf; ++face)
5476  if ((face->at_boundary() ||
5477  face->manifold_id() != numbers::flat_manifold_id) &&
5478  faces_to_remove[face->index()] == false)
5479  {
5480  for (unsigned int l = 0; l < face->n_lines(); ++l)
5481  {
5482  CellData<1> line;
5483  if (dim == 2)
5484  {
5485  for (const unsigned int v : face->vertex_indices())
5486  line.vertices[v] = face->vertex_index(v);
5487  line.boundary_id = face->boundary_id();
5488  line.manifold_id = face->manifold_id();
5489  }
5490  else
5491  {
5492  for (const unsigned int v : face->line(l)->vertex_indices())
5493  line.vertices[v] = face->line(l)->vertex_index(v);
5494  line.boundary_id = face->line(l)->boundary_id();
5495  line.manifold_id = face->line(l)->manifold_id();
5496  }
5497  subcelldata_to_add.boundary_lines.push_back(line);
5498  }
5499  if (dim == 3)
5500  {
5501  CellData<2> quad(face->n_vertices());
5502  for (const unsigned int v : face->vertex_indices())
5503  quad.vertices[v] = face->vertex_index(v);
5504  quad.boundary_id = face->boundary_id();
5505  quad.manifold_id = face->manifold_id();
5506  subcelldata_to_add.boundary_quads.push_back(quad);
5507  }
5508  }
5510  cells_to_add,
5511  subcelldata_to_add);
5513 
5514  // Save manifolds
5515  auto manifold_ids = tria.get_manifold_ids();
5516  std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
5517  manifolds;
5518  // Set manifolds in new Triangulation
5519  for (const auto manifold_id : manifold_ids)
5521  manifolds[manifold_id] = tria.get_manifold(manifold_id).clone();
5522 
5523  tria.clear();
5524 
5525  tria.create_triangulation(vertices, cells_to_add, subcelldata_to_add);
5526 
5527  // Restore manifolds
5528  for (const auto manifold_id : manifold_ids)
5530  tria.set_manifold(manifold_id, *manifolds[manifold_id]);
5531  }
5532 
5533 
5534 
5535  template <int dim, int spacedim>
5536 #ifndef DOXYGEN
5537  std::tuple<
5538  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5539  std::vector<std::vector<Point<dim>>>,
5540  std::vector<std::vector<unsigned int>>>
5541 #else
5542  return_type
5543 #endif
5545  const Cache<dim, spacedim> & cache,
5546  const std::vector<Point<spacedim>> &points,
5548  &cell_hint)
5549  {
5550  const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
5551  // Splitting the tuple's components
5552  auto &cells = std::get<0>(cqmp);
5553  auto &qpoints = std::get<1>(cqmp);
5554  auto &maps = std::get<2>(cqmp);
5555 
5556  return std::make_tuple(std::move(cells),
5557  std::move(qpoints),
5558  std::move(maps));
5559  }
5560 
5561 
5562 
5563  template <int dim, int spacedim>
5564 #ifndef DOXYGEN
5565  std::tuple<
5566  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5567  std::vector<std::vector<Point<dim>>>,
5568  std::vector<std::vector<unsigned int>>,
5569  std::vector<unsigned int>>
5570 #else
5571  return_type
5572 #endif
5574  const Cache<dim, spacedim> & cache,
5575  const std::vector<Point<spacedim>> &points,
5577  &cell_hint)
5578  {
5579  Assert((dim == spacedim),
5580  ExcMessage("Only implemented for dim==spacedim."));
5581 
5582  // Alias
5583  namespace bgi = boost::geometry::index;
5584 
5585  // Get the mapping
5586  const auto &mapping = cache.get_mapping();
5587 
5588  // How many points are here?
5589  const unsigned int np = points.size();
5590 
5591  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5592  cells_out;
5593  std::vector<std::vector<Point<dim>>> qpoints_out;
5594  std::vector<std::vector<unsigned int>> maps_out;
5595  std::vector<unsigned int> missing_points_out;
5596 
5597  // Now the easy case.
5598  if (np == 0)
5599  return std::make_tuple(std::move(cells_out),
5600  std::move(qpoints_out),
5601  std::move(maps_out),
5602  std::move(missing_points_out));
5603 
5604  // For the search we shall use the following tree
5605  const auto &b_tree = cache.get_cell_bounding_boxes_rtree();
5606 
5607  // Now make a tree of indices for the points
5608  // [TODO] This would work better with pack_rtree_of_indices, but
5609  // windows does not like it. Build a tree with pairs of point and id
5610  std::vector<std::pair<Point<spacedim>, unsigned int>> points_and_ids(np);
5611  for (unsigned int i = 0; i < np; ++i)
5612  points_and_ids[i] = std::make_pair(points[i], i);
5613  const auto p_tree = pack_rtree(points_and_ids);
5614 
5615  // Keep track of all found points
5616  std::vector<bool> found_points(points.size(), false);
5617 
5618  // Check if a point was found
5619  const auto already_found = [&found_points](const auto &id) {
5620  AssertIndexRange(id.second, found_points.size());
5621  return found_points[id.second];
5622  };
5623 
5624  // check if the given cell was already in the vector of cells before. If so,
5625  // insert in the corresponding vectors the reference point and the id.
5626  // Otherwise append a new entry to all vectors.
5627  const auto store_cell_point_and_id =
5628  [&](
5630  const Point<dim> & ref_point,
5631  const unsigned int &id) {
5632  const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
5633  if (it != cells_out.rend())
5634  {
5635  const auto cell_id =
5636  (cells_out.size() - 1 - (it - cells_out.rbegin()));
5637  qpoints_out[cell_id].emplace_back(ref_point);
5638  maps_out[cell_id].emplace_back(id);
5639  }
5640  else
5641  {
5642  cells_out.emplace_back(cell);
5643  qpoints_out.emplace_back(std::vector<Point<dim>>({ref_point}));
5644  maps_out.emplace_back(std::vector<unsigned int>({id}));
5645  }
5646  };
5647 
5648  // Check all points within a given pair of box and cell
5649  const auto check_all_points_within_box = [&](const auto &leaf) {
5650  const auto &box = leaf.first;
5651  const auto &cell_hint = leaf.second;
5652 
5653  for (const auto &point_and_id :
5654  p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
5655  bgi::intersects(box)))
5656  {
5657  const auto id = point_and_id.second;
5658  const auto cell_and_ref =
5660  points[id],
5661  cell_hint);
5662  const auto &cell = cell_and_ref.first;
5663  const auto &ref_point = cell_and_ref.second;
5664 
5665  if (cell.state() == IteratorState::valid)
5666  store_cell_point_and_id(cell, ref_point, id);
5667  else
5668  missing_points_out.emplace_back(id);
5669 
5670  // Don't look anymore for this point
5671  found_points[id] = true;
5672  }
5673  };
5674 
5675  // If a hint cell was given, use it
5676  if (cell_hint.state() == IteratorState::valid)
5677  check_all_points_within_box(
5678  std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
5679 
5680  // Now loop over all points that have not been found yet
5681  for (unsigned int i = 0; i < np; ++i)
5682  if (found_points[i] == false)
5683  {
5684  // Get the closest cell to this point
5685  const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
5686  // Now checks all points that fall within this box
5687  if (leaf != b_tree.qend())
5688  check_all_points_within_box(*leaf);
5689  else
5690  {
5691  // We should not get here. Throw an error.
5692  Assert(false, ExcInternalError());
5693  }
5694  }
5695  // Now make sure we send out the rest of the points that we did not find.
5696  for (unsigned int i = 0; i < np; ++i)
5697  if (found_points[i] == false)
5698  missing_points_out.emplace_back(i);
5699 
5700  // Debug Checking
5701  AssertDimension(cells_out.size(), maps_out.size());
5702  AssertDimension(cells_out.size(), qpoints_out.size());
5703 
5704 #ifdef DEBUG
5705  unsigned int c = cells_out.size();
5706  unsigned int qps = 0;
5707  // The number of points in all
5708  // the cells must be the same as
5709  // the number of points we
5710  // started off from,
5711  // plus the points which were ignored
5712  for (unsigned int n = 0; n < c; ++n)
5713  {
5714  AssertDimension(qpoints_out[n].size(), maps_out[n].size());
5715  qps += qpoints_out[n].size();
5716  }
5717 
5718  Assert(qps + missing_points_out.size() == np,
5719  ExcDimensionMismatch(qps + missing_points_out.size(), np));
5720 #endif
5721 
5722  return std::make_tuple(std::move(cells_out),
5723  std::move(qpoints_out),
5724  std::move(maps_out),
5725  std::move(missing_points_out));
5726  }
5727 
5728 
5729 
5730  template <int dim, int spacedim>
5731 #ifndef DOXYGEN
5732  std::tuple<
5733  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5734  std::vector<std::vector<Point<dim>>>,
5735  std::vector<std::vector<unsigned int>>,
5736  std::vector<std::vector<Point<spacedim>>>,
5737  std::vector<std::vector<unsigned int>>>
5738 #else
5739  return_type
5740 #endif
5742  const GridTools::Cache<dim, spacedim> & cache,
5743  const std::vector<Point<spacedim>> & points,
5744  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5745  const double tolerance)
5746  {
5747  // run internal function ...
5749  cache, points, global_bboxes, {}, tolerance, false, true)
5750  .send_components;
5751 
5752  // ... and reshuffle the data
5753  std::tuple<
5754  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5755  std::vector<std::vector<Point<dim>>>,
5756  std::vector<std::vector<unsigned int>>,
5757  std::vector<std::vector<Point<spacedim>>>,
5758  std::vector<std::vector<unsigned int>>>
5759  result;
5760 
5761  std::pair<int, int> dummy{-1, -1};
5762 
5763  for (unsigned int i = 0; i < all.size(); ++i)
5764  {
5765  if (dummy != std::get<0>(all[i]))
5766  {
5767  std::get<0>(result).push_back(
5769  &cache.get_triangulation(),
5770  std::get<0>(all[i]).first,
5771  std::get<0>(all[i]).second});
5772 
5773  const unsigned int new_size = std::get<0>(result).size();
5774 
5775  std::get<1>(result).resize(new_size);
5776  std::get<2>(result).resize(new_size);
5777  std::get<3>(result).resize(new_size);
5778  std::get<4>(result).resize(new_size);
5779 
5780  dummy = std::get<0>(all[i]);
5781  }
5782 
5783  std::get<1>(result).back().push_back(
5784  std::get<3>(all[i])); // reference point
5785  std::get<2>(result).back().push_back(std::get<2>(all[i])); // index
5786  std::get<3>(result).back().push_back(std::get<4>(all[i])); // real point
5787  std::get<4>(result).back().push_back(std::get<1>(all[i])); // rank
5788  }
5789 
5790  return result;
5791  }
5792 
5793 
5794 
5795  namespace internal
5796  {
5797  template <int spacedim>
5798  std::tuple<std::vector<unsigned int>,
5799  std::vector<unsigned int>,
5800  std::vector<unsigned int>>
5802  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5803  const std::vector<Point<spacedim>> & points,
5804  const double tolerance)
5805  {
5806  std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
5807  ranks_and_indices.reserve(points.size());
5808 
5809  for (unsigned int i = 0; i < points.size(); ++i)
5810  {
5811  const auto &point = points[i];
5812  for (unsigned rank = 0; rank < global_bboxes.size(); ++rank)
5813  for (const auto &box : global_bboxes[rank])
5814  if (box.point_inside(point, tolerance))
5815  {
5816  ranks_and_indices.emplace_back(rank, i);
5817  break;
5818  }
5819  }
5820 
5821  // convert to CRS
5822  std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
5823 
5824  std::vector<unsigned int> ranks;
5825  std::vector<unsigned int> ptr;
5826  std::vector<unsigned int> indices;
5827 
5828  unsigned int dummy_rank = numbers::invalid_unsigned_int;
5829 
5830  for (const auto &i : ranks_and_indices)
5831  {
5832  if (dummy_rank != i.first)
5833  {
5834  dummy_rank = i.first;
5835  ranks.push_back(dummy_rank);
5836  ptr.push_back(indices.size());
5837  }
5838 
5839  indices.push_back(i.second);
5840  }
5841  ptr.push_back(indices.size());
5842 
5843  return std::make_tuple(std::move(ranks),
5844  std::move(ptr),
5845  std::move(indices));
5846  }
5847 
5848 
5849 
5850  template <int dim, int spacedim>
5851  std::vector<
5852  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5853  Point<dim>>>
5855  const Cache<dim, spacedim> & cache,
5856  const Point<spacedim> & point,
5858  const std::vector<bool> &marked_vertices,
5859  const double tolerance,
5860  const bool enforce_unique_mapping)
5861  {
5862  std::vector<
5863  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5864  Point<dim>>>
5865  locally_owned_active_cells_around_point;
5866 
5867  const auto first_cell = GridTools::find_active_cell_around_point(
5868  cache.get_mapping(),
5869  cache.get_triangulation(),
5870  point,
5871  cache.get_vertex_to_cell_map(),
5873  cell_hint,
5874  marked_vertices,
5875  cache.get_used_vertices_rtree(),
5876  tolerance,
5878 
5879  const unsigned int my_rank = Utilities::MPI::this_mpi_process(
5881 
5882  cell_hint = first_cell.first;
5883  if (cell_hint.state() == IteratorState::valid)
5884  {
5885  const auto active_cells_around_point =
5887  cache.get_mapping(),
5888  cache.get_triangulation(),
5889  point,
5890  tolerance,
5891  first_cell);
5892 
5893  if (enforce_unique_mapping)
5894  {
5895  // check if the rank of this process is the lowest of all cells
5896  // if not, the other process will handle this cell and we don't
5897  // have to do here anything in the case of unique mapping
5898  unsigned int lowes_rank = numbers::invalid_unsigned_int;
5899 
5900  for (const auto &cell : active_cells_around_point)
5901  lowes_rank = std::min(lowes_rank, cell.first->subdomain_id());
5902 
5903  if (lowes_rank != my_rank)
5904  return {};
5905  }
5906 
5907  locally_owned_active_cells_around_point.reserve(
5908  active_cells_around_point.size());
5909 
5910  for (const auto &cell : active_cells_around_point)
5911  if (cell.first->is_locally_owned())
5912  locally_owned_active_cells_around_point.push_back(cell);
5913  }
5914 
5915  std::sort(locally_owned_active_cells_around_point.begin(),
5916  locally_owned_active_cells_around_point.end(),
5917  [](const auto &a, const auto &b) { return a.first < b.first; });
5918 
5919  if (enforce_unique_mapping &&
5920  locally_owned_active_cells_around_point.size() > 1)
5921  // in the case of unique mapping, we only need a single cell
5922  return {locally_owned_active_cells_around_point.front()};
5923  else
5924  return locally_owned_active_cells_around_point;
5925  }
5926 
5927 
5928 
5929  template <int dim, int spacedim>
5930  DistributedComputePointLocationsInternal<dim, spacedim>
5932  const GridTools::Cache<dim, spacedim> & cache,
5933  const std::vector<Point<spacedim>> & points,
5934  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5935  const std::vector<bool> & marked_vertices,
5936  const double tolerance,
5937  const bool perform_handshake,
5938  const bool enforce_unique_mapping)
5939  {
5941 
5942  auto &send_components = result.send_components;
5943  auto &send_ranks = result.send_ranks;
5944  auto &send_ptrs = result.send_ptrs;
5945  auto &recv_components = result.recv_components;
5946  auto &recv_ranks = result.recv_ranks;
5947  auto &recv_ptrs = result.recv_ptrs;
5948 
5949  const auto potential_owners =
5950  internal::guess_point_owner(global_bboxes, points, tolerance);
5951 
5952  const auto &potential_owners_ranks = std::get<0>(potential_owners);
5953  const auto &potential_owners_ptrs = std::get<1>(potential_owners);
5954  const auto &potential_owners_indices = std::get<2>(potential_owners);
5955 
5956  auto cell_hint = cache.get_triangulation().begin_active();
5957 
5958  const auto translate = [&](const unsigned int other_rank) {
5959  const auto ptr = std::find(potential_owners_ranks.begin(),
5960  potential_owners_ranks.end(),
5961  other_rank);
5962 
5963  Assert(ptr != potential_owners_ranks.end(), ExcInternalError());
5964 
5965  const auto other_rank_index =
5966  std::distance(potential_owners_ranks.begin(), ptr);
5967 
5968  return other_rank_index;
5969  };
5970 
5971  Assert(
5972  (marked_vertices.size() == 0) ||
5973  (marked_vertices.size() == cache.get_triangulation().n_vertices()),
5974  ExcMessage(
5975  "The marked_vertices vector has to be either empty or its size has "
5976  "to equal the number of vertices of the triangulation."));
5977 
5978  using RequestType = std::vector<std::pair<unsigned int, Point<spacedim>>>;
5979  using AnswerType = std::vector<unsigned int>;
5980 
5981  // In the case that a marked_vertices vector has been given and none
5982  // of its entries is true, we know that this process does not own
5983  // any of the incoming points (and it will not send any data) so
5984  // that we can take a short cut.
5985  const bool has_relevant_vertices =
5986  (marked_vertices.size() == 0) ||
5987  (std::find(marked_vertices.begin(), marked_vertices.end(), true) !=
5988  marked_vertices.end());
5989 
5990  const auto create_request = [&](const unsigned int other_rank) {
5991  const auto other_rank_index = translate(other_rank);
5992 
5993  RequestType request;
5994  request.reserve(potential_owners_ptrs[other_rank_index + 1] -
5995  potential_owners_ptrs[other_rank_index]);
5996 
5997  for (unsigned int i = potential_owners_ptrs[other_rank_index];
5998  i < potential_owners_ptrs[other_rank_index + 1];
5999  ++i)
6000  request.emplace_back(potential_owners_indices[i],
6001  points[potential_owners_indices[i]]);
6002 
6003  return request;
6004  };
6005 
6006  const auto answer_request =
6007  [&](const unsigned int &other_rank,
6008  const RequestType & request) -> AnswerType {
6009  AnswerType answer(request.size(), 0);
6010 
6011  if (has_relevant_vertices)
6012  {
6013  cell_hint = cache.get_triangulation().begin_active();
6014 
6015  for (unsigned int i = 0; i < request.size(); ++i)
6016  {
6017  const auto &index_and_point = request[i];
6018 
6019  const auto cells_and_reference_positions =
6021  cache,
6022  index_and_point.second,
6023  cell_hint,
6024  marked_vertices,
6025  tolerance,
6026  enforce_unique_mapping);
6027 
6028  for (const auto &cell_and_reference_position :
6029  cells_and_reference_positions)
6030  {
6031  send_components.emplace_back(
6032  std::pair<int, int>(
6033  cell_and_reference_position.first->level(),
6034  cell_and_reference_position.first->index()),
6035  other_rank,
6036  index_and_point.first,
6037  cell_and_reference_position.second,
6038  index_and_point.second,
6040  }
6041 
6042  answer[i] = cells_and_reference_positions.size();
6043  }
6044  }
6045 
6046  if (perform_handshake)
6047  return answer;
6048  else
6049  return {};
6050  };
6051 
6052  const auto process_answer = [&](const unsigned int other_rank,
6053  const AnswerType & answer) {
6054  if (perform_handshake)
6055  {
6056  const auto other_rank_index = translate(other_rank);
6057 
6058  for (unsigned int i = 0; i < answer.size(); ++i)
6059  for (unsigned int j = 0; j < answer[i]; ++j)
6060  recv_components.emplace_back(
6061  other_rank,
6062  potential_owners_indices
6063  [i + potential_owners_ptrs[other_rank_index]],
6065  }
6066  };
6067 
6068  Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
6069  potential_owners_ranks,
6070  create_request,
6071  answer_request,
6072  process_answer,
6074 
6075  if (true)
6076  {
6077  // sort according to rank (and point index and cell) -> make
6078  // deterministic
6079  std::sort(send_components.begin(),
6080  send_components.end(),
6081  [&](const auto &a, const auto &b) {
6082  if (std::get<1>(a) != std::get<1>(b)) // rank
6083  return std::get<1>(a) < std::get<1>(b);
6084 
6085  if (std::get<2>(a) != std::get<2>(b)) // point index
6086  return std::get<2>(a) < std::get<2>(b);
6087 
6088  return std::get<0>(a) < std::get<0>(b); // cell
6089  });
6090 
6091  // perform enumeration and extract rank information
6092  for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
6093  i < send_components.size();
6094  ++i)
6095  {
6096  std::get<5>(send_components[i]) = i;
6097 
6098  if (dummy != std::get<1>(send_components[i]))
6099  {
6100  dummy = std::get<1>(send_components[i]);
6101  send_ranks.push_back(dummy);
6102  send_ptrs.push_back(i);
6103  }
6104  }
6105  send_ptrs.push_back(send_components.size());
6106 
6107  // sort according to cell, rank, point index (while keeping
6108  // partial ordering)
6109  std::sort(send_components.begin(),
6110  send_components.end(),
6111  [&](const auto &a, const auto &b) {
6112  if (std::get<0>(a) != std::get<0>(b))
6113  return std::get<0>(a) < std::get<0>(b); // cell
6114 
6115  if (std::get<1>(a) != std::get<1>(b))
6116  return std::get<1>(a) < std::get<1>(b); // rank
6117 
6118  if (std::get<2>(a) != std::get<2>(b))
6119  return std::get<2>(a) < std::get<2>(b); // point index
6120 
6121  return std::get<5>(a) < std::get<5>(b); // enumeration
6122  });
6123  }
6124 
6125  if (perform_handshake)
6126  {
6127  // sort according to rank (and point index) -> make deterministic
6128  std::sort(recv_components.begin(),
6129  recv_components.end(),
6130  [&](const auto &a, const auto &b) {
6131  if (std::get<0>(a) != std::get<0>(b))
6132  return std::get<0>(a) < std::get<0>(b); // rank
6133 
6134  return std::get<1>(a) < std::get<1>(b); // point index
6135  });
6136 
6137  // perform enumeration and extract rank information
6138  for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
6139  i < recv_components.size();
6140  ++i)
6141  {
6142  std::get<2>(recv_components[i]) = i;
6143 
6144  if (dummy != std::get<0>(recv_components[i]))
6145  {
6146  dummy = std::get<0>(recv_components[i]);
6147  recv_ranks.push_back(dummy);
6148  recv_ptrs.push_back(i);
6149  }
6150  }
6151  recv_ptrs.push_back(recv_components.size());
6152 
6153  // sort according to point index and rank (while keeping partial
6154  // ordering)
6155  std::sort(recv_components.begin(),
6156  recv_components.end(),
6157  [&](const auto &a, const auto &b) {
6158  if (std::get<1>(a) != std::get<1>(b))
6159  return std::get<1>(a) < std::get<1>(b); // point index
6160 
6161  if (std::get<0>(a) != std::get<0>(b))
6162  return std::get<0>(a) < std::get<0>(b); // rank
6163 
6164  return std::get<2>(a) < std::get<2>(b); // enumeration
6165  });
6166  }
6167 
6168  return result;
6169  }
6170  } // namespace internal
6171 
6172 
6173 
6174  template <int dim, int spacedim>
6175  std::map<unsigned int, Point<spacedim>>
6177  const Mapping<dim, spacedim> & mapping)
6178  {
6179  std::map<unsigned int, Point<spacedim>> result;
6180  for (const auto &cell : container.active_cell_iterators())
6181  {
6182  if (!cell->is_artificial())
6183  {
6184  const auto vs = mapping.get_vertices(cell);
6185  for (unsigned int i = 0; i < vs.size(); ++i)
6186  result[cell->vertex_index(i)] = vs[i];
6187  }
6188  }
6189  return result;
6190  }
6191 
6192 
6193  template <int spacedim>
6194  unsigned int
6195  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
6196  const Point<spacedim> & p)
6197  {
6198  auto id_and_v = std::min_element(
6199  vertices.begin(),
6200  vertices.end(),
6201  [&](const std::pair<const unsigned int, Point<spacedim>> &p1,
6202  const std::pair<const unsigned int, Point<spacedim>> &p2) -> bool {
6203  return p1.second.distance(p) < p2.second.distance(p);
6204  });
6205  return id_and_v->first;
6206  }
6207 
6208 
6209  template <int dim, int spacedim>
6210  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
6211  Point<dim>>
6213  const Cache<dim, spacedim> &cache,
6214  const Point<spacedim> & p,
6216  & cell_hint,
6217  const std::vector<bool> &marked_vertices,
6218  const double tolerance)
6219  {
6220  const auto &mesh = cache.get_triangulation();
6221  const auto &mapping = cache.get_mapping();
6222  const auto &vertex_to_cells = cache.get_vertex_to_cell_map();
6223  const auto &vertex_to_cell_centers =
6225  const auto &used_vertices_rtree = cache.get_used_vertices_rtree();
6226 
6227  return find_active_cell_around_point(mapping,
6228  mesh,
6229  p,
6230  vertex_to_cells,
6231  vertex_to_cell_centers,
6232  cell_hint,
6233  marked_vertices,
6234  used_vertices_rtree,
6235  tolerance);
6236  }
6237 
6238  template <int spacedim>
6239  std::vector<std::vector<BoundingBox<spacedim>>>
6241  const std::vector<BoundingBox<spacedim>> &local_bboxes,
6242  const MPI_Comm & mpi_communicator)
6243  {
6244 #ifndef DEAL_II_WITH_MPI
6245  (void)local_bboxes;
6246  (void)mpi_communicator;
6247  Assert(false,
6248  ExcMessage(
6249  "GridTools::exchange_local_bounding_boxes() requires MPI."));
6250  return {};
6251 #else
6252  // Step 1: preparing data to be sent
6253  unsigned int n_bboxes = local_bboxes.size();
6254  // Dimension of the array to be exchanged (number of double)
6255  int n_local_data = 2 * spacedim * n_bboxes;
6256  // data array stores each entry of each point describing the bounding boxes
6257  std::vector<double> loc_data_array(n_local_data);
6258  for (unsigned int i = 0; i < n_bboxes; ++i)
6259  for (unsigned int d = 0; d < spacedim; ++d)
6260  {
6261  // Extracting the coordinates of each boundary point
6262  loc_data_array[2 * i * spacedim + d] =
6263  local_bboxes[i].get_boundary_points().first[d];
6264  loc_data_array[2 * i * spacedim + spacedim + d] =
6265  local_bboxes[i].get_boundary_points().second[d];
6266  }
6267 
6268  // Step 2: exchanging the size of local data
6269  unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_communicator);
6270 
6271  // Vector to store the size of loc_data_array for every process
6272  std::vector<int> size_all_data(n_procs);
6273 
6274  // Exchanging the number of bboxes
6275  int ierr = MPI_Allgather(&n_local_data,
6276  1,
6277  MPI_INT,
6278  size_all_data.data(),
6279  1,
6280  MPI_INT,
6281  mpi_communicator);
6282  AssertThrowMPI(ierr);
6283 
6284  // Now computing the displacement, relative to recvbuf,
6285  // at which to store the incoming data
6286  std::vector<int> rdispls(n_procs);
6287  rdispls[0] = 0;
6288  for (unsigned int i = 1; i < n_procs; ++i)
6289  rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
6290 
6291  // Step 3: exchange the data and bounding boxes:
6292  // Allocating a vector to contain all the received data
6293  std::vector<double> data_array(rdispls.back() + size_all_data.back());
6294 
6295  ierr = MPI_Allgatherv(loc_data_array.data(),
6296  n_local_data,
6297  MPI_DOUBLE,
6298  data_array.data(),
6299  size_all_data.data(),
6300  rdispls.data(),
6301  MPI_DOUBLE,
6302  mpi_communicator);
6303  AssertThrowMPI(ierr);
6304 
6305  // Step 4: create the array of bboxes for output
6306  std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
6307  unsigned int begin_idx = 0;
6308  for (unsigned int i = 0; i < n_procs; ++i)
6309  {
6310  // Number of local bounding boxes
6311  unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
6312  global_bboxes[i].resize(n_bbox_i);
6313  for (unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
6314  {
6315  Point<spacedim> p1, p2; // boundary points for bbox
6316  for (unsigned int d = 0; d < spacedim; ++d)
6317  {
6318  p1[d] = data_array[begin_idx + 2 * bbox * spacedim + d];
6319  p2[d] =
6320  data_array[begin_idx + 2 * bbox * spacedim + spacedim + d];
6321  }
6322  BoundingBox<spacedim> loc_bbox(std::make_pair(p1, p2));
6323  global_bboxes[i][bbox] = loc_bbox;
6324  }
6325  // Shifting the first index to the start of the next vector
6326  begin_idx += size_all_data[i];
6327  }
6328  return global_bboxes;
6329 #endif // DEAL_II_WITH_MPI
6330  }
6331 
6332 
6333 
6334  template <int spacedim>
6337  const std::vector<BoundingBox<spacedim>> &local_description,
6338  const MPI_Comm & mpi_communicator)
6339  {
6340 #ifndef DEAL_II_WITH_MPI
6341  (void)mpi_communicator;
6342  // Building a tree with the only boxes available without MPI
6343  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> boxes_index(
6344  local_description.size());
6345  // Adding to each box the rank of the process owning it
6346  for (unsigned int i = 0; i < local_description.size(); ++i)
6347  boxes_index[i] = std::make_pair(local_description[i], 0u);
6348  return pack_rtree(boxes_index);
6349 #else
6350  // Exchanging local bounding boxes
6351  const std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes =
6352  Utilities::MPI::all_gather(mpi_communicator, local_description);
6353 
6354  // Preparing to flatten the vector
6355  const unsigned int n_procs =
6356  Utilities::MPI::n_mpi_processes(mpi_communicator);
6357  // The i'th element of the following vector contains the index of the first
6358  // local bounding box from the process of rank i