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