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