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