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