Reference documentation for deal.II version Git d3aed38b93 2021-10-28 13:33:27 +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>
36 
41 #include <deal.II/grid/manifold.h>
42 #include <deal.II/grid/tria.h>
45 
49 #include <deal.II/lac/solver_cg.h>
53 #include <deal.II/lac/vector.h>
55 
58 
60 
62 #include <boost/random/mersenne_twister.hpp>
63 #include <boost/random/uniform_real_distribution.hpp>
65 
66 #include <array>
67 #include <cmath>
68 #include <iostream>
69 #include <list>
70 #include <numeric>
71 #include <set>
72 #include <tuple>
73 #include <unordered_map>
74 
76 
77 
78 namespace GridTools
79 {
80  template <int dim, int spacedim>
81  double
83  {
84  // we can't deal with distributed meshes since we don't have all
85  // vertices locally. there is one exception, however: if the mesh has
86  // never been refined. the way to test this is not to ask
87  // tria.n_levels()==1, since this is something that can happen on one
88  // processor without being true on all. however, we can ask for the
89  // global number of active cells and use that
90 #if defined(DEAL_II_WITH_P4EST) && defined(DEBUG)
92  dynamic_cast<
94  Assert(p_tria->n_global_active_cells() == tria.n_cells(0),
96 #endif
97 
98  // the algorithm used simply traverses all cells and picks out the
99  // boundary vertices. it may or may not be faster to simply get all
100  // vectors, don't mark boundary vertices, and compute the distances
101  // thereof, but at least as the mesh is refined, it seems better to
102  // first mark boundary nodes, as marking is O(N) in the number of
103  // cells/vertices, while computing the maximal distance is O(N*N)
104  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
105  std::vector<bool> boundary_vertices(vertices.size(), false);
106 
108  tria.begin_active();
110  tria.end();
111  for (; cell != endc; ++cell)
112  for (const unsigned int face : cell->face_indices())
113  if (cell->face(face)->at_boundary())
114  for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
115  boundary_vertices[cell->face(face)->vertex_index(i)] = true;
116 
117  // now traverse the list of boundary vertices and check distances.
118  // since distances are symmetric, we only have to check one half
119  double max_distance_sqr = 0;
120  std::vector<bool>::const_iterator pi = boundary_vertices.begin();
121  const unsigned int N = boundary_vertices.size();
122  for (unsigned int i = 0; i < N; ++i, ++pi)
123  {
124  std::vector<bool>::const_iterator pj = pi + 1;
125  for (unsigned int j = i + 1; j < N; ++j, ++pj)
126  if ((*pi == true) && (*pj == true) &&
127  ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
128  max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
129  }
130 
131  return std::sqrt(max_distance_sqr);
132  }
133 
134 
135 
136  template <int dim, int spacedim>
137  double
139  const Mapping<dim, spacedim> & mapping)
140  {
141  // get the degree of the mapping if possible. if not, just assume 1
142  unsigned int mapping_degree = 1;
143  if (const auto *p = dynamic_cast<const MappingQ<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
1916  reorient(std::vector<CellData<1>> &)
1917  {}
1918  } // namespace
1919 
1920  template <int dim>
1921  void
1923  {
1924  Assert(cells.size() != 0,
1925  ExcMessage(
1926  "List of elements to orient must have at least one cell"));
1927 
1928  // there is nothing for us to do in 1d
1929  if (dim == 1)
1930  return;
1931 
1932  // check if grids are already consistent. if so, do
1933  // nothing. if not, then do the reordering
1934  if (!is_consistent(cells))
1935  try
1936  {
1937  reorient(cells);
1938  }
1939  catch (const ExcMeshNotOrientable &)
1940  {
1941  // the mesh is not orientable. this is acceptable if we are in 3d,
1942  // as class Triangulation knows how to handle this, but it is
1943  // not in 2d; in that case, re-throw the exception
1944  if (dim < 3)
1945  throw;
1946  }
1947  }
1948 
1949 
1950  // define some transformations
1951  namespace internal
1952  {
1953  template <int spacedim>
1954  class Shift
1955  {
1956  public:
1958  : shift(shift)
1959  {}
1962  {
1963  return p + shift;
1964  }
1965 
1966  private:
1968  };
1969 
1970 
1971  // Transformation to rotate around one of the cartesian axes.
1972  class Rotate3d
1973  {
1974  public:
1975  Rotate3d(const Tensor<1, 3, double> &axis, const double angle)
1976  : rotation_matrix(
1977  Physics::Transformations::Rotations::rotation_matrix_3d(axis,
1978  angle))
1979  {}
1980 
1981  Point<3>
1982  operator()(const Point<3> &p) const
1983  {
1984  return static_cast<Point<3>>(rotation_matrix * p);
1985  }
1986 
1987  private:
1989  };
1990 
1991 
1992  template <int spacedim>
1993  class Scale
1994  {
1995  public:
1996  explicit Scale(const double factor)
1997  : factor(factor)
1998  {}
2001  {
2002  return p * factor;
2003  }
2004 
2005  private:
2006  const double factor;
2007  };
2008  } // namespace internal
2009 
2010 
2011  template <int dim, int spacedim>
2012  void
2013  shift(const Tensor<1, spacedim> & shift_vector,
2015  {
2016  transform(internal::Shift<spacedim>(shift_vector), triangulation);
2017  }
2018 
2019 
2020  template <int dim>
2021  void
2023  const double angle,
2025  {
2026  transform(internal::Rotate3d(axis, angle), triangulation);
2027  }
2028 
2029 
2030  template <int dim>
2031  void
2032  rotate(const double angle,
2033  const unsigned int axis,
2035  {
2036  Assert(axis < 3, ExcMessage("Invalid axis given!"));
2037 
2038  Tensor<1, 3, double> vector;
2039  vector[axis] = 1.;
2040 
2041  transform(internal::Rotate3d(vector, angle), triangulation);
2042  }
2043 
2044 
2045  template <int dim, int spacedim>
2046  void
2047  scale(const double scaling_factor,
2049  {
2050  Assert(scaling_factor > 0, ExcScalingFactorNotPositive(scaling_factor));
2051  transform(internal::Scale<spacedim>(scaling_factor), triangulation);
2052  }
2053 
2054 
2055  namespace internal
2056  {
2062  inline void
2064  const AffineConstraints<double> &constraints,
2065  Vector<double> & u)
2066  {
2067  const unsigned int n_dofs = S.n();
2068  const auto op = linear_operator(S);
2069  const auto SF = constrained_linear_operator(constraints, op);
2071  prec.initialize(S, 1.2);
2072 
2073  SolverControl control(n_dofs, 1.e-10, false, false);
2075  SolverCG<Vector<double>> solver(control, mem);
2076 
2077  Vector<double> f(n_dofs);
2078 
2079  const auto constrained_rhs =
2080  constrained_right_hand_side(constraints, op, f);
2081  solver.solve(SF, u, constrained_rhs, prec);
2082 
2083  constraints.distribute(u);
2084  }
2085  } // namespace internal
2086 
2087 
2088  // Implementation for dimensions except 1
2089  template <int dim>
2090  void
2091  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
2093  const Function<dim> * coefficient,
2094  const bool solve_for_absolute_positions)
2095  {
2096  if (dim == 1)
2097  Assert(false, ExcNotImplemented());
2098 
2099  // first provide everything that is needed for solving a Laplace
2100  // equation.
2101  FE_Q<dim> q1(1);
2102 
2103  DoFHandler<dim> dof_handler(triangulation);
2104  dof_handler.distribute_dofs(q1);
2105 
2106  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
2107  DoFTools::make_sparsity_pattern(dof_handler, dsp);
2108  dsp.compress();
2109 
2110  SparsityPattern sparsity_pattern;
2111  sparsity_pattern.copy_from(dsp);
2112  sparsity_pattern.compress();
2113 
2114  SparseMatrix<double> S(sparsity_pattern);
2115 
2116  QGauss<dim> quadrature(4);
2117 
2119  StaticMappingQ1<dim>::mapping, dof_handler, quadrature, S, coefficient);
2120 
2121  // set up the boundary values for the laplace problem
2122  std::array<AffineConstraints<double>, dim> constraints;
2123  typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
2124  new_points.end();
2125 
2126  // fill these maps using the data given by new_points
2127  for (const auto &cell : dof_handler.active_cell_iterators())
2128  {
2129  // loop over all vertices of the cell and see if it is listed in the map
2130  // given as first argument of the function
2131  for (const unsigned int vertex_no : cell->vertex_indices())
2132  {
2133  const unsigned int vertex_index = cell->vertex_index(vertex_no);
2134  const Point<dim> & vertex_point = cell->vertex(vertex_no);
2135 
2136  const typename std::map<unsigned int, Point<dim>>::const_iterator
2137  map_iter = new_points.find(vertex_index);
2138 
2139  if (map_iter != map_end)
2140  for (unsigned int i = 0; i < dim; ++i)
2141  {
2142  constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
2143  constraints[i].set_inhomogeneity(
2144  cell->vertex_dof_index(vertex_no, 0),
2145  (solve_for_absolute_positions ?
2146  map_iter->second(i) :
2147  map_iter->second(i) - vertex_point[i]));
2148  }
2149  }
2150  }
2151 
2152  for (unsigned int i = 0; i < dim; ++i)
2153  constraints[i].close();
2154 
2155  // solve the dim problems with different right hand sides.
2156  Vector<double> us[dim];
2157  for (unsigned int i = 0; i < dim; ++i)
2158  us[i].reinit(dof_handler.n_dofs());
2159 
2160  // solve linear systems in parallel
2161  Threads::TaskGroup<> tasks;
2162  for (unsigned int i = 0; i < dim; ++i)
2163  tasks +=
2164  Threads::new_task(&internal::laplace_solve, S, constraints[i], us[i]);
2165  tasks.join_all();
2166 
2167  // change the coordinates of the points of the triangulation
2168  // according to the computed values
2169  std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2170  for (const auto &cell : dof_handler.active_cell_iterators())
2171  for (const unsigned int vertex_no : cell->vertex_indices())
2172  if (vertex_touched[cell->vertex_index(vertex_no)] == false)
2173  {
2174  Point<dim> &v = cell->vertex(vertex_no);
2175 
2176  const types::global_dof_index dof_index =
2177  cell->vertex_dof_index(vertex_no, 0);
2178  for (unsigned int i = 0; i < dim; ++i)
2179  if (solve_for_absolute_positions)
2180  v(i) = us[i](dof_index);
2181  else
2182  v(i) += us[i](dof_index);
2183 
2184  vertex_touched[cell->vertex_index(vertex_no)] = true;
2185  }
2186  }
2187 
2188  template <int dim, int spacedim>
2189  std::map<unsigned int, Point<spacedim>>
2191  {
2192  std::map<unsigned int, Point<spacedim>> vertex_map;
2194  cell = tria.begin_active(),
2195  endc = tria.end();
2196  for (; cell != endc; ++cell)
2197  {
2198  for (unsigned int i : cell->face_indices())
2199  {
2200  const typename Triangulation<dim, spacedim>::face_iterator &face =
2201  cell->face(i);
2202  if (face->at_boundary())
2203  {
2204  for (unsigned j = 0; j < face->n_vertices(); ++j)
2205  {
2206  const Point<spacedim> &vertex = face->vertex(j);
2207  const unsigned int vertex_index = face->vertex_index(j);
2208  vertex_map[vertex_index] = vertex;
2209  }
2210  }
2211  }
2212  }
2213  return vertex_map;
2214  }
2215 
2220  template <int dim, int spacedim>
2221  void
2222  distort_random(const double factor,
2224  const bool keep_boundary,
2225  const unsigned int seed)
2226  {
2227  // if spacedim>dim we need to make sure that we perturb
2228  // points but keep them on
2229  // the manifold. however, this isn't implemented right now
2230  Assert(spacedim == dim, ExcNotImplemented());
2231 
2232 
2233  // find the smallest length of the
2234  // lines adjacent to the
2235  // vertex. take the initial value
2236  // to be larger than anything that
2237  // might be found: the diameter of
2238  // the triangulation, here
2239  // estimated by adding up the
2240  // diameters of the coarse grid
2241  // cells.
2242  double almost_infinite_length = 0;
2243  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
2244  triangulation.begin(0);
2245  cell != triangulation.end(0);
2246  ++cell)
2247  almost_infinite_length += cell->diameter();
2248 
2249  std::vector<double> minimal_length(triangulation.n_vertices(),
2250  almost_infinite_length);
2251 
2252  // also note if a vertex is at the boundary
2253  std::vector<bool> at_boundary(keep_boundary ? triangulation.n_vertices() :
2254  0,
2255  false);
2256  // for parallel::shared::Triangulation we need to work on all vertices,
2257  // not just the ones related to locally owned cells;
2258  const bool is_parallel_shared =
2260  &triangulation) != nullptr);
2261  for (const auto &cell : triangulation.active_cell_iterators())
2262  if (is_parallel_shared || cell->is_locally_owned())
2263  {
2264  if (dim > 1)
2265  {
2266  for (unsigned int i = 0; i < cell->n_lines(); ++i)
2267  {
2269  line = cell->line(i);
2270 
2271  if (keep_boundary && line->at_boundary())
2272  {
2273  at_boundary[line->vertex_index(0)] = true;
2274  at_boundary[line->vertex_index(1)] = true;
2275  }
2276 
2277  minimal_length[line->vertex_index(0)] =
2278  std::min(line->diameter(),
2279  minimal_length[line->vertex_index(0)]);
2280  minimal_length[line->vertex_index(1)] =
2281  std::min(line->diameter(),
2282  minimal_length[line->vertex_index(1)]);
2283  }
2284  }
2285  else // dim==1
2286  {
2287  if (keep_boundary)
2288  for (unsigned int vertex = 0; vertex < 2; ++vertex)
2289  if (cell->at_boundary(vertex) == true)
2290  at_boundary[cell->vertex_index(vertex)] = true;
2291 
2292  minimal_length[cell->vertex_index(0)] =
2293  std::min(cell->diameter(),
2294  minimal_length[cell->vertex_index(0)]);
2295  minimal_length[cell->vertex_index(1)] =
2296  std::min(cell->diameter(),
2297  minimal_length[cell->vertex_index(1)]);
2298  }
2299  }
2300 
2301  // create a random number generator for the interval [-1,1]
2302  boost::random::mt19937 rng(seed);
2303  boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
2304 
2305  // If the triangulation is distributed, we need to
2306  // exchange the moved vertices across mpi processes
2307  if (auto distributed_triangulation =
2309  &triangulation))
2310  {
2311  const std::vector<bool> locally_owned_vertices =
2312  get_locally_owned_vertices(triangulation);
2313  std::vector<bool> vertex_moved(triangulation.n_vertices(), false);
2314 
2315  // Next move vertices on locally owned cells
2316  for (const auto &cell : triangulation.active_cell_iterators())
2317  if (cell->is_locally_owned())
2318  {
2319  for (const unsigned int vertex_no : cell->vertex_indices())
2320  {
2321  const unsigned global_vertex_no =
2322  cell->vertex_index(vertex_no);
2323 
2324  // ignore this vertex if we shall keep the boundary and
2325  // this vertex *is* at the boundary, if it is already moved
2326  // or if another process moves this vertex
2327  if ((keep_boundary && at_boundary[global_vertex_no]) ||
2328  vertex_moved[global_vertex_no] ||
2329  !locally_owned_vertices[global_vertex_no])
2330  continue;
2331 
2332  // first compute a random shift vector
2333  Point<spacedim> shift_vector;
2334  for (unsigned int d = 0; d < spacedim; ++d)
2335  shift_vector(d) = uniform_distribution(rng);
2336 
2337  shift_vector *= factor * minimal_length[global_vertex_no] /
2338  std::sqrt(shift_vector.square());
2339 
2340  // finally move the vertex
2341  cell->vertex(vertex_no) += shift_vector;
2342  vertex_moved[global_vertex_no] = true;
2343  }
2344  }
2345 
2346  distributed_triangulation->communicate_locally_moved_vertices(
2347  locally_owned_vertices);
2348  }
2349  else
2350  // if this is a sequential triangulation, we could in principle
2351  // use the algorithm above, but we'll use an algorithm that we used
2352  // before the parallel::distributed::Triangulation was introduced
2353  // in order to preserve backward compatibility
2354  {
2355  // loop over all vertices and compute their new locations
2356  const unsigned int n_vertices = triangulation.n_vertices();
2357  std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
2358  const std::vector<Point<spacedim>> &old_vertex_locations =
2359  triangulation.get_vertices();
2360 
2361  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2362  {
2363  // ignore this vertex if we will keep the boundary and
2364  // this vertex *is* at the boundary
2365  if (keep_boundary && at_boundary[vertex])
2366  new_vertex_locations[vertex] = old_vertex_locations[vertex];
2367  else
2368  {
2369  // compute a random shift vector
2370  Point<spacedim> shift_vector;
2371  for (unsigned int d = 0; d < spacedim; ++d)
2372  shift_vector(d) = uniform_distribution(rng);
2373 
2374  shift_vector *= factor * minimal_length[vertex] /
2375  std::sqrt(shift_vector.square());
2376 
2377  // record new vertex location
2378  new_vertex_locations[vertex] =
2379  old_vertex_locations[vertex] + shift_vector;
2380  }
2381  }
2382 
2383  // now do the actual move of the vertices
2384  for (const auto &cell : triangulation.active_cell_iterators())
2385  for (const unsigned int vertex_no : cell->vertex_indices())
2386  cell->vertex(vertex_no) =
2387  new_vertex_locations[cell->vertex_index(vertex_no)];
2388  }
2389 
2390  // Correct hanging nodes if necessary
2391  if (dim >= 2)
2392  {
2393  // We do the same as in GridTools::transform
2394  //
2395  // exclude hanging nodes at the boundaries of artificial cells:
2396  // these may belong to ghost cells for which we know the exact
2397  // location of vertices, whereas the artificial cell may or may
2398  // not be further refined, and so we cannot know whether
2399  // the location of the hanging node is correct or not
2401  cell = triangulation.begin_active(),
2402  endc = triangulation.end();
2403  for (; cell != endc; ++cell)
2404  if (!cell->is_artificial())
2405  for (const unsigned int face : cell->face_indices())
2406  if (cell->face(face)->has_children() &&
2407  !cell->face(face)->at_boundary())
2408  {
2409  // this face has hanging nodes
2410  if (dim == 2)
2411  cell->face(face)->child(0)->vertex(1) =
2412  (cell->face(face)->vertex(0) +
2413  cell->face(face)->vertex(1)) /
2414  2;
2415  else if (dim == 3)
2416  {
2417  cell->face(face)->child(0)->vertex(1) =
2418  .5 * (cell->face(face)->vertex(0) +
2419  cell->face(face)->vertex(1));
2420  cell->face(face)->child(0)->vertex(2) =
2421  .5 * (cell->face(face)->vertex(0) +
2422  cell->face(face)->vertex(2));
2423  cell->face(face)->child(1)->vertex(3) =
2424  .5 * (cell->face(face)->vertex(1) +
2425  cell->face(face)->vertex(3));
2426  cell->face(face)->child(2)->vertex(3) =
2427  .5 * (cell->face(face)->vertex(2) +
2428  cell->face(face)->vertex(3));
2429 
2430  // center of the face
2431  cell->face(face)->child(0)->vertex(3) =
2432  .25 * (cell->face(face)->vertex(0) +
2433  cell->face(face)->vertex(1) +
2434  cell->face(face)->vertex(2) +
2435  cell->face(face)->vertex(3));
2436  }
2437  }
2438  }
2439  }
2440 
2441 
2442 
2443  template <int dim, template <int, int> class MeshType, int spacedim>
2444  unsigned int
2445  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
2446  const Point<spacedim> & p,
2447  const std::vector<bool> & marked_vertices)
2448  {
2449  // first get the underlying triangulation from the mesh and determine
2450  // vertices and used vertices
2451  const Triangulation<dim, spacedim> &tria = mesh.get_triangulation();
2452 
2453  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
2454 
2455  Assert(tria.get_vertices().size() == marked_vertices.size() ||
2456  marked_vertices.size() == 0,
2457  ExcDimensionMismatch(tria.get_vertices().size(),
2458  marked_vertices.size()));
2459 
2460  // marked_vertices is expected to be a subset of used_vertices. Thus,
2461  // comparing the range marked_vertices.begin() to marked_vertices.end() with
2462  // the range used_vertices.begin() to used_vertices.end() the element in the
2463  // second range must be valid if the element in the first range is valid.
2464  Assert(
2465  marked_vertices.size() == 0 ||
2466  std::equal(marked_vertices.begin(),
2467  marked_vertices.end(),
2468  tria.get_used_vertices().begin(),
2469  [](bool p, bool q) { return !p || q; }),
2470  ExcMessage(
2471  "marked_vertices should be a subset of used vertices in the triangulation "
2472  "but marked_vertices contains one or more vertices that are not used vertices!"));
2473 
2474  // If marked_indices is empty, consider all used_vertices for finding the
2475  // closest vertex to the point. Otherwise, marked_indices is used.
2476  const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
2477  tria.get_used_vertices() :
2478  marked_vertices;
2479 
2480  // At the beginning, the first used vertex is considered to be the closest
2481  // one.
2482  std::vector<bool>::const_iterator first =
2483  std::find(vertices_to_use.begin(), vertices_to_use.end(), true);
2484 
2485  // Assert that at least one vertex is actually used
2486  Assert(first != vertices_to_use.end(), ExcInternalError());
2487 
2488  unsigned int best_vertex = std::distance(vertices_to_use.begin(), first);
2489  double best_dist = (p - vertices[best_vertex]).norm_square();
2490 
2491  // For all remaining vertices, test
2492  // whether they are any closer
2493  for (unsigned int j = best_vertex + 1; j < vertices.size(); ++j)
2494  if (vertices_to_use[j])
2495  {
2496  const double dist = (p - vertices[j]).norm_square();
2497  if (dist < best_dist)
2498  {
2499  best_vertex = j;
2500  best_dist = dist;
2501  }
2502  }
2503 
2504  return best_vertex;
2505  }
2506 
2507 
2508 
2509  template <int dim, template <int, int> class MeshType, int spacedim>
2510  unsigned int
2512  const MeshType<dim, spacedim> &mesh,
2513  const Point<spacedim> & p,
2514  const std::vector<bool> & marked_vertices)
2515  {
2516  // Take a shortcut in the simple case.
2517  if (mapping.preserves_vertex_locations() == true)
2518  return find_closest_vertex(mesh, p, marked_vertices);
2519 
2520  // first get the underlying triangulation from the mesh and determine
2521  // vertices and used vertices
2522  const Triangulation<dim, spacedim> &tria = mesh.get_triangulation();
2523 
2524  auto vertices = extract_used_vertices(tria, mapping);
2525 
2526  Assert(tria.get_vertices().size() == marked_vertices.size() ||
2527  marked_vertices.size() == 0,
2528  ExcDimensionMismatch(tria.get_vertices().size(),
2529  marked_vertices.size()));
2530 
2531  // marked_vertices is expected to be a subset of used_vertices. Thus,
2532  // comparing the range marked_vertices.begin() to marked_vertices.end()
2533  // with the range used_vertices.begin() to used_vertices.end() the element
2534  // in the second range must be valid if the element in the first range is
2535  // valid.
2536  Assert(
2537  marked_vertices.size() == 0 ||
2538  std::equal(marked_vertices.begin(),
2539  marked_vertices.end(),
2540  tria.get_used_vertices().begin(),
2541  [](bool p, bool q) { return !p || q; }),
2542  ExcMessage(
2543  "marked_vertices should be a subset of used vertices in the triangulation "
2544  "but marked_vertices contains one or more vertices that are not used vertices!"));
2545 
2546  // Remove from the map unwanted elements.
2547  if (marked_vertices.size() != 0)
2548  for (auto it = vertices.begin(); it != vertices.end();)
2549  {
2550  if (marked_vertices[it->first] == false)
2551  {
2552  it = vertices.erase(it);
2553  }
2554  else
2555  {
2556  ++it;
2557  }
2558  }
2559 
2560  return find_closest_vertex(vertices, p);
2561  }
2562 
2563 
2564 
2565  template <int dim, template <int, int> class MeshType, int spacedim>
2566 #ifndef _MSC_VER
2567  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
2568 #else
2569  std::vector<
2570  typename ::internal::
2571  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2572 #endif
2573  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &mesh,
2574  const unsigned int vertex)
2575  {
2576  // make sure that the given vertex is
2577  // an active vertex of the underlying
2578  // triangulation
2579  AssertIndexRange(vertex, mesh.get_triangulation().n_vertices());
2580  Assert(mesh.get_triangulation().get_used_vertices()[vertex],
2581  ExcVertexNotUsed(vertex));
2582 
2583  // use a set instead of a vector
2584  // to ensure that cells are inserted only
2585  // once
2586  std::set<typename ::internal::
2587  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2589 
2590  // go through all active cells and look if the vertex is part of that cell
2591  //
2592  // in 1d, this is all we need to care about. in 2d/3d we also need to worry
2593  // that the vertex might be a hanging node on a face or edge of a cell; in
2594  // this case, we would want to add those cells as well on whose faces the
2595  // vertex is located but for which it is not a vertex itself.
2596  //
2597  // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging
2598  // node can only be in the middle of a face and we can query the neighboring
2599  // cell from the current cell. on the other hand, in 3d a hanging node
2600  // vertex can also be on an edge but there can be many other cells on
2601  // this edge and we can not access them from the cell we are currently
2602  // on.
2603  //
2604  // so, in the 3d case, if we run the algorithm as in 2d, we catch all
2605  // those cells for which the vertex we seek is on a *subface*, but we
2606  // miss the case of cells for which the vertex we seek is on a
2607  // sub-edge for which there is no corresponding sub-face (because the
2608  // immediate neighbor behind this face is not refined), see for example
2609  // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we
2610  // haven't yet found the vertex for the current cell we also need to
2611  // look at the mid-points of edges
2612  //
2613  // as a final note, deciding whether a neighbor is actually coarser is
2614  // simple in the case of isotropic refinement (we just need to look at
2615  // the level of the current and the neighboring cell). however, this
2616  // isn't so simple if we have used anisotropic refinement since then
2617  // the level of a cell is not indicative of whether it is coarser or
2618  // not than the current cell. ultimately, we want to add all cells on
2619  // which the vertex is, independent of whether they are coarser or
2620  // finer and so in the 2d case below we simply add *any* *active* neighbor.
2621  // in the worst case, we add cells multiple times to the adjacent_cells
2622  // list, but std::set throws out those cells already entered
2623  for (const auto &cell : mesh.active_cell_iterators())
2624  {
2625  for (const unsigned int v : cell->vertex_indices())
2626  if (cell->vertex_index(v) == vertex)
2627  {
2628  // OK, we found a cell that contains
2629  // the given vertex. We add it
2630  // to the list.
2631  adjacent_cells.insert(cell);
2632 
2633  // as explained above, in 2+d we need to check whether
2634  // this vertex is on a face behind which there is a
2635  // (possibly) coarser neighbor. if this is the case,
2636  // then we need to also add this neighbor
2637  if (dim >= 2)
2638  for (const auto face :
2639  cell->reference_cell().faces_for_given_vertex(v))
2640  if (!cell->at_boundary(face) &&
2641  cell->neighbor(face)->is_active())
2642  {
2643  // there is a (possibly) coarser cell behind a
2644  // face to which the vertex belongs. the
2645  // vertex we are looking at is then either a
2646  // vertex of that coarser neighbor, or it is a
2647  // hanging node on one of the faces of that
2648  // cell. in either case, it is adjacent to the
2649  // vertex, so add it to the list as well (if
2650  // the cell was already in the list then the
2651  // std::set makes sure that we get it only
2652  // once)
2653  adjacent_cells.insert(cell->neighbor(face));
2654  }
2655 
2656  // in any case, we have found a cell, so go to the next cell
2657  goto next_cell;
2658  }
2659 
2660  // in 3d also loop over the edges
2661  if (dim >= 3)
2662  {
2663  for (unsigned int e = 0; e < cell->n_lines(); ++e)
2664  if (cell->line(e)->has_children())
2665  // the only place where this vertex could have been
2666  // hiding is on the mid-edge point of the edge we
2667  // are looking at
2668  if (cell->line(e)->child(0)->vertex_index(1) == vertex)
2669  {
2670  adjacent_cells.insert(cell);
2671 
2672  // jump out of this tangle of nested loops
2673  goto next_cell;
2674  }
2675  }
2676 
2677  // in more than 3d we would probably have to do the same as
2678  // above also for even lower-dimensional objects
2679  Assert(dim <= 3, ExcNotImplemented());
2680 
2681  // move on to the next cell if we have found the
2682  // vertex on the current one
2683  next_cell:;
2684  }
2685 
2686  // if this was an active vertex then there needs to have been
2687  // at least one cell to which it is adjacent!
2688  Assert(adjacent_cells.size() > 0, ExcInternalError());
2689 
2690  // return the result as a vector, rather than the set we built above
2691  return std::vector<
2692  typename ::internal::
2693  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
2694  adjacent_cells.begin(), adjacent_cells.end());
2695  }
2696 
2697 
2698 
2699  template <int dim, int spacedim>
2700  std::vector<std::vector<Tensor<1, spacedim>>>
2702  const Triangulation<dim, spacedim> &mesh,
2703  const std::vector<
2705  &vertex_to_cells)
2706  {
2707  const std::vector<Point<spacedim>> &vertices = mesh.get_vertices();
2708  const unsigned int n_vertices = vertex_to_cells.size();
2709 
2710  AssertDimension(vertices.size(), n_vertices);
2711 
2712 
2713  std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
2714  n_vertices);
2715  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2716  if (mesh.vertex_used(vertex))
2717  {
2718  const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
2719  vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
2720 
2721  typename std::set<typename Triangulation<dim, spacedim>::
2722  active_cell_iterator>::iterator it =
2723  vertex_to_cells[vertex].begin();
2724  for (unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
2725  {
2726  vertex_to_cell_centers[vertex][cell] =
2727  (*it)->center() - vertices[vertex];
2728  vertex_to_cell_centers[vertex][cell] /=
2729  vertex_to_cell_centers[vertex][cell].norm();
2730  }
2731  }
2732  return vertex_to_cell_centers;
2733  }
2734 
2735 
2736  namespace internal
2737  {
2738  template <int spacedim>
2739  bool
2741  const unsigned int a,
2742  const unsigned int b,
2743  const Tensor<1, spacedim> & point_direction,
2744  const std::vector<Tensor<1, spacedim>> &center_directions)
2745  {
2746  const double scalar_product_a = center_directions[a] * point_direction;
2747  const double scalar_product_b = center_directions[b] * point_direction;
2748 
2749  // The function is supposed to return if a is before b. We are looking
2750  // for the alignment of point direction and center direction, therefore
2751  // return if the scalar product of a is larger.
2752  return (scalar_product_a > scalar_product_b);
2753  }
2754  } // namespace internal
2755 
2756  template <int dim, template <int, int> class MeshType, int spacedim>
2757 #ifndef _MSC_VER
2758  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
2759 #else
2760  std::pair<typename ::internal::
2761  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
2762  Point<dim>>
2763 #endif
2765  const Mapping<dim, spacedim> & mapping,
2766  const MeshType<dim, spacedim> &mesh,
2767  const Point<spacedim> & p,
2768  const std::vector<
2769  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
2770  & vertex_to_cells,
2771  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
2772  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
2773  const std::vector<bool> & marked_vertices,
2774  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree,
2775  const double tolerance,
2776  const RTree<
2777  std::pair<BoundingBox<spacedim>,
2779  *relevant_cell_bounding_boxes_rtree)
2780  {
2781  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2782  Point<dim>>
2783  cell_and_position;
2784  cell_and_position.first = mesh.end();
2785 
2786  // To handle points at the border we keep track of points which are close to
2787  // the unit cell:
2788  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2789  Point<dim>>
2790  cell_and_position_approx;
2791 
2792  if (relevant_cell_bounding_boxes_rtree != nullptr &&
2793  !relevant_cell_bounding_boxes_rtree->empty())
2794  {
2795  // create a bounding box around point p with 2*tolerance as side length.
2796  auto p1 = p;
2797  auto p2 = p;
2798 
2799  for (int d = 0; d < spacedim; ++d)
2800  {
2801  p1[d] = p1[d] - tolerance;
2802  p2[d] = p2[d] + tolerance;
2803  }
2804 
2805  BoundingBox<spacedim> bb({p1, p2});
2806 
2807  if (relevant_cell_bounding_boxes_rtree->qbegin(
2808  boost::geometry::index::intersects(bb)) ==
2809  relevant_cell_bounding_boxes_rtree->qend())
2810  return cell_and_position;
2811  }
2812 
2813  bool found_cell = false;
2814  bool approx_cell = false;
2815 
2816  unsigned int closest_vertex_index = 0;
2817  Tensor<1, spacedim> vertex_to_point;
2818  auto current_cell = cell_hint;
2819 
2820  while (found_cell == false)
2821  {
2822  // First look at the vertices of the cell cell_hint. If it's an
2823  // invalid cell, then query for the closest global vertex
2824  if (current_cell.state() == IteratorState::valid)
2825  {
2826  const auto cell_vertices = mapping.get_vertices(current_cell);
2827  const unsigned int closest_vertex =
2828  find_closest_vertex_of_cell<dim, spacedim>(current_cell,
2829  p,
2830  mapping);
2831  vertex_to_point = p - cell_vertices[closest_vertex];
2832  closest_vertex_index = current_cell->vertex_index(closest_vertex);
2833  }
2834  else
2835  {
2836  if (!used_vertices_rtree.empty())
2837  {
2838  // If we have an rtree at our disposal, use it.
2839  using ValueType = std::pair<Point<spacedim>, unsigned int>;
2840  std::function<bool(const ValueType &)> marked;
2841  if (marked_vertices.size() == mesh.n_vertices())
2842  marked = [&marked_vertices](const ValueType &value) -> bool {
2843  return marked_vertices[value.second];
2844  };
2845  else
2846  marked = [](const ValueType &) -> bool { return true; };
2847 
2848  std::vector<std::pair<Point<spacedim>, unsigned int>> res;
2849  used_vertices_rtree.query(
2850  boost::geometry::index::nearest(p, 1) &&
2851  boost::geometry::index::satisfies(marked),
2852  std::back_inserter(res));
2853 
2854  // We should have one and only one result
2855  AssertDimension(res.size(), 1);
2856  closest_vertex_index = res[0].second;
2857  }
2858  else
2859  {
2860  closest_vertex_index = GridTools::find_closest_vertex(
2861  mapping, mesh, p, marked_vertices);
2862  }
2863  vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
2864  }
2865 
2866  const double vertex_point_norm = vertex_to_point.norm();
2867  if (vertex_point_norm > 0)
2868  vertex_to_point /= vertex_point_norm;
2869 
2870  const unsigned int n_neighbor_cells =
2871  vertex_to_cells[closest_vertex_index].size();
2872 
2873  // Create a corresponding map of vectors from vertex to cell center
2874  std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
2875 
2876  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
2877  neighbor_permutation[i] = i;
2878 
2879  auto comp = [&](const unsigned int a, const unsigned int b) -> bool {
2880  return internal::compare_point_association<spacedim>(
2881  a,
2882  b,
2883  vertex_to_point,
2884  vertex_to_cell_centers[closest_vertex_index]);
2885  };
2886 
2887  std::sort(neighbor_permutation.begin(),
2888  neighbor_permutation.end(),
2889  comp);
2890  // It is possible the vertex is close
2891  // to an edge, thus we add a tolerance
2892  // to keep also the "best" cell
2893  double best_distance = tolerance;
2894 
2895  // Search all of the cells adjacent to the closest vertex of the cell
2896  // hint Most likely we will find the point in them.
2897  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
2898  {
2899  try
2900  {
2901  auto cell = vertex_to_cells[closest_vertex_index].begin();
2902  std::advance(cell, neighbor_permutation[i]);
2903 
2904  if (!(*cell)->is_artificial())
2905  {
2906  const Point<dim> p_unit =
2907  mapping.transform_real_to_unit_cell(*cell, p);
2909  tolerance))
2910  {
2911  cell_and_position.first = *cell;
2912  cell_and_position.second = p_unit;
2913  found_cell = true;
2914  approx_cell = false;
2915  break;
2916  }
2917  // The point is not inside this cell: checking how far
2918  // outside it is and whether we want to use this cell as a
2919  // backup if we can't find a cell within which the point
2920  // lies.
2921  const double dist =
2923  if (dist < best_distance)
2924  {
2925  best_distance = dist;
2926  cell_and_position_approx.first = *cell;
2927  cell_and_position_approx.second = p_unit;
2928  approx_cell = true;
2929  }
2930  }
2931  }
2932  catch (typename Mapping<dim>::ExcTransformationFailed &)
2933  {}
2934  }
2935 
2936  if (found_cell == true)
2937  return cell_and_position;
2938  else if (approx_cell == true)
2939  return cell_and_position_approx;
2940 
2941  // The first time around, we check for vertices in the hint_cell. If
2942  // that does not work, we set the cell iterator to an invalid one, and
2943  // look for a global vertex close to the point. If that does not work,
2944  // we are in trouble, and just throw an exception.
2945  //
2946  // If we got here, then we did not find the point. If the
2947  // current_cell.state() here is not IteratorState::valid, it means that
2948  // the user did not provide a hint_cell, and at the beginning of the
2949  // while loop we performed an actual global search on the mesh
2950  // vertices. Not finding the point then means the point is outside the
2951  // domain, or that we've had problems with the algorithm above. Try as a
2952  // last resort the other (simpler) algorithm.
2953  if (current_cell.state() != IteratorState::valid)
2955  mapping, mesh, p, marked_vertices, tolerance);
2956 
2957  current_cell = typename MeshType<dim, spacedim>::active_cell_iterator();
2958  }
2959  return cell_and_position;
2960  }
2961 
2962 
2963 
2964  template <int dim, int spacedim>
2965  unsigned int
2968  const Point<spacedim> & position,
2969  const Mapping<dim, spacedim> & mapping)
2970  {
2971  const auto vertices = mapping.get_vertices(cell);
2972  double minimum_distance = position.distance_square(vertices[0]);
2973  unsigned int closest_vertex = 0;
2974 
2975  for (unsigned int v = 1; v < cell->n_vertices(); ++v)
2976  {
2977  const double vertex_distance = position.distance_square(vertices[v]);
2978  if (vertex_distance < minimum_distance)
2979  {
2980  closest_vertex = v;
2981  minimum_distance = vertex_distance;
2982  }
2983  }
2984  return closest_vertex;
2985  }
2986 
2987 
2988 
2989  namespace internal
2990  {
2991  namespace BoundingBoxPredicate
2992  {
2993  template <class MeshType>
2994  std::tuple<BoundingBox<MeshType::space_dimension>, bool>
2996  const typename MeshType::cell_iterator &parent_cell,
2997  const std::function<
2998  bool(const typename MeshType::active_cell_iterator &)> &predicate)
2999  {
3000  bool has_predicate =
3001  false; // Start assuming there's no cells with predicate inside
3002  std::vector<typename MeshType::active_cell_iterator> active_cells;
3003  if (parent_cell->is_active())
3004  active_cells = {parent_cell};
3005  else
3006  // Finding all active cells descendants of the current one (or the
3007  // current one if it is active)
3008  active_cells = get_active_child_cells<MeshType>(parent_cell);
3009 
3010  const unsigned int spacedim = MeshType::space_dimension;
3011 
3012  // Looking for the first active cell which has the property predicate
3013  unsigned int i = 0;
3014  while (i < active_cells.size() && !predicate(active_cells[i]))
3015  ++i;
3016 
3017  // No active cells or no active cells with property
3018  if (active_cells.size() == 0 || i == active_cells.size())
3019  {
3020  BoundingBox<spacedim> bbox;
3021  return std::make_tuple(bbox, has_predicate);
3022  }
3023 
3024  // The two boundary points defining the boundary box
3025  Point<spacedim> maxp = active_cells[i]->vertex(0);
3026  Point<spacedim> minp = active_cells[i]->vertex(0);
3027 
3028  for (; i < active_cells.size(); ++i)
3029  if (predicate(active_cells[i]))
3030  for (const unsigned int v : active_cells[i]->vertex_indices())
3031  for (unsigned int d = 0; d < spacedim; ++d)
3032  {
3033  minp[d] = std::min(minp[d], active_cells[i]->vertex(v)[d]);
3034  maxp[d] = std::max(maxp[d], active_cells[i]->vertex(v)[d]);
3035  }
3036 
3037  has_predicate = true;
3038  BoundingBox<spacedim> bbox(std::make_pair(minp, maxp));
3039  return std::make_tuple(bbox, has_predicate);
3040  }
3041  } // namespace BoundingBoxPredicate
3042  } // namespace internal
3043 
3044 
3045 
3046  template <class MeshType>
3047  std::vector<BoundingBox<MeshType::space_dimension>>
3049  const MeshType &mesh,
3050  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3051  & predicate,
3052  const unsigned int refinement_level,
3053  const bool allow_merge,
3054  const unsigned int max_boxes)
3055  {
3056  // Algorithm brief description: begin with creating bounding boxes of all
3057  // cells at refinement_level (and coarser levels if there are active cells)
3058  // which have the predicate property. These are then merged
3059 
3060  Assert(
3061  refinement_level <= mesh.n_levels(),
3062  ExcMessage(
3063  "Error: refinement level is higher then total levels in the triangulation!"));
3064 
3065  const unsigned int spacedim = MeshType::space_dimension;
3066  std::vector<BoundingBox<spacedim>> bounding_boxes;
3067 
3068  // Creating a bounding box for all active cell on coarser level
3069 
3070  for (unsigned int i = 0; i < refinement_level; ++i)
3071  for (const typename MeshType::cell_iterator &cell :
3072  mesh.active_cell_iterators_on_level(i))
3073  {
3074  bool has_predicate = false;
3075  BoundingBox<spacedim> bbox;
3076  std::tie(bbox, has_predicate) =
3078  MeshType>(cell, predicate);
3079  if (has_predicate)
3080  bounding_boxes.push_back(bbox);
3081  }
3082 
3083  // Creating a Bounding Box for all cells on the chosen refinement_level
3084  for (const typename MeshType::cell_iterator &cell :
3085  mesh.cell_iterators_on_level(refinement_level))
3086  {
3087  bool has_predicate = false;
3088  BoundingBox<spacedim> bbox;
3089  std::tie(bbox, has_predicate) =
3091  MeshType>(cell, predicate);
3092  if (has_predicate)
3093  bounding_boxes.push_back(bbox);
3094  }
3095 
3096  if (!allow_merge)
3097  // If merging is not requested return the created bounding_boxes
3098  return bounding_boxes;
3099  else
3100  {
3101  // Merging part of the algorithm
3102  // Part 1: merging neighbors
3103  // This array stores the indices of arrays we have already merged
3104  std::vector<unsigned int> merged_boxes_idx;
3105  bool found_neighbors = true;
3106 
3107  // We merge only neighbors which can be expressed by a single bounding
3108  // box e.g. in 1d [0,1] and [1,2] can be described with [0,2] without
3109  // losing anything
3110  while (found_neighbors)
3111  {
3112  found_neighbors = false;
3113  for (unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
3114  {
3115  if (std::find(merged_boxes_idx.begin(),
3116  merged_boxes_idx.end(),
3117  i) == merged_boxes_idx.end())
3118  for (unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
3119  if (std::find(merged_boxes_idx.begin(),
3120  merged_boxes_idx.end(),
3121  j) == merged_boxes_idx.end() &&
3122  bounding_boxes[i].get_neighbor_type(
3123  bounding_boxes[j]) ==
3125  {
3126  bounding_boxes[i].merge_with(bounding_boxes[j]);
3127  merged_boxes_idx.push_back(j);
3128  found_neighbors = true;
3129  }
3130  }
3131  }
3132 
3133  // Copying the merged boxes into merged_b_boxes
3134  std::vector<BoundingBox<spacedim>> merged_b_boxes;
3135  for (unsigned int i = 0; i < bounding_boxes.size(); ++i)
3136  if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
3137  merged_boxes_idx.end())
3138  merged_b_boxes.push_back(bounding_boxes[i]);
3139 
3140  // Part 2: if there are too many bounding boxes, merging smaller boxes
3141  // This has sense only in dimension 2 or greater, since in dimension 1,
3142  // neighboring intervals can always be merged without problems
3143  if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
3144  {
3145  std::vector<double> volumes;
3146  for (unsigned int i = 0; i < merged_b_boxes.size(); ++i)
3147  volumes.push_back(merged_b_boxes[i].volume());
3148 
3149  while (merged_b_boxes.size() > max_boxes)
3150  {
3151  unsigned int min_idx =
3152  std::min_element(volumes.begin(), volumes.end()) -
3153  volumes.begin();
3154  volumes.erase(volumes.begin() + min_idx);
3155  // Finding a neighbor
3156  bool not_removed = true;
3157  for (unsigned int i = 0;
3158  i < merged_b_boxes.size() && not_removed;
3159  ++i)
3160  // We merge boxes if we have "attached" or "mergeable"
3161  // neighbors, even though mergeable should be dealt with in
3162  // Part 1
3163  if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
3164  merged_b_boxes[min_idx]) ==
3166  merged_b_boxes[i].get_neighbor_type(
3167  merged_b_boxes[min_idx]) ==
3169  {
3170  merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
3171  merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
3172  not_removed = false;
3173  }
3174  Assert(!not_removed,
3175  ExcMessage("Error: couldn't merge bounding boxes!"));
3176  }
3177  }
3178  Assert(merged_b_boxes.size() <= max_boxes,
3179  ExcMessage(
3180  "Error: couldn't reach target number of bounding boxes!"));
3181  return merged_b_boxes;
3182  }
3183  }
3184 
3185 
3186 
3187  template <int spacedim>
3188 #ifndef DOXYGEN
3189  std::tuple<std::vector<std::vector<unsigned int>>,
3190  std::map<unsigned int, unsigned int>,
3191  std::map<unsigned int, std::vector<unsigned int>>>
3192 #else
3193  return_type
3194 #endif
3196  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3197  const std::vector<Point<spacedim>> & points)
3198  {
3199  unsigned int n_procs = global_bboxes.size();
3200  std::vector<std::vector<unsigned int>> point_owners(n_procs);
3201  std::map<unsigned int, unsigned int> map_owners_found;
3202  std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3203 
3204  unsigned int n_points = points.size();
3205  for (unsigned int pt = 0; pt < n_points; ++pt)
3206  {
3207  // Keep track of how many processes we guess to own the point
3208  std::vector<unsigned int> owners_found;
3209  // Check in which other processes the point might be
3210  for (unsigned int rk = 0; rk < n_procs; ++rk)
3211  {
3212  for (const BoundingBox<spacedim> &bbox : global_bboxes[rk])
3213  if (bbox.point_inside(points[pt]))
3214  {
3215  point_owners[rk].emplace_back(pt);
3216  owners_found.emplace_back(rk);
3217  break; // We can check now the next process
3218  }
3219  }
3220  Assert(owners_found.size() > 0,
3221  ExcMessage("No owners found for the point " +
3222  std::to_string(pt)));
3223  if (owners_found.size() == 1)
3224  map_owners_found[pt] = owners_found[0];
3225  else
3226  // Multiple owners
3227  map_owners_guessed[pt] = owners_found;
3228  }
3229 
3230  return std::make_tuple(std::move(point_owners),
3231  std::move(map_owners_found),
3232  std::move(map_owners_guessed));
3233  }
3234 
3235  template <int spacedim>
3236 #ifndef DOXYGEN
3237  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
3238  std::map<unsigned int, unsigned int>,
3239  std::map<unsigned int, std::vector<unsigned int>>>
3240 #else
3241  return_type
3242 #endif
3244  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
3245  const std::vector<Point<spacedim>> & points)
3246  {
3247  std::map<unsigned int, std::vector<unsigned int>> point_owners;
3248  std::map<unsigned int, unsigned int> map_owners_found;
3249  std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3250  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> search_result;
3251 
3252  unsigned int n_points = points.size();
3253  for (unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
3254  {
3255  search_result.clear(); // clearing last output
3256 
3257  // Running tree search
3258  covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
3259  std::back_inserter(search_result));
3260 
3261  // Keep track of how many processes we guess to own the point
3262  std::set<unsigned int> owners_found;
3263  // Check in which other processes the point might be
3264  for (const auto &rank_bbox : search_result)
3265  {
3266  // Try to add the owner to the owners found,
3267  // and check if it was already present
3268  const bool pt_inserted = owners_found.insert(pt_n).second;
3269  if (pt_inserted)
3270  point_owners[rank_bbox.second].emplace_back(pt_n);
3271  }
3272  Assert(owners_found.size() > 0,
3273  ExcMessage("No owners found for the point " +
3274  std::to_string(pt_n)));
3275  if (owners_found.size() == 1)
3276  map_owners_found[pt_n] = *owners_found.begin();
3277  else
3278  // Multiple owners
3279  std::copy(owners_found.begin(),
3280  owners_found.end(),
3281  std::back_inserter(map_owners_guessed[pt_n]));
3282  }
3283 
3284  return std::make_tuple(std::move(point_owners),
3285  std::move(map_owners_found),
3286  std::move(map_owners_guessed));
3287  }
3288 
3289 
3290  template <int dim, int spacedim>
3291  std::vector<
3292  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3294  {
3295  std::vector<
3296  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3297  vertex_to_cell_map(triangulation.n_vertices());
3299  cell = triangulation.begin_active(),
3300  endc = triangulation.end();
3301  for (; cell != endc; ++cell)
3302  for (const unsigned int i : cell->vertex_indices())
3303  vertex_to_cell_map[cell->vertex_index(i)].insert(cell);
3304 
3305  // Take care of hanging nodes
3306  cell = triangulation.begin_active();
3307  for (; cell != endc; ++cell)
3308  {
3309  for (unsigned int i : cell->face_indices())
3310  {
3311  if ((cell->at_boundary(i) == false) &&
3312  (cell->neighbor(i)->is_active()))
3313  {
3315  adjacent_cell = cell->neighbor(i);
3316  for (unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
3317  vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert(
3318  adjacent_cell);
3319  }
3320  }
3321 
3322  // in 3d also loop over the edges
3323  if (dim == 3)
3324  {
3325  for (unsigned int i = 0; i < cell->n_lines(); ++i)
3326  if (cell->line(i)->has_children())
3327  // the only place where this vertex could have been
3328  // hiding is on the mid-edge point of the edge we
3329  // are looking at
3330  vertex_to_cell_map[cell->line(i)->child(0)->vertex_index(1)]
3331  .insert(cell);
3332  }
3333  }
3334 
3335  return vertex_to_cell_map;
3336  }
3337 
3338 
3339 
3340  template <int dim, int spacedim>
3341  std::map<unsigned int, types::global_vertex_index>
3344  {
3345  std::map<unsigned int, types::global_vertex_index>
3346  local_to_global_vertex_index;
3347 
3348 #ifndef DEAL_II_WITH_MPI
3349 
3350  // without MPI, this function doesn't make sense because on cannot
3351  // use parallel::distributed::Triangulation in any meaningful
3352  // way
3353  (void)triangulation;
3354  Assert(false,
3355  ExcMessage("This function does not make any sense "
3356  "for parallel::distributed::Triangulation "
3357  "objects if you do not have MPI enabled."));
3358 
3359 #else
3360 
3361  using active_cell_iterator =
3363  const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
3364  vertex_to_cell_map(triangulation);
3365 
3366  // Create a local index for the locally "owned" vertices
3367  types::global_vertex_index next_index = 0;
3368  unsigned int max_cellid_size = 0;
3369  std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
3370  vertices_added;
3371  std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
3372  std::map<types::subdomain_id,
3373  std::vector<std::tuple<types::global_vertex_index,
3375  std::string>>>
3376  vertices_to_send;
3377  active_cell_iterator cell = triangulation.begin_active(),
3378  endc = triangulation.end();
3379  std::set<active_cell_iterator> missing_vert_cells;
3380  std::set<unsigned int> used_vertex_index;
3381  for (; cell != endc; ++cell)
3382  {
3383  if (cell->is_locally_owned())
3384  {
3385  for (const unsigned int i : cell->vertex_indices())
3386  {
3387  types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
3388  typename std::set<active_cell_iterator>::iterator
3389  adjacent_cell = vertex_to_cell[cell->vertex_index(i)].begin(),
3390  end_adj_cell = vertex_to_cell[cell->vertex_index(i)].end();
3391  for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
3392  lowest_subdomain_id =
3393  std::min(lowest_subdomain_id,
3394  (*adjacent_cell)->subdomain_id());
3395 
3396  // See if I "own" this vertex
3397  if (lowest_subdomain_id == cell->subdomain_id())
3398  {
3399  // Check that the vertex we are working on a vertex that has
3400  // not be dealt with yet
3401  if (used_vertex_index.find(cell->vertex_index(i)) ==
3402  used_vertex_index.end())
3403  {
3404  // Set the local index
3405  local_to_global_vertex_index[cell->vertex_index(i)] =
3406  next_index++;
3407 
3408  // Store the information that will be sent to the
3409  // adjacent cells on other subdomains
3410  adjacent_cell =
3411  vertex_to_cell[cell->vertex_index(i)].begin();
3412  for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
3413  if ((*adjacent_cell)->subdomain_id() !=
3414  cell->subdomain_id())
3415  {
3416  std::pair<types::subdomain_id,
3417  types::global_vertex_index>
3418  tmp((*adjacent_cell)->subdomain_id(),
3419  cell->vertex_index(i));
3420  if (vertices_added.find(tmp) ==
3421  vertices_added.end())
3422  {
3423  vertices_to_send[(*adjacent_cell)
3424  ->subdomain_id()]
3425  .emplace_back(i,
3426  cell->vertex_index(i),
3427  cell->id().to_string());
3428  if (cell->id().to_string().size() >
3429  max_cellid_size)
3430  max_cellid_size =
3431  cell->id().to_string().size();
3432  vertices_added.insert(tmp);
3433  }
3434  }
3435  used_vertex_index.insert(cell->vertex_index(i));
3436  }
3437  }
3438  else
3439  {
3440  // We don't own the vertex so we will receive its global
3441  // index
3442  vertices_to_recv[lowest_subdomain_id].insert(
3443  cell->vertex_index(i));
3444  missing_vert_cells.insert(cell);
3445  }
3446  }
3447  }
3448 
3449  // Some hanging nodes are vertices of ghost cells. They need to be
3450  // received.
3451  if (cell->is_ghost())
3452  {
3453  for (unsigned int i : cell->face_indices())
3454  {
3455  if (cell->at_boundary(i) == false)
3456  {
3457  if (cell->neighbor(i)->is_active())
3458  {
3459  typename Triangulation<dim,
3460  spacedim>::active_cell_iterator
3461  adjacent_cell = cell->neighbor(i);
3462  if ((adjacent_cell->is_locally_owned()))
3463  {
3464  types::subdomain_id adj_subdomain_id =
3465  adjacent_cell->subdomain_id();
3466  if (cell->subdomain_id() < adj_subdomain_id)
3467  for (unsigned int j = 0;
3468  j < cell->face(i)->n_vertices();
3469  ++j)
3470  {
3471  vertices_to_recv[cell->subdomain_id()].insert(
3472  cell->face(i)->vertex_index(j));
3473  missing_vert_cells.insert(cell);
3474  }
3475  }
3476  }
3477  }
3478  }
3479  }
3480  }
3481 
3482  // Get the size of the largest CellID string
3483  max_cellid_size =
3484  Utilities::MPI::max(max_cellid_size, triangulation.get_communicator());
3485 
3486  // Make indices global by getting the number of vertices owned by each
3487  // processors and shifting the indices accordingly
3488  types::global_vertex_index shift = 0;
3489  int ierr = MPI_Exscan(&next_index,
3490  &shift,
3491  1,
3493  MPI_SUM,
3494  triangulation.get_communicator());
3495  AssertThrowMPI(ierr);
3496 
3497  std::map<unsigned int, types::global_vertex_index>::iterator
3498  global_index_it = local_to_global_vertex_index.begin(),
3499  global_index_end = local_to_global_vertex_index.end();
3500  for (; global_index_it != global_index_end; ++global_index_it)
3501  global_index_it->second += shift;
3502 
3503 
3504  const int mpi_tag = Utilities::MPI::internal::Tags::
3506  const int mpi_tag2 = Utilities::MPI::internal::Tags::
3508 
3509 
3510  // In a first message, send the global ID of the vertices and the local
3511  // positions in the cells. In a second messages, send the cell ID as a
3512  // resize string. This is done in two messages so that types are not mixed
3513 
3514  // Send the first message
3515  std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
3516  vertices_to_send.size());
3517  std::vector<MPI_Request> first_requests(vertices_to_send.size());
3518  typename std::map<types::subdomain_id,
3519  std::vector<std::tuple<types::global_vertex_index,
3521  std::string>>>::iterator
3522  vert_to_send_it = vertices_to_send.begin(),
3523  vert_to_send_end = vertices_to_send.end();
3524  for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3525  ++vert_to_send_it, ++i)
3526  {
3527  int destination = vert_to_send_it->first;
3528  const unsigned int n_vertices = vert_to_send_it->second.size();
3529  const int buffer_size = 2 * n_vertices;
3530  vertices_send_buffers[i].resize(buffer_size);
3531 
3532  // fill the buffer
3533  for (unsigned int j = 0; j < n_vertices; ++j)
3534  {
3535  vertices_send_buffers[i][2 * j] =
3536  std::get<0>(vert_to_send_it->second[j]);
3537  vertices_send_buffers[i][2 * j + 1] =
3538  local_to_global_vertex_index[std::get<1>(
3539  vert_to_send_it->second[j])];
3540  }
3541 
3542  // Send the message
3543  ierr = MPI_Isend(vertices_send_buffers[i].data(),
3544  buffer_size,
3546  destination,
3547  mpi_tag,
3548  triangulation.get_communicator(),
3549  &first_requests[i]);
3550  AssertThrowMPI(ierr);
3551  }
3552 
3553  // Receive the first message
3554  std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
3555  vertices_to_recv.size());
3556  typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
3557  vert_to_recv_it = vertices_to_recv.begin(),
3558  vert_to_recv_end = vertices_to_recv.end();
3559  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3560  ++vert_to_recv_it, ++i)
3561  {
3562  int source = vert_to_recv_it->first;
3563  const unsigned int n_vertices = vert_to_recv_it->second.size();
3564  const int buffer_size = 2 * n_vertices;
3565  vertices_recv_buffers[i].resize(buffer_size);
3566 
3567  // Receive the message
3568  ierr = MPI_Recv(vertices_recv_buffers[i].data(),
3569  buffer_size,
3571  source,
3572  mpi_tag,
3573  triangulation.get_communicator(),
3574  MPI_STATUS_IGNORE);
3575  AssertThrowMPI(ierr);
3576  }
3577 
3578 
3579  // Send second message
3580  std::vector<std::vector<char>> cellids_send_buffers(
3581  vertices_to_send.size());
3582  std::vector<MPI_Request> second_requests(vertices_to_send.size());
3583  vert_to_send_it = vertices_to_send.begin();
3584  for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3585  ++vert_to_send_it, ++i)
3586  {
3587  int destination = vert_to_send_it->first;
3588  const unsigned int n_vertices = vert_to_send_it->second.size();
3589  const int buffer_size = max_cellid_size * n_vertices;
3590  cellids_send_buffers[i].resize(buffer_size);
3591 
3592  // fill the buffer
3593  unsigned int pos = 0;
3594  for (unsigned int j = 0; j < n_vertices; ++j)
3595  {
3596  std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
3597  for (unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
3598  {
3599  if (k < cell_id.size())
3600  cellids_send_buffers[i][pos] = cell_id[k];
3601  // if necessary fill up the reserved part of the buffer with an
3602  // invalid value
3603  else
3604  cellids_send_buffers[i][pos] = '-';
3605  }
3606  }
3607 
3608  // Send the message
3609  ierr = MPI_Isend(cellids_send_buffers[i].data(),
3610  buffer_size,
3611  MPI_CHAR,
3612  destination,
3613  mpi_tag2,
3614  triangulation.get_communicator(),
3615  &second_requests[i]);
3616  AssertThrowMPI(ierr);
3617  }
3618 
3619  // Receive the second message
3620  std::vector<std::vector<char>> cellids_recv_buffers(
3621  vertices_to_recv.size());
3622  vert_to_recv_it = vertices_to_recv.begin();
3623  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3624  ++vert_to_recv_it, ++i)
3625  {
3626  int source = vert_to_recv_it->first;
3627  const unsigned int n_vertices = vert_to_recv_it->second.size();
3628  const int buffer_size = max_cellid_size * n_vertices;
3629  cellids_recv_buffers[i].resize(buffer_size);
3630 
3631  // Receive the message
3632  ierr = MPI_Recv(cellids_recv_buffers[i].data(),
3633  buffer_size,
3634  MPI_CHAR,
3635  source,
3636  mpi_tag2,
3637  triangulation.get_communicator(),
3638  MPI_STATUS_IGNORE);
3639  AssertThrowMPI(ierr);
3640  }
3641 
3642 
3643  // Match the data received with the required vertices
3644  vert_to_recv_it = vertices_to_recv.begin();
3645  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3646  ++i, ++vert_to_recv_it)
3647  {
3648  for (unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
3649  {
3650  const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
3651  const types::global_vertex_index global_id_recv =
3652  vertices_recv_buffers[i][2 * j + 1];
3653  const std::string cellid_recv(
3654  &cellids_recv_buffers[i][max_cellid_size * j],
3655  &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
3656  bool found = false;
3657  typename std::set<active_cell_iterator>::iterator
3658  cell_set_it = missing_vert_cells.begin(),
3659  end_cell_set = missing_vert_cells.end();
3660  for (; (found == false) && (cell_set_it != end_cell_set);
3661  ++cell_set_it)
3662  {
3663  typename std::set<active_cell_iterator>::iterator
3664  candidate_cell =
3665  vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
3666  end_cell =
3667  vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
3668  for (; candidate_cell != end_cell; ++candidate_cell)
3669  {
3670  std::string current_cellid =
3671  (*candidate_cell)->id().to_string();
3672  current_cellid.resize(max_cellid_size, '-');
3673  if (current_cellid.compare(cellid_recv) == 0)
3674  {
3675  local_to_global_vertex_index
3676  [(*candidate_cell)->vertex_index(local_pos_recv)] =
3677  global_id_recv;
3678  found = true;
3679 
3680  break;
3681  }
3682  }
3683  }
3684  }
3685  }
3686 #endif
3687 
3688  return local_to_global_vertex_index;
3689  }
3690 
3691 
3692 
3693  template <int dim, int spacedim>
3694  void
3697  DynamicSparsityPattern & cell_connectivity)
3698  {
3699  cell_connectivity.reinit(triangulation.n_active_cells(),
3700  triangulation.n_active_cells());
3701 
3702  // loop over all cells and their neighbors to build the sparsity
3703  // pattern. note that it's a bit hard to enter all the connections when a
3704  // neighbor has children since we would need to find out which of its
3705  // children is adjacent to the current cell. this problem can be omitted
3706  // if we only do something if the neighbor has no children -- in that case
3707  // it is either on the same or a coarser level than we are. in return, we
3708  // have to add entries in both directions for both cells
3709  for (const auto &cell : triangulation.active_cell_iterators())
3710  {
3711  const unsigned int index = cell->active_cell_index();
3712  cell_connectivity.add(index, index);
3713  for (auto f : cell->face_indices())
3714  if ((cell->at_boundary(f) == false) &&
3715  (cell->neighbor(f)->has_children() == false))
3716  {
3717  const unsigned int other_index =
3718  cell->neighbor(f)->active_cell_index();
3719  cell_connectivity.add(index, other_index);
3720  cell_connectivity.add(other_index, index);
3721  }
3722  }
3723  }
3724 
3725 
3726 
3727  template <int dim, int spacedim>
3728  void
3731  DynamicSparsityPattern & cell_connectivity)
3732  {
3733  std::vector<std::vector<unsigned int>> vertex_to_cell(
3734  triangulation.n_vertices());
3735  for (const auto &cell : triangulation.active_cell_iterators())
3736  {
3737  for (const unsigned int v : cell->vertex_indices())
3738  vertex_to_cell[cell->vertex_index(v)].push_back(
3739  cell->active_cell_index());
3740  }
3741 
3742  cell_connectivity.reinit(triangulation.n_active_cells(),
3743  triangulation.n_active_cells());
3744  for (const auto &cell : triangulation.active_cell_iterators())
3745  {
3746  for (const unsigned int v : cell->vertex_indices())
3747  for (unsigned int n = 0;
3748  n < vertex_to_cell[cell->vertex_index(v)].size();
3749  ++n)
3750  cell_connectivity.add(cell->active_cell_index(),
3751  vertex_to_cell[cell->vertex_index(v)][n]);
3752  }
3753  }
3754 
3755 
3756  template <int dim, int spacedim>
3757  void
3760  const unsigned int level,
3761  DynamicSparsityPattern & cell_connectivity)
3762  {
3763  std::vector<std::vector<unsigned int>> vertex_to_cell(
3764  triangulation.n_vertices());
3765  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3766  triangulation.begin(level);
3767  cell != triangulation.end(level);
3768  ++cell)
3769  {
3770  for (const unsigned int v : cell->vertex_indices())
3771  vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
3772  }
3773 
3774  cell_connectivity.reinit(triangulation.n_cells(level),
3775  triangulation.n_cells(level));
3776  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3777  triangulation.begin(level);
3778  cell != triangulation.end(level);
3779  ++cell)
3780  {
3781  for (const unsigned int v : cell->vertex_indices())
3782  for (unsigned int n = 0;
3783  n < vertex_to_cell[cell->vertex_index(v)].size();
3784  ++n)
3785  cell_connectivity.add(cell->index(),
3786  vertex_to_cell[cell->vertex_index(v)][n]);
3787  }
3788  }
3789 
3790 
3791 
3792  template <int dim, int spacedim>
3793  void
3794  partition_triangulation(const unsigned int n_partitions,
3796  const SparsityTools::Partitioner partitioner)
3797  {
3799  &triangulation) == nullptr),
3800  ExcMessage("Objects of type parallel::distributed::Triangulation "
3801  "are already partitioned implicitly and can not be "
3802  "partitioned again explicitly."));
3803 
3804  std::vector<unsigned int> cell_weights;
3805 
3806  // Get cell weighting if a signal has been attached to the triangulation
3807  if (!triangulation.signals.cell_weight.empty())
3808  {
3809  cell_weights.resize(triangulation.n_active_cells(), 0U);
3810 
3811  // In a first step, obtain the weights of the locally owned
3812  // cells. For all others, the weight remains at the zero the
3813  // vector was initialized with above.
3814  for (const auto &cell : triangulation.active_cell_iterators())
3815  if (cell->is_locally_owned())
3816  cell_weights[cell->active_cell_index()] =
3817  triangulation.signals.cell_weight(
3819 
3820  // If this is a parallel triangulation, we then need to also
3821  // get the weights for all other cells. We have asserted above
3822  // that this function can't be used for
3823  // parallel::distribute::Triangulation objects, so the only
3824  // ones we have to worry about here are
3825  // parallel::shared::Triangulation
3826  if (const auto shared_tria =
3828  &triangulation))
3829  Utilities::MPI::sum(cell_weights,
3830  shared_tria->get_communicator(),
3831  cell_weights);
3832  }
3833 
3834  // Call the other more general function
3835  partition_triangulation(n_partitions,
3836  cell_weights,
3837  triangulation,
3838  partitioner);
3839  }
3840 
3841 
3842 
3843  template <int dim, int spacedim>
3844  void
3845  partition_triangulation(const unsigned int n_partitions,
3846  const std::vector<unsigned int> &cell_weights,
3848  const SparsityTools::Partitioner partitioner)
3849  {
3851  &triangulation) == nullptr),
3852  ExcMessage("Objects of type parallel::distributed::Triangulation "
3853  "are already partitioned implicitly and can not be "
3854  "partitioned again explicitly."));
3855  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
3856 
3857  // check for an easy return
3858  if (n_partitions == 1)
3859  {
3860  for (const auto &cell : triangulation.active_cell_iterators())
3861  cell->set_subdomain_id(0);
3862  return;
3863  }
3864 
3865  // we decompose the domain by first
3866  // generating the connection graph of all
3867  // cells with their neighbors, and then
3868  // passing this graph off to METIS.
3869  // finally defer to the other function for
3870  // partitioning and assigning subdomain ids
3871  DynamicSparsityPattern cell_connectivity;
3872  get_face_connectivity_of_cells(triangulation, cell_connectivity);
3873 
3874  SparsityPattern sp_cell_connectivity;
3875  sp_cell_connectivity.copy_from(cell_connectivity);
3876  partition_triangulation(n_partitions,
3877  cell_weights,
3878  sp_cell_connectivity,
3879  triangulation,
3880  partitioner);
3881  }
3882 
3883 
3884 
3885  template <int dim, int spacedim>
3886  void
3887  partition_triangulation(const unsigned int n_partitions,
3888  const SparsityPattern & cell_connection_graph,
3890  const SparsityTools::Partitioner partitioner)
3891  {
3893  &triangulation) == nullptr),
3894  ExcMessage("Objects of type parallel::distributed::Triangulation "
3895  "are already partitioned implicitly and can not be "
3896  "partitioned again explicitly."));
3897 
3898  std::vector<unsigned int> cell_weights;
3899 
3900  // Get cell weighting if a signal has been attached to the triangulation
3901  if (!triangulation.signals.cell_weight.empty())
3902  {
3903  cell_weights.resize(triangulation.n_active_cells(), 0U);
3904 
3905  // In a first step, obtain the weights of the locally owned
3906  // cells. For all others, the weight remains at the zero the
3907  // vector was initialized with above.
3908  for (const auto &cell : triangulation.active_cell_iterators())
3909  if (cell->is_locally_owned())
3910  cell_weights[cell->active_cell_index()] =
3911  triangulation.signals.cell_weight(
3913 
3914  // If this is a parallel triangulation, we then need to also
3915  // get the weights for all other cells. We have asserted above
3916  // that this function can't be used for
3917  // parallel::distribute::Triangulation objects, so the only
3918  // ones we have to worry about here are
3919  // parallel::shared::Triangulation
3920  if (const auto shared_tria =
3922  &triangulation))
3923  Utilities::MPI::sum(cell_weights,
3924  shared_tria->get_communicator(),
3925  cell_weights);
3926  }
3927 
3928  // Call the other more general function
3929  partition_triangulation(n_partitions,
3930  cell_weights,
3931  cell_connection_graph,
3932  triangulation,
3933  partitioner);
3934  }
3935 
3936 
3937 
3938  template <int dim, int spacedim>
3939  void
3940  partition_triangulation(const unsigned int n_partitions,
3941  const std::vector<unsigned int> &cell_weights,
3942  const SparsityPattern & cell_connection_graph,
3944  const SparsityTools::Partitioner partitioner)
3945  {
3947  &triangulation) == nullptr),
3948  ExcMessage("Objects of type parallel::distributed::Triangulation "
3949  "are already partitioned implicitly and can not be "
3950  "partitioned again explicitly."));
3951  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
3952  Assert(cell_connection_graph.n_rows() == triangulation.n_active_cells(),
3953  ExcMessage("Connectivity graph has wrong size"));
3954  Assert(cell_connection_graph.n_cols() == triangulation.n_active_cells(),
3955  ExcMessage("Connectivity graph has wrong size"));
3956 
3957  // signal that partitioning is going to happen
3958  triangulation.signals.pre_partition();
3959 
3960  // check for an easy return
3961  if (n_partitions == 1)
3962  {
3963  for (const auto &cell : triangulation.active_cell_iterators())
3964  cell->set_subdomain_id(0);
3965  return;
3966  }
3967 
3968  // partition this connection graph and get
3969  // back a vector of indices, one per degree
3970  // of freedom (which is associated with a
3971  // cell)
3972  std::vector<unsigned int> partition_indices(triangulation.n_active_cells());
3973  SparsityTools::partition(cell_connection_graph,
3974  cell_weights,
3975  n_partitions,
3976  partition_indices,
3977  partitioner);
3978 
3979  // finally loop over all cells and set the subdomain ids
3980  for (const auto &cell : triangulation.active_cell_iterators())
3981  cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
3982  }
3983 
3984 
3985  namespace internal
3986  {
3990  template <class IT>
3991  void
3993  unsigned int & current_proc_idx,
3994  unsigned int & current_cell_idx,
3995  const unsigned int n_active_cells,
3996  const unsigned int n_partitions)
3997  {
3998  if (cell->is_active())
3999  {
4000  while (current_cell_idx >=
4001  std::floor(static_cast<uint_least64_t>(n_active_cells) *
4002  (current_proc_idx + 1) / n_partitions))
4003  ++current_proc_idx;
4004  cell->set_subdomain_id(current_proc_idx);
4005  ++current_cell_idx;
4006  }
4007  else
4008  {
4009  for (unsigned int n = 0; n < cell->n_children(); ++n)
4011  current_proc_idx,
4012  current_cell_idx,
4014  n_partitions);
4015  }
4016  }
4017  } // namespace internal
4018 
4019  template <int dim, int spacedim>
4020  void
4021  partition_triangulation_zorder(const unsigned int n_partitions,
4023  const bool group_siblings)
4024  {
4026  &triangulation) == nullptr),
4027  ExcMessage("Objects of type parallel::distributed::Triangulation "
4028  "are already partitioned implicitly and can not be "
4029  "partitioned again explicitly."));
4030  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
4031 
4032  // signal that partitioning is going to happen
4033  triangulation.signals.pre_partition();
4034 
4035  // check for an easy return
4036  if (n_partitions == 1)
4037  {
4038  for (const auto &cell : triangulation.active_cell_iterators())
4039  cell->set_subdomain_id(0);
4040  return;
4041  }
4042 
4043  // Duplicate the coarse cell reordoring
4044  // as done in p4est
4045  std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
4046  std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
4047 
4048  DynamicSparsityPattern cell_connectivity;
4050  0,
4051  cell_connectivity);
4052  coarse_cell_to_p4est_tree_permutation.resize(triangulation.n_cells(0));
4053  SparsityTools::reorder_hierarchical(cell_connectivity,
4054  coarse_cell_to_p4est_tree_permutation);
4055 
4056  p4est_tree_to_coarse_cell_permutation =
4057  Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
4058 
4059  unsigned int current_proc_idx = 0;
4060  unsigned int current_cell_idx = 0;
4061  const unsigned int n_active_cells = triangulation.n_active_cells();
4062 
4063  // set subdomain id for active cell descendants
4064  // of each coarse cell in permuted order
4065  for (unsigned int idx = 0; idx < triangulation.n_cells(0); ++idx)
4066  {
4067  const unsigned int coarse_cell_idx =
4068  p4est_tree_to_coarse_cell_permutation[idx];
4069  typename Triangulation<dim, spacedim>::cell_iterator coarse_cell(
4070  &triangulation, 0, coarse_cell_idx);
4071 
4073  current_proc_idx,
4074  current_cell_idx,
4075  n_active_cells,
4076  n_partitions);
4077  }
4078 
4079  // if all children of a cell are active (e.g. we
4080  // have a cell that is refined once and no part
4081  // is refined further), p4est places all of them
4082  // on the same processor. The new owner will be
4083  // the processor with the largest number of children
4084  // (ties are broken by picking the lower rank).
4085  // Duplicate this logic here.
4086  if (group_siblings)
4087  {
4089  cell = triangulation.begin(),
4090  endc = triangulation.end();
4091  for (; cell != endc; ++cell)
4092  {
4093  if (cell->is_active())
4094  continue;
4095  bool all_children_active = true;
4096  std::map<unsigned int, unsigned int> map_cpu_n_cells;
4097  for (unsigned int n = 0; n < cell->n_children(); ++n)
4098  if (!cell->child(n)->is_active())
4099  {
4100  all_children_active = false;
4101  break;
4102  }
4103  else
4104  ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
4105 
4106  if (!all_children_active)
4107  continue;
4108 
4109  unsigned int new_owner = cell->child(0)->subdomain_id();
4110  for (std::map<unsigned int, unsigned int>::iterator it =
4111  map_cpu_n_cells.begin();
4112  it != map_cpu_n_cells.end();
4113  ++it)
4114  if (it->second > map_cpu_n_cells[new_owner])
4115  new_owner = it->first;
4116 
4117  for (unsigned int n = 0; n < cell->n_children(); ++n)
4118  cell->child(n)->set_subdomain_id(new_owner);
4119  }
4120  }
4121  }
4122 
4123 
4124  template <int dim, int spacedim>
4125  void
4127  {
4128  unsigned int n_levels = triangulation.n_levels();
4129  for (int lvl = n_levels - 1; lvl >= 0; --lvl)
4130  {
4132  cell = triangulation.begin(lvl),
4133  endc = triangulation.end(lvl);
4134  for (; cell != endc; ++cell)
4135  {
4136  if (cell->is_active())
4137  cell->set_level_subdomain_id(cell->subdomain_id());
4138  else
4139  {
4140  Assert(cell->child(0)->level_subdomain_id() !=
4142  ExcInternalError());
4143  cell->set_level_subdomain_id(
4144  cell->child(0)->level_subdomain_id());
4145  }
4146  }
4147  }
4148  }
4149 
4150  namespace internal
4151  {
4152  namespace
4153  {
4154  // Split get_subdomain_association() for p::d::T since we want to compile
4155  // it in 1D but none of the p4est stuff is available in 1D.
4156  template <int dim, int spacedim>
4157  void
4160  & triangulation,
4161  const std::vector<CellId> & cell_ids,
4162  std::vector<types::subdomain_id> &subdomain_ids)
4163  {
4164 #ifndef DEAL_II_WITH_P4EST
4165  (void)triangulation;
4166  (void)cell_ids;
4167  (void)subdomain_ids;
4168  Assert(
4169  false,
4170  ExcMessage(
4171  "You are attempting to use a functionality that is only available "
4172  "if deal.II was configured to use p4est, but cmake did not find a "
4173  "valid p4est library."));
4174 #else
4175  // for parallel distributed triangulations, we will ask the p4est oracle
4176  // about the global partitioning of active cells since this information
4177  // is stored on every process
4178  for (const auto &cell_id : cell_ids)
4179  {
4180  // find descendent from coarse quadrant
4181  typename ::internal::p4est::types<dim>::quadrant p4est_cell,
4183 
4184  ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
4185  for (const auto &child_index : cell_id.get_child_indices())
4186  {
4187  ::internal::p4est::init_quadrant_children<dim>(
4188  p4est_cell, p4est_children);
4189  p4est_cell =
4190  p4est_children[static_cast<unsigned int>(child_index)];
4191  }
4192 
4193  // find owning process, i.e., the subdomain id
4194  const int owner =
4196  const_cast<typename ::internal::p4est::types<dim>::forest
4197  *>(triangulation.get_p4est()),
4198  cell_id.get_coarse_cell_id(),
4199  &p4est_cell,
4201  triangulation.get_communicator()));
4202 
4203  Assert(owner >= 0, ExcMessage("p4est should know the owner."));
4204 
4205  subdomain_ids.push_back(owner);
4206  }
4207 #endif
4208  }
4209 
4210 
4211 
4212  template <int spacedim>
4213  void
4216  const std::vector<CellId> &,
4217  std::vector<types::subdomain_id> &)
4218  {
4219  Assert(false, ExcNotImplemented());
4220  }
4221  } // anonymous namespace
4222  } // namespace internal
4223 
4224 
4225 
4226  template <int dim, int spacedim>
4227  std::vector<types::subdomain_id>
4229  const std::vector<CellId> & cell_ids)
4230  {
4231  std::vector<types::subdomain_id> subdomain_ids;
4232  subdomain_ids.reserve(cell_ids.size());
4233 
4234  if (dynamic_cast<
4236  &triangulation) != nullptr)
4237  {
4238  Assert(false, ExcNotImplemented());
4239  }
4241  *parallel_tria = dynamic_cast<
4243  &triangulation))
4244  {
4245  internal::get_subdomain_association(*parallel_tria,
4246  cell_ids,
4247  subdomain_ids);
4248  }
4249  else if (const parallel::shared::Triangulation<dim, spacedim> *shared_tria =
4251  *>(&triangulation))
4252  {
4253  // for parallel shared triangulations, we need to access true subdomain
4254  // ids which are also valid for artificial cells
4255  const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
4256  shared_tria->get_true_subdomain_ids_of_cells();
4257 
4258  for (const auto &cell_id : cell_ids)
4259  {
4260  const unsigned int active_cell_index =
4261  shared_tria->create_cell_iterator(cell_id)->active_cell_index();
4262  subdomain_ids.push_back(
4263  true_subdomain_ids_of_cells[active_cell_index]);
4264  }
4265  }
4266  else
4267  {
4268  // the most general type of triangulation is the serial one. here, all
4269  // subdomain information is directly available
4270  for (const auto &cell_id : cell_ids)
4271  {
4272  subdomain_ids.push_back(
4273  triangulation.create_cell_iterator(cell_id)->subdomain_id());
4274  }
4275  }
4276 
4277  return subdomain_ids;
4278  }
4279 
4280 
4281 
4282  template <int dim, int spacedim>
4283  void
4285  std::vector<types::subdomain_id> & subdomain)
4286  {
4287  Assert(subdomain.size() == triangulation.n_active_cells(),
4288  ExcDimensionMismatch(subdomain.size(),
4289  triangulation.n_active_cells()));
4290  for (const auto &cell : triangulation.active_cell_iterators())
4291  subdomain[cell->active_cell_index()] = cell->subdomain_id();
4292  }
4293 
4294 
4295 
4296  template <int dim, int spacedim>
4297  unsigned int
4300  const types::subdomain_id subdomain)
4301  {
4302  unsigned int count = 0;
4303  for (const auto &cell : triangulation.active_cell_iterators())
4304  if (cell->subdomain_id() == subdomain)
4305  ++count;
4306 
4307  return count;
4308  }
4309 
4310 
4311 
4312  template <int dim, int spacedim>
4313  std::vector<bool>
4315  {
4316  // start with all vertices
4317  std::vector<bool> locally_owned_vertices =
4318  triangulation.get_used_vertices();
4319 
4320  // if the triangulation is distributed, eliminate those that
4321  // are owned by other processors -- either because the vertex is
4322  // on an artificial cell, or because it is on a ghost cell with
4323  // a smaller subdomain
4324  if (const auto *tr = dynamic_cast<
4326  &triangulation))
4327  for (const auto &cell : triangulation.active_cell_iterators())
4328  if (cell->is_artificial() ||
4329  (cell->is_ghost() &&
4330  (cell->subdomain_id() < tr->locally_owned_subdomain())))
4331  for (const unsigned int v : cell->vertex_indices())
4332  locally_owned_vertices[cell->vertex_index(v)] = false;
4333 
4334  return locally_owned_vertices;
4335  }
4336 
4337 
4338 
4339  template <int dim, int spacedim>
4340  double
4342  const Mapping<dim, spacedim> & mapping)
4343  {
4344  double min_diameter = std::numeric_limits<double>::max();
4345  for (const auto &cell : triangulation.active_cell_iterators())
4346  if (!cell->is_artificial())
4347  min_diameter = std::min(min_diameter, cell->diameter(mapping));
4348 
4349  double global_min_diameter = 0;
4350 
4351 #ifdef DEAL_II_WITH_MPI
4352  if (const parallel::TriangulationBase<dim, spacedim> *p_tria =
4353  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4354  &triangulation))
4355  global_min_diameter =
4356  Utilities::MPI::min(min_diameter, p_tria->get_communicator());
4357  else
4358 #endif
4359  global_min_diameter = min_diameter;
4360 
4361  return global_min_diameter;
4362  }
4363 
4364 
4365 
4366  template <int dim, int spacedim>
4367  double
4369  const Mapping<dim, spacedim> & mapping)
4370  {
4371  double max_diameter = 0.;
4372  for (const auto &cell : triangulation.active_cell_iterators())
4373  if (!cell->is_artificial())
4374  max_diameter = std::max(max_diameter, cell->diameter(mapping));
4375 
4376  double global_max_diameter = 0;
4377 
4378 #ifdef DEAL_II_WITH_MPI
4379  if (const parallel::TriangulationBase<dim, spacedim> *p_tria =
4380  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4381  &triangulation))
4382  global_max_diameter =
4383  Utilities::MPI::max(max_diameter, p_tria->get_communicator());
4384  else
4385 #endif
4386  global_max_diameter = max_diameter;
4387 
4388  return global_max_diameter;
4389  }
4390 
4391 
4392 
4393  namespace internal
4394  {
4395  namespace FixUpDistortedChildCells
4396  {
4397  // compute the mean square
4398  // deviation of the alternating
4399  // forms of the children of the
4400  // given object from that of
4401  // the object itself. for
4402  // objects with
4403  // structdim==spacedim, the
4404  // alternating form is the
4405  // determinant of the jacobian,
4406  // whereas for faces with
4407  // structdim==spacedim-1, the
4408  // alternating form is the
4409  // (signed and scaled) normal
4410  // vector
4411  //
4412  // this average square
4413  // deviation is computed for an
4414  // object where the center node
4415  // has been replaced by the
4416  // second argument to this
4417  // function
4418  template <typename Iterator, int spacedim>
4419  double
4420  objective_function(const Iterator & object,
4421  const Point<spacedim> &object_mid_point)
4422  {
4423  const unsigned int structdim =
4424  Iterator::AccessorType::structure_dimension;
4425  Assert(spacedim == Iterator::AccessorType::dimension,
4426  ExcInternalError());
4427 
4428  // everything below is wrong
4429  // if not for the following
4430  // condition
4431  Assert(object->refinement_case() ==
4433  ExcNotImplemented());
4434  // first calculate the
4435  // average alternating form
4436  // for the parent cell/face
4439  Tensor<spacedim - structdim, spacedim>
4440  parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
4441 
4442  for (const unsigned int i : object->vertex_indices())
4443  parent_vertices[i] = object->vertex(i);
4444 
4446  parent_vertices, parent_alternating_forms);
4447 
4448  const Tensor<spacedim - structdim, spacedim>
4449  average_parent_alternating_form =
4450  std::accumulate(parent_alternating_forms,
4451  parent_alternating_forms +
4454 
4455  // now do the same
4456  // computation for the
4457  // children where we use the
4458  // given location for the
4459  // object mid point instead of
4460  // the one the triangulation
4461  // currently reports
4465  Tensor<spacedim - structdim, spacedim> child_alternating_forms
4468 
4469  for (unsigned int c = 0; c < object->n_children(); ++c)
4470  for (const unsigned int i : object->child(c)->vertex_indices())
4471  child_vertices[c][i] = object->child(c)->vertex(i);
4472 
4473  // replace mid-object
4474  // vertex. note that for
4475  // child i, the mid-object
4476  // vertex happens to have the
4477  // number
4478  // max_children_per_cell-i
4479  for (unsigned int c = 0; c < object->n_children(); ++c)
4480  child_vertices[c][GeometryInfo<structdim>::max_children_per_cell - c -
4481  1] = object_mid_point;
4482 
4483  for (unsigned int c = 0; c < object->n_children(); ++c)
4485  child_vertices[c], child_alternating_forms[c]);
4486 
4487  // on a uniformly refined
4488  // hypercube object, the child
4489  // alternating forms should
4490  // all be smaller by a factor
4491  // of 2^structdim than the
4492  // ones of the parent. as a
4493  // consequence, we'll use the
4494  // squared deviation from
4495  // this ideal value as an
4496  // objective function
4497  double objective = 0;
4498  for (unsigned int c = 0; c < object->n_children(); ++c)
4499  for (const unsigned int i : object->child(c)->vertex_indices())
4500  objective +=
4501  (child_alternating_forms[c][i] -
4502  average_parent_alternating_form / std::pow(2., 1. * structdim))
4503  .norm_square();
4504 
4505  return objective;
4506  }
4507 
4508 
4514  template <typename Iterator>
4516  get_face_midpoint(const Iterator & object,
4517  const unsigned int f,
4518  std::integral_constant<int, 1>)
4519  {
4520  return object->vertex(f);
4521  }
4522 
4523 
4524 
4530  template <typename Iterator>
4532  get_face_midpoint(const Iterator & object,
4533  const unsigned int f,
4534  std::integral_constant<int, 2>)
4535  {
4536  return object->line(f)->center();
4537  }
4538 
4539 
4540 
4546  template <typename Iterator>
4548  get_face_midpoint(const Iterator & object,
4549  const unsigned int f,
4550  std::integral_constant<int, 3>)
4551  {
4552  return object->face(f)->center();
4553  }
4554 
4555 
4556 
4579  template <typename Iterator>
4580  double
4581  minimal_diameter(const Iterator &object)
4582  {
4583  const unsigned int structdim =
4584  Iterator::AccessorType::structure_dimension;
4585 
4586  double diameter = object->diameter();
4587  for (const unsigned int f : object->face_indices())
4588  for (unsigned int e = f + 1; e < object->n_faces(); ++e)
4589  diameter = std::min(
4590  diameter,
4591  get_face_midpoint(object,
4592  f,
4593  std::integral_constant<int, structdim>())
4594  .distance(get_face_midpoint(
4595  object, e, std::integral_constant<int, structdim>())));
4596 
4597  return diameter;
4598  }
4599 
4600 
4601 
4606  template <typename Iterator>
4607  bool
4608  fix_up_object(const Iterator &object)
4609  {
4610  const unsigned int structdim =
4611  Iterator::AccessorType::structure_dimension;
4612  const unsigned int spacedim = Iterator::AccessorType::space_dimension;
4613 
4614  // right now we can only deal with cells that have been refined
4615  // isotropically because that is the only case where we have a cell
4616  // mid-point that can be moved around without having to consider
4617  // boundary information
4618  Assert(object->has_children(), ExcInternalError());
4619  Assert(object->refinement_case() ==
4621  ExcNotImplemented());
4622 
4623  // get the current location of the object mid-vertex:
4624  Point<spacedim> object_mid_point = object->child(0)->vertex(
4626 
4627  // now do a few steepest descent steps to reduce the objective
4628  // function. compute the diameter in the helper function above
4629  unsigned int iteration = 0;
4630  const double diameter = minimal_diameter(object);
4631 
4632  // current value of objective function and initial delta
4633  double current_value = objective_function(object, object_mid_point);
4634  double initial_delta = 0;
4635 
4636  do
4637  {
4638  // choose a step length that is initially 1/4 of the child
4639  // objects' diameter, and a sequence whose sum does not converge
4640  // (to avoid premature termination of the iteration)
4641  const double step_length = diameter / 4 / (iteration + 1);
4642 
4643  // compute the objective function's derivative using a two-sided
4644  // difference formula with eps=step_length/10
4645  Tensor<1, spacedim> gradient;
4646  for (unsigned int d = 0; d < spacedim; ++d)
4647  {
4648  const double eps = step_length / 10;
4649 
4651  h[d] = eps / 2;
4652 
4653  gradient[d] =
4655  object, project_to_object(object, object_mid_point + h)) -
4657  object, project_to_object(object, object_mid_point - h))) /
4658  eps;
4659  }
4660 
4661  // there is nowhere to go
4662  if (gradient.norm() == 0)
4663  break;
4664 
4665  // We need to go in direction -gradient. the optimal value of the
4666  // objective function is zero, so assuming that the model is
4667  // quadratic we would have to go -2*val/||gradient|| in this
4668  // direction, make sure we go at most step_length into this
4669  // direction
4670  object_mid_point -=
4671  std::min(2 * current_value / (gradient * gradient),
4672  step_length / gradient.norm()) *
4673  gradient;
4674  object_mid_point = project_to_object(object, object_mid_point);
4675 
4676  // compute current value of the objective function
4677  const double previous_value = current_value;
4678  current_value = objective_function(object, object_mid_point);
4679 
4680  if (iteration == 0)
4681  initial_delta = (previous_value - current_value);
4682 
4683  // stop if we aren't moving much any more
4684  if ((iteration >= 1) &&
4685  ((previous_value - current_value < 0) ||
4686  (std::fabs(previous_value - current_value) <
4687  0.001 * initial_delta)))
4688  break;
4689 
4690  ++iteration;
4691  }
4692  while (iteration < 20);
4693 
4694  // verify that the new
4695  // location is indeed better
4696  // than the one before. check
4697  // this by comparing whether
4698  // the minimum value of the
4699  // products of parent and
4700  // child alternating forms is
4701  // positive. for cells this
4702  // means that the
4703  // determinants have the same
4704  // sign, for faces that the
4705  // face normals of parent and
4706  // children point in the same
4707  // general direction
4708  double old_min_product, new_min_product;
4709 
4712  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
4713  parent_vertices[i] = object->vertex(i);
4714 
4715  Tensor<spacedim - structdim, spacedim>
4716  parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
4718  parent_vertices, parent_alternating_forms);
4719 
4723 
4724  for (unsigned int c = 0; c < object->n_children(); ++c)
4725  for (const unsigned int i : object->child(c)->vertex_indices())
4726  child_vertices[c][i] = object->child(c)->vertex(i);
4727 
4728  Tensor<spacedim - structdim, spacedim> child_alternating_forms
4731 
4732  for (unsigned int c = 0; c < object->n_children(); ++c)
4734  child_vertices[c], child_alternating_forms[c]);
4735 
4736  old_min_product =
4737  child_alternating_forms[0][0] * parent_alternating_forms[0];
4738  for (unsigned int c = 0; c < object->n_children(); ++c)
4739  for (const unsigned int i : object->child(c)->vertex_indices())
4740  for (const unsigned int j : object->vertex_indices())
4741  old_min_product = std::min<double>(old_min_product,
4742  child_alternating_forms[c][i] *
4743  parent_alternating_forms[j]);
4744 
4745  // for the new minimum value,
4746  // replace mid-object
4747  // vertex. note that for child
4748  // i, the mid-object vertex
4749  // happens to have the number
4750  // max_children_per_cell-i
4751  for (unsigned int c = 0; c < object->n_children(); ++c)
4752  child_vertices[c][GeometryInfo<structdim>::max_children_per_cell - c -
4753  1] = object_mid_point;
4754 
4755  for (unsigned int c = 0; c < object->n_children(); ++c)
4757  child_vertices[c], child_alternating_forms[c]);
4758 
4759  new_min_product =
4760  child_alternating_forms[0][0] * parent_alternating_forms[0];
4761  for (unsigned int c = 0; c < object->n_children(); ++c)
4762  for (const unsigned int i : object->child(c)->vertex_indices())
4763  for (const unsigned int j : object->vertex_indices())
4764  new_min_product = std::min<double>(new_min_product,
4765  child_alternating_forms[c][i] *
4766  parent_alternating_forms[j]);
4767 
4768  // if new minimum value is
4769  // better than before, then set the
4770  // new mid point. otherwise
4771  // return this object as one of
4772  // those that can't apparently
4773  // be fixed
4774  if (new_min_product >= old_min_product)
4775  object->child(0)->vertex(
4777  object_mid_point;
4778 
4779  // return whether after this
4780  // operation we have an object that
4781  // is well oriented
4782  return (std::max(new_min_product, old_min_product) > 0);
4783  }
4784 
4785 
4786 
4787  // possibly fix up the faces of a cell by moving around its mid-points
4788  template <int dim, int spacedim>
4789  void
4791  const typename ::Triangulation<dim, spacedim>::cell_iterator
4792  &cell,
4793  std::integral_constant<int, dim>,
4794  std::integral_constant<int, spacedim>)
4795  {
4796  // see if we first can fix up some of the faces of this object. We can
4797  // mess with faces if and only if the neighboring cell is not even
4798  // more refined than we are (since in that case the sub-faces have
4799  // themselves children that we can't move around any more). however,
4800  // the latter case shouldn't happen anyway: if the current face is
4801  // distorted but the neighbor is even more refined, then the face had
4802  // been deformed before already, and had been ignored at the time; we
4803  // should then also be able to ignore it this time as well
4804  for (auto f : cell->face_indices())
4805  {
4806  Assert(cell->face(f)->has_children(), ExcInternalError());
4807  Assert(cell->face(f)->refinement_case() ==
4808  RefinementCase<dim - 1>::isotropic_refinement,
4809  ExcInternalError());
4810 
4811  bool subface_is_more_refined = false;
4812  for (unsigned int g = 0;
4813  g < GeometryInfo<dim>::max_children_per_face;
4814  ++g)
4815  if (cell->face(f)->child(g)->has_children())
4816  {
4817  subface_is_more_refined = true;
4818  break;
4819  }
4820 
4821  if (subface_is_more_refined == true)
4822  continue;
4823 
4824  // we finally know that we can do something about this face
4825  fix_up_object(cell->face(f));
4826  }
4827  }
4828  } /* namespace FixUpDistortedChildCells */
4829  } /* namespace internal */
4830 
4831 
4832  template <int dim, int spacedim>
4836  &distorted_cells,
4837  Triangulation<dim, spacedim> & /*triangulation*/)
4838  {
4839  static_assert(
4840  dim != 1 && spacedim != 1,
4841  "This function is only valid when dim != 1 or spacedim != 1.");
4842  typename Triangulation<dim, spacedim>::DistortedCellList unfixable_subset;
4843 
4844  // loop over all cells that we have to fix up
4845  for (typename std::list<
4846  typename Triangulation<dim, spacedim>::cell_iterator>::const_iterator
4847  cell_ptr = distorted_cells.distorted_cells.begin();
4848  cell_ptr != distorted_cells.distorted_cells.end();
4849  ++cell_ptr)
4850  {
4851  const typename Triangulation<dim, spacedim>::cell_iterator cell =
4852  *cell_ptr;
4853 
4854  Assert(!cell->is_active(),
4855  ExcMessage(
4856  "This function is only valid for a list of cells that "
4857  "have children (i.e., no cell in the list may be active)."));
4858 
4860  cell,
4861  std::integral_constant<int, dim>(),
4862  std::integral_constant<int, spacedim>());
4863 
4864  // If possible, fix up the object.
4866  unfixable_subset.distorted_cells.push_back(cell);
4867  }
4868 
4869  return unfixable_subset;
4870  }
4871 
4872 
4873 
4874  template <int dim, int spacedim>
4875  void
4877  const bool reset_boundary_ids)
4878  {
4879  const auto src_boundary_ids = tria.get_boundary_ids();
4880  std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
4881  auto m_it = dst_manifold_ids.begin();
4882  for (const auto b : src_boundary_ids)
4883  {
4884  *m_it = static_cast<types::manifold_id>(b);
4885  ++m_it;
4886  }
4887  const std::vector<types::boundary_id> reset_boundary_id =
4888  reset_boundary_ids ?
4889  std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
4890  src_boundary_ids;
4891  map_boundary_to_manifold_ids(src_boundary_ids,
4892  dst_manifold_ids,
4893  tria,
4894  reset_boundary_id);
4895  }
4896 
4897 
4898 
4899  template <int dim, int spacedim>
4900  void
4902  const std::vector<types::boundary_id> &src_boundary_ids,
4903  const std::vector<types::manifold_id> &dst_manifold_ids,
4905  const std::vector<types::boundary_id> &reset_boundary_ids_)
4906  {
4907  AssertDimension(src_boundary_ids.size(), dst_manifold_ids.size());
4908  const auto reset_boundary_ids =
4909  reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
4910  AssertDimension(reset_boundary_ids.size(), src_boundary_ids.size());
4911 
4912  // in 3d, we not only have to copy boundary ids of faces, but also of edges
4913  // because we see them twice (once from each adjacent boundary face),
4914  // we cannot immediately reset their boundary ids. thus, copy first
4915  // and reset later
4916  if (dim >= 3)
4917  for (const auto &cell : tria.active_cell_iterators())
4918  for (auto f : cell->face_indices())
4919  if (cell->face(f)->at_boundary())
4920  for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
4921  {
4922  const auto bid = cell->face(f)->line(e)->boundary_id();
4923  const unsigned int ind = std::find(src_boundary_ids.begin(),
4924  src_boundary_ids.end(),
4925  bid) -
4926  src_boundary_ids.begin();
4927  if (ind < src_boundary_ids.size())
4928  cell->face(f)->line(e)->set_manifold_id(
4929  dst_manifold_ids[ind]);
4930  }
4931 
4932  // now do cells
4933  for (const auto &cell : tria.active_cell_iterators())
4934  for (auto f : cell->face_indices())
4935  if (cell->face(f)->at_boundary())
4936  {
4937  const auto bid = cell->face(f)->boundary_id();
4938  const unsigned int ind =
4939  std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
4940  src_boundary_ids.begin();
4941 
4942  if (ind < src_boundary_ids.size())
4943  {
4944  // assign the manifold id
4945  cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
4946  // then reset boundary id
4947  cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
4948  }
4949 
4950  if (dim >= 3)
4951  for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
4952  {
4953  const auto bid = cell->face(f)->line(e)->boundary_id();
4954  const unsigned int ind = std::find(src_boundary_ids.begin(),
4955  src_boundary_ids.end(),
4956  bid) -
4957  src_boundary_ids.begin();
4958  if (ind < src_boundary_ids.size())
4959  cell->face(f)->line(e)->set_boundary_id(
4960  reset_boundary_ids[ind]);
4961  }
4962  }
4963  }
4964 
4965 
4966  template <int dim, int spacedim>
4967  void
4969  const bool compute_face_ids)
4970  {
4972  cell = tria.begin_active(),
4973  endc = tria.end();
4974 
4975  for (; cell != endc; ++cell)
4976  {
4977  cell->set_manifold_id(cell->material_id());
4978  if (compute_face_ids == true)
4979  {
4980  for (auto f : cell->face_indices())
4981  {
4982  if (cell->at_boundary(f) == false)
4983  cell->face(f)->set_manifold_id(
4984  std::min(cell->material_id(),
4985  cell->neighbor(f)->material_id()));
4986  else
4987  cell->face(f)->set_manifold_id(cell->material_id());
4988  }
4989  }
4990  }
4991  }
4992 
4993 
4994  template <int dim, int spacedim>
4995  void
4998  const std::function<types::manifold_id(
4999  const std::set<types::manifold_id> &)> &disambiguation_function,
5000  bool overwrite_only_flat_manifold_ids)
5001  {
5002  // Easy case first:
5003  if (dim == 1)
5004  return;
5005  const unsigned int n_subobjects =
5006  dim == 2 ? tria.n_lines() : tria.n_lines() + tria.n_quads();
5007 
5008  // If user index is zero, then it has not been set.
5009  std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
5010  std::vector<unsigned int> backup;
5011  tria.save_user_indices(backup);
5012  tria.clear_user_data();
5013 
5014  unsigned next_index = 1;
5015  for (auto &cell : tria.active_cell_iterators())
5016  {
5017  if (dim > 1)
5018  for (unsigned int l = 0; l < cell->n_lines(); ++l)
5019  {
5020  if (cell->line(l)->user_index() == 0)
5021  {
5022  AssertIndexRange(next_index, n_subobjects + 1);
5023  manifold_ids[next_index].insert(cell->manifold_id());
5024  cell->line(l)->set_user_index(next_index++);
5025  }
5026  else
5027  manifold_ids[cell->line(l)->user_index()].insert(
5028  cell->manifold_id());
5029  }
5030  if (dim > 2)
5031  for (unsigned int l = 0; l < cell->n_faces(); ++l)
5032  {
5033  if (cell->quad(l)->user_index() == 0)
5034  {
5035  AssertIndexRange(next_index, n_subobjects + 1);
5036  manifold_ids[next_index].insert(cell->manifold_id());
5037  cell->quad(l)->set_user_index(next_index++);
5038  }
5039  else
5040  manifold_ids[cell->quad(l)->user_index()].insert(
5041  cell->manifold_id());
5042  }
5043  }
5044  for (auto &cell : tria.active_cell_iterators())
5045  {
5046  if (dim > 1)
5047  for (unsigned int l = 0; l < cell->n_lines(); ++l)
5048  {
5049  const auto id = cell->line(l)->user_index();
5050  // Make sure we change the manifold indicator only once
5051  if (id != 0)
5052  {
5053  if (cell->line(l)->manifold_id() ==
5055  overwrite_only_flat_manifold_ids == false)
5056  cell->line(l)->set_manifold_id(
5057  disambiguation_function(manifold_ids[id]));
5058  cell->line(l)->set_user_index(0);
5059  }
5060  }
5061  if (dim > 2)
5062  for (unsigned int l = 0; l < cell->n_faces(); ++l)
5063  {
5064  const auto id = cell->quad(l)->user_index();
5065  // Make sure we change the manifold indicator only once
5066  if (id != 0)
5067  {
5068  if (cell->quad(l)->manifold_id() ==
5070  overwrite_only_flat_manifold_ids == false)
5071  cell->quad(l)->set_manifold_id(
5072  disambiguation_function(manifold_ids[id]));
5073  cell->quad(l)->set_user_index(0);
5074  }
5075  }
5076  }
5077  tria.load_user_indices(backup);
5078  }
5079 
5080 
5081 
5082  template <int dim, int spacedim>
5083  std::pair<unsigned int, double>
5086  {
5087  double max_ratio = 1;
5088  unsigned int index = 0;
5089 
5090  for (unsigned int i = 0; i < dim; ++i)
5091  for (unsigned int j = i + 1; j < dim; ++j)
5092  {
5093  unsigned int ax = i % dim;
5094  unsigned int next_ax = j % dim;
5095 
5096  double ratio =
5097  cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
5098 
5099  if (ratio > max_ratio)
5100  {
5101  max_ratio = ratio;
5102  index = ax;
5103  }
5104  else if (1.0 / ratio > max_ratio)
5105  {
5106  max_ratio = 1.0 / ratio;
5107  index = next_ax;
5108  }
5109  }
5110  return std::make_pair(index, max_ratio);
5111  }
5112 
5113 
5114  template <int dim, int spacedim>
5115  void
5117  const bool isotropic,
5118  const unsigned int max_iterations)
5119  {
5120  unsigned int iter = 0;
5121  bool continue_refinement = true;
5122 
5123  while (continue_refinement && (iter < max_iterations))
5124  {
5125  if (max_iterations != numbers::invalid_unsigned_int)
5126  iter++;
5127  continue_refinement = false;
5128 
5129  for (const auto &cell : tria.active_cell_iterators())
5130  for (const unsigned int j : cell->face_indices())
5131  if (cell->at_boundary(j) == false &&
5132  cell->neighbor(j)->has_children())
5133  {
5134  if (isotropic)
5135  {
5136  cell->set_refine_flag();
5137  continue_refinement = true;
5138  }
5139  else
5140  continue_refinement |= cell->flag_for_face_refinement(j);
5141  }
5142 
5144  }
5145  }
5146 
5147  template <int dim, int spacedim>
5148  void
5150  const double max_ratio,
5151  const unsigned int max_iterations)
5152  {
5153  unsigned int iter = 0;
5154  bool continue_refinement = true;
5155 
5156  while (continue_refinement && (iter < max_iterations))
5157  {
5158  iter++;
5159  continue_refinement = false;
5160  for (const auto &cell : tria.active_cell_iterators())
5161  {
5162  std::pair<unsigned int, double> info =
5163  GridTools::get_longest_direction<dim, spacedim>(cell);
5164  if (info.second > max_ratio)
5165  {
5166  cell->set_refine_flag(
5167  RefinementCase<dim>::cut_axis(info.first));
5168  continue_refinement = true;
5169  }
5170  }
5172  }
5173  }
5174 
5175 
5176  template <int dim, int spacedim>
5177  void
5179  const double limit_angle_fraction)
5180  {
5181  if (dim == 1)
5182  return; // Nothing to do
5183 
5184  // Check that we don't have hanging nodes
5186  ExcMessage("The input Triangulation cannot "
5187  "have hanging nodes."));
5188 
5189 
5190  bool has_cells_with_more_than_dim_faces_on_boundary = true;
5191  bool has_cells_with_dim_faces_on_boundary = false;
5192 
5193  unsigned int refinement_cycles = 0;
5194 
5195  while (has_cells_with_more_than_dim_faces_on_boundary)
5196  {
5197  has_cells_with_more_than_dim_faces_on_boundary = false;
5198 
5199  for (const auto &cell : tria.active_cell_iterators())
5200  {
5201  unsigned int boundary_face_counter = 0;
5202  for (auto f : cell->face_indices())
5203  if (cell->face(f)->at_boundary())
5204  boundary_face_counter++;
5205  if (boundary_face_counter > dim)
5206  {
5207  has_cells_with_more_than_dim_faces_on_boundary = true;
5208  break;
5209  }
5210  else if (boundary_face_counter == dim)
5211  has_cells_with_dim_faces_on_boundary = true;
5212  }
5213  if (has_cells_with_more_than_dim_faces_on_boundary)
5214  {
5215  tria.refine_global(1);
5216  refinement_cycles++;
5217  }
5218  }
5219 
5220  if (has_cells_with_dim_faces_on_boundary)
5221  {
5222  tria.refine_global(1);
5223  refinement_cycles++;
5224  }
5225  else
5226  {
5227  while (refinement_cycles > 0)
5228  {
5229  for (const auto &cell : tria.active_cell_iterators())
5230  cell->set_coarsen_flag();
5232  refinement_cycles--;
5233  }
5234  return;
5235  }
5236 
5237  std::vector<bool> cells_to_remove(tria.n_active_cells(), false);
5238  std::vector<Point<spacedim>> vertices = tria.get_vertices();
5239 
5240  std::vector<bool> faces_to_remove(tria.n_raw_faces(), false);
5241 
5242  std::vector<CellData<dim>> cells_to_add;
5243  SubCellData subcelldata_to_add;
5244 
5245  // Trick compiler for dimension independent things
5246  const unsigned int v0 = 0, v1 = 1, v2 = (dim > 1 ? 2 : 0),
5247  v3 = (dim > 1 ? 3 : 0);
5248 
5249  for (const auto &cell : tria.active_cell_iterators())
5250  {
5251  double angle_fraction = 0;
5252  unsigned int vertex_at_corner = numbers::invalid_unsigned_int;
5253 
5254  if (dim == 2)
5255  {
5257  p0[spacedim > 1 ? 1 : 0] = 1;
5259  p1[0] = 1;
5260 
5261  if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
5262  {
5263  p0 = cell->vertex(v0) - cell->vertex(v2);
5264  p1 = cell->vertex(v3) - cell->vertex(v2);
5265  vertex_at_corner = v2;
5266  }
5267  else if (cell->face(v3)->at_boundary() &&
5268  cell->face(v1)->at_boundary())
5269  {
5270  p0 = cell->vertex(v2) - cell->vertex(v3);
5271  p1 = cell->vertex(v1) - cell->vertex(v3);
5272  vertex_at_corner = v3;
5273  }
5274  else if (cell->face(1)->at_boundary() &&
5275  cell->face(2)->at_boundary())
5276  {
5277  p0 = cell->vertex(v0) - cell->vertex(v1);
5278  p1 = cell->vertex(v3) - cell->vertex(v1);
5279  vertex_at_corner = v1;
5280  }
5281  else if (cell->face(2)->at_boundary() &&
5282  cell->face(0)->at_boundary())
5283  {
5284  p0 = cell->vertex(v2) - cell->vertex(v0);
5285  p1 = cell->vertex(v1) - cell->vertex(v0);
5286  vertex_at_corner = v0;
5287  }
5288  p0 /= p0.norm();
5289  p1 /= p1.norm();
5290  angle_fraction = std::acos(p0 * p1) / numbers::PI;
5291  }
5292  else
5293  {
5294  Assert(false, ExcNotImplemented());
5295  }
5296 
5297  if (angle_fraction > limit_angle_fraction)
5298  {
5299  auto flags_removal = [&](unsigned int f1,
5300  unsigned int f2,
5301  unsigned int n1,
5302  unsigned int n2) -> void {
5303  cells_to_remove[cell->active_cell_index()] = true;
5304  cells_to_remove[cell->neighbor(n1)->active_cell_index()] = true;
5305  cells_to_remove[cell->neighbor(n2)->active_cell_index()] = true;
5306 
5307  faces_to_remove[cell->face(f1)->index()] = true;
5308  faces_to_remove[cell->face(f2)->index()] = true;
5309 
5310  faces_to_remove[cell->neighbor(n1)->face(f1)->index()] = true;
5311  faces_to_remove[cell->neighbor(n2)->face(f2)->index()] = true;
5312  };
5313 
5314  auto cell_creation = [&](const unsigned int vv0,
5315  const unsigned int vv1,
5316  const unsigned int f0,
5317  const unsigned int f1,
5318 
5319  const unsigned int n0,
5320  const unsigned int v0n0,
5321  const unsigned int v1n0,
5322 
5323  const unsigned int n1,
5324  const unsigned int v0n1,
5325  const unsigned int v1n1) {
5326  CellData<dim> c1, c2;
5327  CellData<1> l1, l2;
5328 
5329  c1.vertices[v0] = cell->vertex_index(vv0);
5330  c1.vertices[v1] = cell->vertex_index(vv1);
5331  c1.vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
5332  c1.vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
5333 
5334  c1.manifold_id = cell->manifold_id();
5335  c1.material_id = cell->material_id();
5336 
5337  c2.vertices[v0] = cell->vertex_index(vv0);
5338  c2.vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
5339  c2.vertices[v2] = cell->vertex_index(vv1);
5340  c2.vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
5341 
5342  c2.manifold_id = cell->manifold_id();
5343  c2.material_id = cell->material_id();
5344 
5345  l1.vertices[0] = cell->vertex_index(vv0);
5346  l1.vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
5347 
5348  l1.boundary_id = cell->line(f0)->boundary_id();
5349  l1.manifold_id = cell->line(f0)->manifold_id();
5350  subcelldata_to_add.boundary_lines.push_back(l1);
5351 
5352  l2.vertices[0] = cell->vertex_index(vv0);
5353  l2.vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
5354 
5355  l2.boundary_id = cell->line(f1)->boundary_id();
5356  l2.manifold_id = cell->line(f1)->manifold_id();
5357  subcelldata_to_add.boundary_lines.push_back(l2);
5358 
5359  cells_to_add.push_back(c1);
5360  cells_to_add.push_back(c2);
5361  };
5362 
5363  if (dim == 2)
5364  {
5365  switch (vertex_at_corner)
5366  {
5367  case 0:
5368  flags_removal(0, 2, 3, 1);
5369  cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
5370  break;
5371  case 1:
5372  flags_removal(1, 2, 3, 0);
5373  cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
5374  break;
5375  case 2:
5376  flags_removal(3, 0, 1, 2);
5377  cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
5378  break;
5379  case 3:
5380  flags_removal(3, 1, 0, 2);
5381  cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
5382  break;
5383  }
5384  }
5385  else
5386  {
5387  Assert(false, ExcNotImplemented());
5388  }
5389  }
5390  }
5391 
5392  // if no cells need to be added, then no regularization is necessary.
5393  // Restore things as they were before this function was called.
5394  if (cells_to_add.size() == 0)
5395  {
5396  while (refinement_cycles > 0)
5397  {
5398  for (const auto &cell : tria.active_cell_iterators())
5399  cell->set_coarsen_flag();
5401  refinement_cycles--;
5402  }
5403  return;
5404  }
5405 
5406  // add the cells that were not marked as skipped
5407  for (const auto &cell : tria.active_cell_iterators())
5408  {
5409  if (cells_to_remove[cell->active_cell_index()] == false)
5410  {
5411  CellData<dim> c;
5412  for (const unsigned int v : cell->vertex_indices())
5413  c.vertices[v] = cell->vertex_index(v);
5414  c.manifold_id = cell->manifold_id();
5415  c.material_id = cell->material_id();
5416  cells_to_add.push_back(c);
5417  }
5418  }
5419 
5420  // Face counter for both dim == 2 and dim == 3
5422  face = tria.begin_active_face(),
5423  endf = tria.end_face();
5424  for (; face != endf; ++face)
5425  if ((face->at_boundary() ||
5426  face->manifold_id() != numbers::flat_manifold_id) &&
5427  faces_to_remove[face->index()] == false)
5428  {
5429  for (unsigned int l = 0; l < face->n_lines(); ++l)
5430  {
5431  CellData<1> line;
5432  if (dim == 2)
5433  {
5434  for (const unsigned int v : face->vertex_indices())
5435  line.vertices[v] = face->vertex_index(v);
5436  line.boundary_id = face->boundary_id();
5437  line.manifold_id = face->manifold_id();
5438  }
5439  else
5440  {
5441  for (const unsigned int v : face->line(l)->vertex_indices())
5442  line.vertices[v] = face->line(l)->vertex_index(v);
5443  line.boundary_id = face->line(l)->boundary_id();
5444  line.manifold_id = face->line(l)->manifold_id();
5445  }
5446  subcelldata_to_add.boundary_lines.push_back(line);
5447  }
5448  if (dim == 3)
5449  {
5450  CellData<2> quad;
5451  for (const unsigned int v : face->vertex_indices())
5452  quad.vertices[v] = face->vertex_index(v);
5453  quad.boundary_id = face->boundary_id();
5454  quad.manifold_id = face->manifold_id();
5455  subcelldata_to_add.boundary_quads.push_back(quad);
5456  }
5457  }
5459  cells_to_add,
5460  subcelldata_to_add);
5462 
5463  // Save manifolds
5464  auto manifold_ids = tria.get_manifold_ids();
5465  std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
5466  manifolds;
5467  // Set manifolds in new Triangulation
5468  for (const auto manifold_id : manifold_ids)
5470  manifolds[manifold_id] = tria.get_manifold(manifold_id).clone();
5471 
5472  tria.clear();
5473 
5474  tria.create_triangulation(vertices, cells_to_add, subcelldata_to_add);
5475 
5476  // Restore manifolds
5477  for (const auto manifold_id : manifold_ids)
5479  tria.set_manifold(manifold_id, *manifolds[manifold_id]);
5480  }
5481 
5482 
5483 
5484  template <int dim, int spacedim>
5485 #ifndef DOXYGEN
5486  std::tuple<
5487  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5488  std::vector<std::vector<Point<dim>>>,
5489  std::vector<std::vector<unsigned int>>>
5490 #else
5491  return_type
5492 #endif
5494  const Cache<dim, spacedim> & cache,
5495  const std::vector<Point<spacedim>> &points,
5497  &cell_hint)
5498  {
5499  const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
5500  // Splitting the tuple's components
5501  auto &cells = std::get<0>(cqmp);
5502  auto &qpoints = std::get<1>(cqmp);
5503  auto &maps = std::get<2>(cqmp);
5504 
5505  return std::make_tuple(std::move(cells),
5506  std::move(qpoints),
5507  std::move(maps));
5508  }
5509 
5510 
5511 
5512  template <int dim, int spacedim>
5513 #ifndef DOXYGEN
5514  std::tuple<
5515  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5516  std::vector<std::vector<Point<dim>>>,
5517  std::vector<std::vector<unsigned int>>,
5518  std::vector<unsigned int>>
5519 #else
5520  return_type
5521 #endif
5523  const Cache<dim, spacedim> & cache,
5524  const std::vector<Point<spacedim>> &points,
5526  &cell_hint)
5527  {
5528  Assert((dim == spacedim),
5529  ExcMessage("Only implemented for dim==spacedim."));
5530 
5531  // Alias
5532  namespace bgi = boost::geometry::index;
5533 
5534  // Get the mapping
5535  const auto &mapping = cache.get_mapping();
5536 
5537  // How many points are here?
5538  const unsigned int np = points.size();
5539 
5540  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5541  cells_out;
5542  std::vector<std::vector<Point<dim>>> qpoints_out;
5543  std::vector<std::vector<unsigned int>> maps_out;
5544  std::vector<unsigned int> missing_points_out;
5545 
5546  // Now the easy case.
5547  if (np == 0)
5548  return std::make_tuple(std::move(cells_out),
5549  std::move(qpoints_out),
5550  std::move(maps_out),
5551  std::move(missing_points_out));
5552 
5553  // For the search we shall use the following tree
5554  const auto &b_tree = cache.get_cell_bounding_boxes_rtree();
5555 
5556  // Now make a tree of indices for the points
5557  // [TODO] This would work better with pack_rtree_of_indices, but
5558  // windows does not like it. Build a tree with pairs of point and id
5559  std::vector<std::pair<Point<spacedim>, unsigned int>> points_and_ids(np);
5560  for (unsigned int i = 0; i < np; ++i)
5561  points_and_ids[i] = std::make_pair(points[i], i);
5562  const auto p_tree = pack_rtree(points_and_ids);
5563 
5564  // Keep track of all found points
5565  std::vector<bool> found_points(points.size(), false);
5566 
5567  // Check if a point was found
5568  const auto already_found = [&found_points](const auto &id) {
5569  AssertIndexRange(id.second, found_points.size());
5570  return found_points[id.second];
5571  };
5572 
5573  // check if the given cell was already in the vector of cells before. If so,
5574  // insert in the corresponding vectors the reference point and the id.
5575  // Otherwise append a new entry to all vectors.
5576  const auto store_cell_point_and_id =
5577  [&](
5579  const Point<dim> & ref_point,
5580  const unsigned int &id) {
5581  const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
5582  if (it != cells_out.rend())
5583  {
5584  const auto cell_id =
5585  (cells_out.size() - 1 - (it - cells_out.rbegin()));
5586  qpoints_out[cell_id].emplace_back(ref_point);
5587  maps_out[cell_id].emplace_back(id);
5588  }
5589  else
5590  {
5591  cells_out.emplace_back(cell);
5592  qpoints_out.emplace_back(std::vector<Point<dim>>({ref_point}));
5593  maps_out.emplace_back(std::vector<unsigned int>({id}));
5594  }
5595  };
5596 
5597  // Check all points within a given pair of box and cell
5598  const auto check_all_points_within_box = [&](const auto &leaf) {
5599  const auto &box = leaf.first;
5600  const auto &cell_hint = leaf.second;
5601 
5602  for (const auto &point_and_id :
5603  p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
5604  bgi::intersects(box)))
5605  {
5606  const auto id = point_and_id.second;
5607  const auto cell_and_ref =
5609  points[id],
5610  cell_hint);
5611  const auto &cell = cell_and_ref.first;
5612  const auto &ref_point = cell_and_ref.second;
5613 
5614  if (cell.state() == IteratorState::valid)
5615  store_cell_point_and_id(cell, ref_point, id);
5616  else
5617  missing_points_out.emplace_back(id);
5618 
5619  // Don't look anymore for this point
5620  found_points[id] = true;
5621  }
5622  };
5623 
5624  // If a hint cell was given, use it
5625  if (cell_hint.state() == IteratorState::valid)
5626  check_all_points_within_box(
5627  std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
5628 
5629  // Now loop over all points that have not been found yet
5630  for (unsigned int i = 0; i < np; ++i)
5631  if (found_points[i] == false)
5632  {
5633  // Get the closest cell to this point
5634  const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
5635  // Now checks all points that fall within this box
5636  if (leaf != b_tree.qend())
5637  check_all_points_within_box(*leaf);
5638  else
5639  {
5640  // We should not get here. Throw an error.
5641  Assert(false, ExcInternalError());
5642  }
5643  }
5644  // Now make sure we send out the rest of the points that we did not find.
5645  for (unsigned int i = 0; i < np; ++i)
5646  if (found_points[i] == false)
5647  missing_points_out.emplace_back(i);
5648 
5649  // Debug Checking
5650  AssertDimension(cells_out.size(), maps_out.size());
5651  AssertDimension(cells_out.size(), qpoints_out.size());
5652 
5653 #ifdef DEBUG
5654  unsigned int c = cells_out.size();
5655  unsigned int qps = 0;
5656  // The number of points in all
5657  // the cells must be the same as
5658  // the number of points we
5659  // started off from,
5660  // plus the points which were ignored
5661  for (unsigned int n = 0; n < c; ++n)
5662  {
5663  AssertDimension(qpoints_out[n].size(), maps_out[n].size());
5664  qps += qpoints_out[n].size();
5665  }
5666 
5667  Assert(qps + missing_points_out.size() == np,
5668  ExcDimensionMismatch(qps + missing_points_out.size(), np));
5669 #endif
5670 
5671  return std::make_tuple(std::move(cells_out),
5672  std::move(qpoints_out),
5673  std::move(maps_out),
5674  std::move(missing_points_out));
5675  }
5676 
5677 
5678 
5679  template <int dim, int spacedim>
5680 #ifndef DOXYGEN
5681  std::tuple<
5682  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5683  std::vector<std::vector<Point<dim>>>,
5684  std::vector<std::vector<unsigned int>>,
5685  std::vector<std::vector<Point<spacedim>>>,
5686  std::vector<std::vector<unsigned int>>>
5687 #else
5688  return_type
5689 #endif
5691  const GridTools::Cache<dim, spacedim> & cache,
5692  const std::vector<Point<spacedim>> & points,
5693  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5694  const double tolerance)
5695  {
5696  // run internal function ...
5698  cache, points, global_bboxes, tolerance, false, true)
5699  .send_components;
5700 
5701  // ... and reshuffle the data
5702  std::tuple<
5703  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5704  std::vector<std::vector<Point<dim>>>,
5705  std::vector<std::vector<unsigned int>>,
5706  std::vector<std::vector<Point<spacedim>>>,
5707  std::vector<std::vector<unsigned int>>>
5708  result;
5709 
5710  std::pair<int, int> dummy{-1, -1};
5711 
5712  for (unsigned int i = 0; i < all.size(); ++i)
5713  {
5714  if (dummy != std::get<0>(all[i]))
5715  {
5716  std::get<0>(result).push_back(
5718  &cache.get_triangulation(),
5719  std::get<0>(all[i]).first,
5720  std::get<0>(all[i]).second});
5721 
5722  const unsigned int new_size = std::get<0>(result).size();
5723 
5724  std::get<1>(result).resize(new_size);
5725  std::get<2>(result).resize(new_size);
5726  std::get<3>(result).resize(new_size);
5727  std::get<4>(result).resize(new_size);
5728 
5729  dummy = std::get<0>(all[i]);
5730  }
5731 
5732  std::get<1>(result).back().push_back(
5733  std::get<3>(all[i])); // reference point
5734  std::get<2>(result).back().push_back(std::get<2>(all[i])); // index
5735  std::get<3>(result).back().push_back(std::get<4>(all[i])); // real point
5736  std::get<4>(result).back().push_back(std::get<1>(all[i])); // rank
5737  }
5738 
5739  return result;
5740  }
5741 
5742 
5743 
5744  namespace internal
5745  {
5746  template <int spacedim>
5747  std::tuple<std::vector<unsigned int>,
5748  std::vector<unsigned int>,
5749  std::vector<unsigned int>>
5751  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5752  const std::vector<Point<spacedim>> & points,
5753  const double tolerance)
5754  {
5755  std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
5756  ranks_and_indices.reserve(points.size());
5757 
5758  for (unsigned int i = 0; i < points.size(); ++i)
5759  {
5760  const auto &point = points[i];
5761  for (unsigned rank = 0; rank < global_bboxes.size(); ++rank)
5762  for (const auto &box : global_bboxes[rank])
5763  if (box.point_inside(point, tolerance))
5764  {
5765  ranks_and_indices.emplace_back(rank, i);
5766  break;
5767  }
5768  }
5769 
5770  // convert to CRS
5771  std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
5772 
5773  std::vector<unsigned int> ranks;
5774  std::vector<unsigned int> ptr;
5775  std::vector<unsigned int> indices;
5776 
5777  unsigned int dummy_rank = numbers::invalid_unsigned_int;
5778 
5779  for (const auto &i : ranks_and_indices)
5780  {
5781  if (dummy_rank != i.first)
5782  {
5783  dummy_rank = i.first;
5784  ranks.push_back(dummy_rank);
5785  ptr.push_back(indices.size());
5786  }
5787 
5788  indices.push_back(i.second);
5789  }
5790  ptr.push_back(indices.size());
5791 
5792  return std::make_tuple(std::move(ranks),
5793  std::move(ptr),
5794  std::move(indices));
5795  }
5796 
5797 
5798 
5799  template <int dim, int spacedim>
5800  std::vector<
5801  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5802  Point<dim>>>
5804  const Cache<dim, spacedim> & cache,
5805  const Point<spacedim> & point,
5807  const std::vector<bool> &marked_vertices,
5808  const double tolerance,
5809  const bool enforce_unique_mapping)
5810  {
5811  std::vector<
5812  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5813  Point<dim>>>
5814  locally_owned_active_cells_around_point;
5815 
5816  const auto first_cell = GridTools::find_active_cell_around_point(
5817  cache.get_mapping(),
5818  cache.get_triangulation(),
5819  point,
5820  cache.get_vertex_to_cell_map(),
5822  cell_hint,
5823  marked_vertices,
5824  cache.get_used_vertices_rtree(),
5825  tolerance,
5827 
5828  const unsigned int my_rank = Utilities::MPI::this_mpi_process(
5829  cache.get_triangulation().get_communicator());
5830 
5831  cell_hint = first_cell.first;
5832  if (cell_hint.state() == IteratorState::valid)
5833  {
5834  const auto active_cells_around_point =
5836  cache.get_mapping(),
5837  cache.get_triangulation(),
5838  point,
5839  tolerance,
5840  first_cell);
5841 
5842  if (enforce_unique_mapping)
5843  {
5844  // check if the rank of this process is the lowest of all cells
5845  // if not, the other process will handle this cell and we don't
5846  // have to do here anything in the case of unique mapping
5847  unsigned int lowes_rank = numbers::invalid_unsigned_int;
5848 
5849  for (const auto &cell : active_cells_around_point)
5850  lowes_rank = std::min(lowes_rank, cell.first->subdomain_id());
5851 
5852  if (lowes_rank != my_rank)
5853  return {};
5854  }
5855 
5856  locally_owned_active_cells_around_point.reserve(
5857  active_cells_around_point.size());
5858 
5859  for (const auto &cell : active_cells_around_point)
5860  if (cell.first->is_locally_owned())
5861  locally_owned_active_cells_around_point.push_back(cell);
5862  }
5863 
5864  std::sort(locally_owned_active_cells_around_point.begin(),
5865  locally_owned_active_cells_around_point.end(),
5866  [](const auto &a, const auto &b) { return a.first < b.first; });
5867 
5868  if (enforce_unique_mapping &&
5869  locally_owned_active_cells_around_point.size() > 1)
5870  // in the case of unique mapping, we only need a single cell
5871  return {locally_owned_active_cells_around_point.front()};
5872  else
5873  return locally_owned_active_cells_around_point;
5874  }
5875 
5876 
5877 
5878  template <int dim, int spacedim>
5881  const GridTools::Cache<dim, spacedim> & cache,
5882  const std::vector<Point<spacedim>> & points,
5883  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5884  const double tolerance,
5885  const bool perform_handshake,
5886  const bool enforce_unique_mapping)
5887  {
5889 
5890  auto &send_components = result.send_components;
5891  auto &send_ranks = result.send_ranks;
5892  auto &send_ptrs = result.send_ptrs;
5893  auto &recv_components = result.recv_components;
5894  auto &recv_ranks = result.recv_ranks;
5895  auto &recv_ptrs = result.recv_ptrs;
5896 
5897  const auto potential_owners =
5898  internal::guess_point_owner(global_bboxes, points, tolerance);
5899 
5900  const auto &potential_owners_ranks = std::get<0>(potential_owners);
5901  const auto &potential_owners_ptrs = std::get<1>(potential_owners);
5902  const auto &potential_owners_indices = std::get<2>(potential_owners);
5903 
5904  const std::vector<bool> marked_vertices;
5905  auto cell_hint = cache.get_triangulation().begin_active();
5906 
5907  const auto translate = [&](const unsigned int other_rank) {
5908  const auto ptr = std::find(potential_owners_ranks.begin(),
5909  potential_owners_ranks.end(),
5910  other_rank);
5911 
5912  Assert(ptr != potential_owners_ranks.end(), ExcInternalError());
5913 
5914  const auto other_rank_index =
5915  std::distance(potential_owners_ranks.begin(), ptr);
5916 
5917  return other_rank_index;
5918  };
5919 
5921  [&]() { return potential_owners_ranks; },
5922  [&](const unsigned int other_rank, std::vector<char> &send_buffer) {
5923  const auto other_rank_index = translate(other_rank);
5924 
5925  std::vector<std::pair<unsigned int, Point<spacedim>>> temp;
5926  temp.reserve(potential_owners_ptrs[other_rank_index + 1] -
5927  potential_owners_ptrs[other_rank_index]);
5928 
5929  for (unsigned int i = potential_owners_ptrs[other_rank_index];
5930  i < potential_owners_ptrs[other_rank_index + 1];
5931  ++i)
5932  temp.emplace_back(potential_owners_indices[i],
5933  points[potential_owners_indices[i]]);
5934 
5935  send_buffer = Utilities::pack(temp, false);
5936  },
5937  [&](const unsigned int & other_rank,
5938  const std::vector<char> &recv_buffer,
5939  std::vector<char> & request_buffer) {
5940  const auto recv_buffer_unpacked = Utilities::unpack<
5941  std::vector<std::pair<unsigned int, Point<spacedim>>>>(recv_buffer,
5942  false);
5943 
5944  std::vector<unsigned int> request_buffer_temp(
5945  recv_buffer_unpacked.size(), 0);
5946 
5947  cell_hint = cache.get_triangulation().begin_active();
5948 
5949  for (unsigned int i = 0; i < recv_buffer_unpacked.size(); ++i)
5950  {
5951  const auto &index_and_point = recv_buffer_unpacked[i];
5952 
5953  const auto cells_and_reference_positions =
5955  cache,
5956  index_and_point.second,
5957  cell_hint,
5958  marked_vertices,
5959  tolerance,
5960  enforce_unique_mapping);
5961 
5962  for (const auto &cell_and_reference_position :
5963  cells_and_reference_positions)
5964  {
5965  send_components.emplace_back(
5966  std::pair<int, int>(
5967  cell_and_reference_position.first->level(),
5968  cell_and_reference_position.first->index()),
5969  other_rank,
5970  index_and_point.first,
5971  cell_and_reference_position.second,
5972  index_and_point.second,
5974  }
5975 
5976  request_buffer_temp[i] = cells_and_reference_positions.size();
5977  }
5978 
5979  if (perform_handshake)
5980  request_buffer = Utilities::pack(request_buffer_temp, false);
5981  },
5982  [&](const unsigned int other_rank, std::vector<char> &recv_buffer) {
5983  if (perform_handshake)
5984  {
5985  const auto other_rank_index = translate(other_rank);
5986 
5987  recv_buffer =
5988  Utilities::pack(std::vector<unsigned int>(
5989  potential_owners_ptrs[other_rank_index + 1] -
5990  potential_owners_ptrs[other_rank_index]),
5991  false);
5992  }
5993  },
5994  [&](const unsigned int other_rank,
5995  const std::vector<char> &recv_buffer) {
5996  if (perform_handshake)
5997  {
5998  const auto recv_buffer_unpacked =
5999  Utilities::unpack<std::vector<unsigned int>>(recv_buffer,
6000  false);
6001 
6002  const auto other_rank_index = translate(other_rank);
6003 
6004  for (unsigned int i = 0; i < recv_buffer_unpacked.size(); ++i)
6005  for (unsigned int j = 0; j < recv_buffer_unpacked[i]; ++j)
6006  recv_components.emplace_back(
6007  other_rank,
6008  potential_owners_indices
6009  [i + potential_owners_ptrs[other_rank_index]],
6011  }
6012  });
6013 
6015  process, cache.get_triangulation().get_communicator())
6016  .run();
6017 
6018  if (true)
6019  {
6020  // sort according to rank (and point index and cell) -> make
6021  // deterministic
6022  std::sort(send_components.begin(),
6023  send_components.end(),
6024  [&](const auto &a, const auto &b) {
6025  if (std::get<1>(a) != std::get<1>(b)) // rank
6026  return std::get<1>(a) < std::get<1>(b);
6027 
6028  if (std::get<2>(a) != std::get<2>(b)) // point index
6029  return std::get<2>(a) < std::get<2>(b);
6030 
6031  return std::get<0>(a) < std::get<0>(b); // cell
6032  });
6033 
6034  // perform enumeration and extract rank information
6035  for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
6036  i < send_components.size();
6037  ++i)
6038  {
6039  std::get<5>(send_components[i]) = i;
6040 
6041  if (dummy != std::get<1>(send_components[i]))
6042  {
6043  dummy = std::get<1>(send_components[i]);
6044  send_ranks.push_back(dummy);
6045  send_ptrs.push_back(i);
6046  }
6047  }
6048  send_ptrs.push_back(send_components.size());
6049 
6050  // sort according to cell, rank, point index (while keeping
6051  // partial ordering)
6052  std::sort(send_components.begin(),
6053  send_components.end(),
6054  [&](const auto &a, const auto &b) {
6055  if (std::get<0>(a) != std::get<0>(b))
6056  return std::get<0>(a) < std::get<0>(b); // cell
6057 
6058  if (std::get<1>(a) != std::get<1>(b))
6059  return std::get<1>(a) < std::get<1>(b); // rank
6060 
6061  if (std::get<2>(a) != std::get<2>(b))
6062  return std::get<2>(a) < std::get<2>(b); // point index
6063 
6064  return std::get<5>(a) < std::get<5>(b); // enumeration
6065  });
6066  }
6067 
6068  if (perform_handshake)
6069  {
6070  // sort according to rank (and point index) -> make deterministic
6071  std::sort(recv_components.begin(),
6072  recv_components.end(),
6073  [&](const auto &a, const auto &b) {
6074  if (std::get<0>(a) != std::get<0>(b))
6075  return std::get<0>(a) < std::get<0>(b); // rank
6076 
6077  return std::get<1>(a) < std::get<1>(b); // point index
6078  });
6079 
6080  // perform enumeration and extract rank information
6081  for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
6082  i < recv_components.size();
6083  ++i)
6084  {
6085  std::get<2>(recv_components[i]) = i;
6086 
6087  if (dummy != std::get<0>(recv_components[i]))
6088  {
6089  dummy = std::get<0>(recv_components[i]);
6090  recv_ranks.push_back(dummy);
6091  recv_ptrs.push_back(i);
6092  }
6093  }
6094  recv_ptrs.push_back(recv_components.size());
6095 
6096  // sort according to point index and rank (while keeping partial
6097  // ordering)
6098  std::sort(recv_components.begin(),
6099  recv_components.end(),
6100  [&](const auto &a, const auto &b) {
6101  if (std::get<1>(a) != std::get<1>(b))
6102  return std::get<1>(a) < std::get<1>(b); // point index
6103 
6104  if (std::get<0>(a) != std::get<0>(b))
6105  return std::get<0>(a) < std::get<0>(b); // rank
6106 
6107  return std::get<2>(a) < std::get<2>(b); // enumeration
6108  });
6109  }
6110 
6111  return result;
6112  }
6113  } // namespace internal
6114 
6115 
6116 
6117  template <int dim, int spacedim>
6118  std::map<unsigned int, Point<spacedim>>
6120  const Mapping<dim, spacedim> & mapping)
6121  {
6122  std::map<unsigned int, Point<spacedim>> result;
6123  for (const auto &cell : container.active_cell_iterators())
6124  {
6125  if (!cell->is_artificial())
6126  {
6127  const auto vs = mapping.get_vertices(cell);
6128  for (unsigned int i = 0; i < vs.size(); ++i)
6129  result[cell->vertex_index(i)] = vs[i];
6130  }
6131  }
6132  return result;
6133  }
6134 
6135 
6136  template <int spacedim>
6137  unsigned int
6138  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
6139  const Point<spacedim> & p)
6140  {
6141  auto id_and_v = std::min_element(
6142  vertices.begin(),
6143  vertices.end(),
6144  [&](const std::pair<const unsigned int, Point<spacedim>> &p1,
6145  const std::pair<const unsigned int, Point<spacedim>> &p2) -> bool {
6146  return p1.second.distance(p) < p2.second.distance(p);
6147  });
6148  return id_and_v->first;
6149  }
6150 
6151 
6152  template <int dim, int spacedim>
6153  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
6154  Point<dim>>
6156  const Cache<dim, spacedim> &cache,
6157  const Point<spacedim> & p,
6159  & cell_hint,
6160  const std::vector<bool> &marked_vertices,
6161  const double tolerance)
6162  {
6163  const auto &mesh = cache.get_triangulation();
6164  const auto &mapping = cache.get_mapping();
6165  const auto &vertex_to_cells = cache.get_vertex_to_cell_map();
6166  const auto &vertex_to_cell_centers =
6168  const auto &used_vertices_rtree = cache.get_used_vertices_rtree();
6169 
6170  return find_active_cell_around_point(mapping,
6171  mesh,
6172  p,
6173  vertex_to_cells,
6174  vertex_to_cell_centers,
6175  cell_hint,
6176  marked_vertices,
6177  used_vertices_rtree,
6178  tolerance);
6179  }
6180 
6181  template <int spacedim>
6182  std::vector<std::vector<BoundingBox<spacedim>>>
6184  const std::vector<BoundingBox<spacedim>> &local_bboxes,
6185  const MPI_Comm & mpi_communicator)
6186  {
6187 #ifndef DEAL_II_WITH_MPI
6188  (void)local_bboxes;
6189  (void)mpi_communicator;
6190  Assert(false,
6191  ExcMessage(
6192  "GridTools::exchange_local_bounding_boxes() requires MPI."));
6193  return {};
6194 #else
6195  // Step 1: preparing data to be sent
6196  unsigned int n_bboxes = local_bboxes.size();
6197  // Dimension of the array to be exchanged (number of double)
6198  int n_local_data = 2 * spacedim * n_bboxes;
6199  // data array stores each entry of each point describing the bounding boxes
6200  std::vector<double> loc_data_array(n_local_data);
6201  for (unsigned int i = 0; i < n_bboxes; ++i)
6202  for (unsigned int d = 0; d < spacedim; ++d)
6203  {
6204  // Extracting the coordinates of each boundary point
6205  loc_data_array[2 * i * spacedim + d] =
6206  local_bboxes[i].get_boundary_points().first[d];
6207  loc_data_array[2 * i * spacedim + spacedim + d] =
6208  local_bboxes[i].get_boundary_points().second[d];
6209  }
6210 
6211  // Step 2: exchanging the size of local data
6212  unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_communicator);
6213 
6214  // Vector to store the size of loc_data_array for every process
6215  std::vector<int> size_all_data(n_procs);
6216 
6217  // Exchanging the number of bboxes
6218  int ierr = MPI_Allgather(&n_local_data,
6219  1,
6220  MPI_INT,
6221  size_all_data.data(),
6222  1,
6223  MPI_INT,
6224  mpi_communicator);
6225  AssertThrowMPI(ierr);
6226 
6227  // Now computing the the displacement, relative to recvbuf,
6228  // at which to store the incoming data
6229  std::vector<int> rdispls(n_procs);
6230  rdispls[0] = 0;
6231  for (unsigned int i = 1; i < n_procs; ++i)
6232  rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
6233 
6234  // Step 3: exchange the data and bounding boxes:
6235  // Allocating a vector to contain all the received data
6236  std::vector<double> data_array(rdispls.back() + size_all_data.back());
6237 
6238  ierr = MPI_Allgatherv(loc_data_array.data(),
6239  n_local_data,
6240  MPI_DOUBLE,
6241  data_array.data(),
6242  size_all_data.data(),
6243  rdispls.data(),
6244  MPI_DOUBLE,
6245  mpi_communicator);
6246  AssertThrowMPI(ierr);
6247 
6248  // Step 4: create the array of bboxes for output
6249  std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
6250  unsigned int begin_idx = 0;
6251  for (unsigned int i = 0; i < n_procs; ++i)
6252  {
6253  // Number of local bounding boxes
6254  unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
6255  global_bboxes[i].resize(n_bbox_i);
6256  for (unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
6257  {
6258  Point<spacedim> p1, p2; // boundary points for bbox
6259  for (unsigned int d = 0; d < spacedim; ++d)
6260  {
6261  p1[d] = data_array[begin_idx + 2 * bbox * spacedim + d];
6262  p2[d] =
6263  data_array[begin_idx + 2 * bbox * spacedim + spacedim + d];
6264  }
6265  BoundingBox<spacedim> loc_bbox(std::make_pair(p1, p2));
6266  global_bboxes[i][bbox] = loc_bbox;
6267  }
6268  // Shifting the first index to the start of the next vector
6269  begin_idx += size_all_data[i];
6270  }
6271  return global_bboxes;
6272 #endif // DEAL_II_WITH_MPI
6273  }
6274 
6275 
6276 
6277  template <int spacedim>
6280  const std::vector<BoundingBox<spacedim>> &local_description,
6281  const MPI_Comm & mpi_communicator)
6282  {
6283 #ifndef DEAL_II_WITH_MPI
6284  (void)mpi_communicator;
6285  // Building a tree with the only boxes available without MPI
6286  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> boxes_index(
6287  local_description.size());
6288  // Adding to each box the rank of the process owning it
6289  for (unsigned int i = 0; i < local_description.size(); ++i)
6290  boxes_index[i] = std::make_pair(local_description[i], 0u);
6291  return pack_rtree(boxes_index);
6292 #else
6293  // Exchanging local bounding boxes
6294  const std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes =
6295  Utilities::MPI::all_gather(mpi_communicator, local_description);
6296 
6297  // Preparing to flatten the vector
6298  const unsigned int n_procs =
6299  Utilities::MPI::n_mpi_processes(mpi_communicator);
6300  // The i'th element of the following vector contains the index of the first
6301  // local bounding box from the process of rank i
6302  std::vector<unsigned int> bboxes_position(n_procs);
6303 
6304  unsigned int tot_bboxes = 0;
6305  for (const auto &process_bboxes : global_bboxes)
6306  tot_bboxes += process_bboxes.size();
6307 
6308  // Now flattening the vector
6309  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
6310  flat_global_bboxes;
6311  flat_global_bboxes.reserve(tot_bboxes);
6312  unsigned int process_index = 0;
6313  for (const auto &process_bboxes : global_bboxes)
6314  {
6315  // Initialize a vector containing bounding boxes and rank of a process
6316  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
6317  boxes_and_indices(process_bboxes.size());
6318 
6319  // Adding to each box the rank of the process owning it
6320  for (unsigned int i = 0; i < process_bboxes.size(); ++i)
6321  boxes_and_indices[i] =
6322  std::make_pair(process_bboxes[i], process_index);
6323 
6324  flat_global_bboxes.insert(flat_global_bboxes.end(),
6325  boxes_and_indices.begin(),
6326  boxes_and_indices.end());
6327 
6328  ++process_index;
6329  }
6330 
6331  // Build a tree out of the bounding boxes. We avoid using the
6332  // insert method so that boost uses the packing algorithm
6333  return RTree<std::pair<BoundingBox<spacedim>, unsigned int>>(
6334  flat_global_bboxes.begin(), flat_global_bboxes.end());
6335 #endif // DEAL_II_WITH_MPI
6336  }
6337 
6338 
6339 
6340  template <int dim, int spacedim>
6341  void
6343  const Triangulation<dim, spacedim> & tria,
6344  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
6345  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group)
6346  {
6347  // 1) determine for each vertex a vertex it concides with and
6348  // put it into a map
6349  {
6350  static const int lookup_table_2d[2][2] =
6351  // flip:
6352  {
6353  {0, 1}, // false
6354  {1, 0} // true
6355  };
6356 
6357  static const int lookup_table_3d[2][2][2][4] =
6358  // orientation flip rotation
6359  {{{
6360  {0, 2, 1, 3}, // false false false
6361  {2, 3, 0, 1} // false false true
6362  },
6363  {
6364  {3, 1, 2, 0}, // false true false
6365  {1, 0, 3, 2} // false true true
6366  }},
6367  {{
6368  {0, 1, 2, 3}, // true false false
6369  {1, 3, 0, 2} // true false true
6370  },
6371  {
6372  {3, 2, 1, 0}, // true true false
6373  {2, 0, 3, 1} // true true true
6374  }}};
6375 
6376  // loop over all periodic face pairs
6377  for (const auto &pair : tria.get_periodic_face_map())
6378  {
6379  if (pair.first.first->level() != pair.second.first.first->level())
6380  continue;
6381 
6382  const auto face_a = pair.first.first->face(pair.first.second);
6383  const auto face_b =
6384  pair.second.first.first->face(pair.second.first.second);
6385  const auto mask = pair.second.second;
6386 
6387  AssertDimension(face_a->n_vertices(), face_b->n_vertices());
6388 
6389  // loop over all vertices on face
6390  for (unsigned int i = 0; i < face_a->n_vertices(); ++i)
6391  {
6392  const bool face_orientation = mask[0];
6393  const bool face_flip = mask[1];
6394  const bool face_rotation = mask[2];
6395 
6396  // find the right local vertex index for the second face
6397  unsigned int j = 0;
6398  switch (dim)
6399  {
6400  case 1:
6401  j = i;
6402  break;
6403  case 2:
6404  j = lookup_table_2d[face_flip][i];
6405  break;
6406  case 3:
6407  j = lookup_table_3d[face_orientation][face_flip]
6408  [face_rotation][i];
6409  break;
6410  default:
6411  AssertThrow(false, ExcNotImplemented());
6412  }
6413 
6414  // get vertex indices and store in map
6415  const auto vertex_a = face_a->vertex_index(i);
6416  const auto vertex_b = face_b->vertex_index(j);
6417  unsigned int temp = std::min(vertex_a, vertex_b);
6418 
6419  auto it_a = vertex_to_coinciding_vertex_group.find(vertex_a);
6420  if (it_a != vertex_to_coinciding_vertex_group.end())
6421  temp = std::min(temp, it_a->second);
6422 
6423  auto it_b = vertex_to_coinciding_vertex_group.find(vertex_b);
6424  if (it_b != vertex_to_coinciding_vertex_group.end())
6425  temp = std::min(temp, it_b->second);
6426 
6427  if (it_a != vertex_to_coinciding_vertex_group.end())
6428  it_a->second = temp;
6429  else
6430  vertex_to_coinciding_vertex_group[vertex_a] = temp;
6431 
6432  if (it_b != vertex_to_coinciding_vertex_group.end())
6433  it_b->second = temp;
6434  else
6435  vertex_to_coinciding_vertex_group[vertex_b] = temp;
6436  }
6437  }
6438 
6439  // 2) compress map: let vertices point to the coinciding vertex with
6440  // the smallest index
6441  for (auto &p : vertex_to_coinciding_vertex_group)
6442  {
6443  if (p.first == p.second)
6444  continue;
6445  unsigned int temp = p.second;
6446  while (temp != vertex_to_coinciding_vertex_group[temp])
6447  temp = vertex_to_coinciding_vertex_group[temp];
6448  p.second = temp;
6449  }
6450 
6451  // 3) create a map: smallest index of coinciding index -> all
6452  // coinciding indices
6453  for (auto p : vertex_to_coinciding_vertex_group)
6454  coinciding_vertex_groups[p.second] = {};
6455 
6456  for (auto p : vertex_to_coinciding_vertex_group)
6457  coinciding_vertex_groups[p.second].push_back(p.first);
6458  }
6459  }
6460 
6461 
6462 
6463  template <int dim, int spacedim>
6464  std::map<unsigned int, std::set<::types::subdomain_id>>
6466  const Triangulation<dim, spacedim> &tria)
6467  {
6468  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
6469  &tria) == nullptr) // nothing to do for a serial triangulation
6470  return {};
6471 
6472  // 1) collect for each vertex on periodic faces all vertices it coincides
6473  // with
6474  std::map<unsigned int, std::vector<unsigned int>> coinciding_vertex_groups;
6475  std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
6476 
6478  coinciding_vertex_groups,
6479  vertex_to_coinciding_vertex_group);
6480 
6481  // 2) collect vertices belonging to local cells
6482  std::vector<bool> vertex_of_own_cell(tria.n_vertices(), false);
6483  for (const auto &cell : tria.active_cell_iterators())
6484  if (cell->is_locally_owned())
6485  for (const unsigned int v : cell->vertex_indices())
6486  vertex_of_own_cell[cell->vertex_index(v)] = true;
6487 
6488  // 3) for each vertex belonging to a locally owned cell all ghost
6489  // neighbors (including the periodic own)
6490  std::map<unsigned int, std::set<types::subdomain_id>> result;
6491 
6492  // loop over all active ghost cells
6493  for (const auto &cell : tria.active_cell_iterators())
6494  if (cell->is_ghost())
6495  {
6496  const types::subdomain_id owner = cell->subdomain_id();
6497 
6498  // loop over all its vertices
6499  for (const unsigned int v : cell->vertex_indices())
6500  {
6501  // set owner if vertex belongs to a local cell
6502  if (vertex_of_own_cell[cell->vertex_index(v)])
6503  result[cel