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