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