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