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