Reference documentation for deal.II version Git 9297d75edf 2020-11-26 18:52:14 +0100
\(\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 - 2020 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>
20 
23 
26 #include <deal.II/dofs/dof_tools.h>
27 
28 #include <deal.II/fe/fe_nothing.h>
29 #include <deal.II/fe/fe_q.h>
30 #include <deal.II/fe/fe_values.h>
31 #include <deal.II/fe/mapping_q.h>
32 #include <deal.II/fe/mapping_q1.h>
34 
39 #include <deal.II/grid/manifold.h>
40 #include <deal.II/grid/tria.h>
43 
47 #include <deal.II/lac/solver_cg.h>
51 #include <deal.II/lac/vector.h>
53 
56 
58 #include <boost/random/mersenne_twister.hpp>
59 #include <boost/random/uniform_real_distribution.hpp>
61 
62 #include <array>
63 #include <cmath>
64 #include <iostream>
65 #include <list>
66 #include <numeric>
67 #include <set>
68 #include <tuple>
69 #include <unordered_map>
70 
72 
73 
74 namespace GridTools
75 {
76  template <int dim, int spacedim>
77  double
79  {
80  // we can't deal with distributed meshes since we don't have all
81  // vertices locally. there is one exception, however: if the mesh has
82  // never been refined. the way to test this is not to ask
83  // tria.n_levels()==1, since this is something that can happen on one
84  // processor without being true on all. however, we can ask for the
85  // global number of active cells and use that
86 #if defined(DEAL_II_WITH_P4EST) && defined(DEBUG)
88  dynamic_cast<
90  Assert(p_tria->n_global_active_cells() == tria.n_cells(0),
92 #endif
93 
94  // the algorithm used simply traverses all cells and picks out the
95  // boundary vertices. it may or may not be faster to simply get all
96  // vectors, don't mark boundary vertices, and compute the distances
97  // thereof, but at least as the mesh is refined, it seems better to
98  // first mark boundary nodes, as marking is O(N) in the number of
99  // cells/vertices, while computing the maximal distance is O(N*N)
100  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
101  std::vector<bool> boundary_vertices(vertices.size(), false);
102 
104  tria.begin_active();
106  tria.end();
107  for (; cell != endc; ++cell)
108  for (const unsigned int face : cell->face_indices())
109  if (cell->face(face)->at_boundary())
110  for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
111  boundary_vertices[cell->face(face)->vertex_index(i)] = true;
112 
113  // now traverse the list of boundary vertices and check distances.
114  // since distances are symmetric, we only have to check one half
115  double max_distance_sqr = 0;
116  std::vector<bool>::const_iterator pi = boundary_vertices.begin();
117  const unsigned int N = boundary_vertices.size();
118  for (unsigned int i = 0; i < N; ++i, ++pi)
119  {
120  std::vector<bool>::const_iterator pj = pi + 1;
121  for (unsigned int j = i + 1; j < N; ++j, ++pj)
122  if ((*pi == true) && (*pj == true) &&
123  ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
124  max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
125  }
126 
127  return std::sqrt(max_distance_sqr);
128  }
129 
130 
131 
132  template <int dim, int spacedim>
133  double
135  const Mapping<dim, spacedim> & mapping)
136  {
137  // get the degree of the mapping if possible. if not, just assume 1
138  unsigned int mapping_degree = 1;
139  if (const auto *p =
140  dynamic_cast<const MappingQGeneric<dim, spacedim> *>(&mapping))
141  mapping_degree = p->get_degree();
142  else if (const auto *p =
143  dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
144  mapping_degree = p->get_degree();
145 
146  // then initialize an appropriate quadrature formula
147  const QGauss<dim> quadrature_formula(mapping_degree + 1);
148  const unsigned int n_q_points = quadrature_formula.size();
149 
150  // we really want the JxW values from the FEValues object, but it
151  // wants a finite element. create a cheap element as a dummy
152  // element
153  FE_Nothing<dim, spacedim> dummy_fe;
154  FEValues<dim, spacedim> fe_values(mapping,
155  dummy_fe,
156  quadrature_formula,
158 
160  cell = triangulation.begin_active(),
161  endc = triangulation.end();
162 
163  double local_volume = 0;
164 
165  // compute the integral quantities by quadrature
166  for (; cell != endc; ++cell)
167  if (cell->is_locally_owned())
168  {
169  fe_values.reinit(cell);
170  for (unsigned int q = 0; q < n_q_points; ++q)
171  local_volume += fe_values.JxW(q);
172  }
173 
174  double global_volume = 0;
175 
176 #ifdef DEAL_II_WITH_MPI
178  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
179  &triangulation))
180  global_volume =
181  Utilities::MPI::sum(local_volume, p_tria->get_communicator());
182  else
183 #endif
184  global_volume = local_volume;
185 
186  return global_volume;
187  }
188 
189 
190 
191  namespace
192  {
207  template <int dim>
208  struct TransformR2UAffine
209  {
210  static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
212  };
213 
214 
215  /*
216  Octave code:
217  M=[0 1; 1 1];
218  K1 = transpose(M) * inverse (M*transpose(M));
219  printf ("{%f, %f},\n", K1' );
220  */
221  template <>
223  [1] = {{-1.000000}, {1.000000}};
224 
225  template <>
227  {1.000000, 0.000000};
228 
229 
230  /*
231  Octave code:
232  M=[0 1 0 1;0 0 1 1;1 1 1 1];
233  K2 = transpose(M) * inverse (M*transpose(M));
234  printf ("{%f, %f, %f},\n", K2' );
235  */
236  template <>
238  [2] = {{-0.500000, -0.500000},
239  {0.500000, -0.500000},
240  {-0.500000, 0.500000},
241  {0.500000, 0.500000}};
242 
243  /*
244  Octave code:
245  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];
246  K3 = transpose(M) * inverse (M*transpose(M))
247  printf ("{%f, %f, %f, %f},\n", K3' );
248  */
249  template <>
251  {0.750000, 0.250000, 0.250000, -0.250000};
252 
253 
254  template <>
256  [3] = {
257  {-0.250000, -0.250000, -0.250000},
258  {0.250000, -0.250000, -0.250000},
259  {-0.250000, 0.250000, -0.250000},
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 
266  };
267 
268 
269  template <>
271  {0.500000,
272  0.250000,
273  0.250000,
274  0.000000,
275  0.250000,
276  0.000000,
277  0.000000,
278  -0.250000};
279  } // namespace
280 
281 
282 
283  template <int dim, int spacedim>
284  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
286  {
288 
289  // A = vertex * KA
291 
292  for (unsigned int d = 0; d < spacedim; ++d)
293  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
294  for (unsigned int e = 0; e < dim; ++e)
295  A[d][e] += vertices[v][d] * TransformR2UAffine<dim>::KA[v][e];
296 
297  // b = vertex * Kb
299  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
301 
302  return std::make_pair(A, b);
303  }
304 
305 
306 
307  template <int dim>
308  Vector<double>
311  const Quadrature<dim> & quadrature)
312  {
313  FE_Nothing<dim> fe;
314  FEValues<dim> fe_values(mapping, fe, quadrature, update_jacobians);
315 
316  Vector<double> aspect_ratio_vector(triangulation.n_active_cells());
317 
318  // loop over cells of processor
319  for (const auto &cell : triangulation.active_cell_iterators())
320  {
321  if (cell->is_locally_owned())
322  {
323  double aspect_ratio_cell = 0.0;
324 
325  fe_values.reinit(cell);
326 
327  // loop over quadrature points
328  for (unsigned int q = 0; q < quadrature.size(); ++q)
329  {
330  const Tensor<2, dim, double> jacobian =
331  Tensor<2, dim, double>(fe_values.jacobian(q));
332 
333  // We intentionally do not want to throw an exception in case of
334  // inverted elements since this is not the task of this
335  // function. Instead, inf is written into the vector in case of
336  // inverted elements.
337  if (determinant(jacobian) <= 0)
338  {
339  aspect_ratio_cell = std::numeric_limits<double>::infinity();
340  }
341  else
342  {
344  for (unsigned int i = 0; i < dim; i++)
345  for (unsigned int j = 0; j < dim; j++)
346  J(i, j) = jacobian[i][j];
347 
348  J.compute_svd();
349 
350  double const max_sv = J.singular_value(0);
351  double const min_sv = J.singular_value(dim - 1);
352  double const ar = max_sv / min_sv;
353 
354  // Take the max between the previous and the current
355  // aspect ratio value; if we had previously encountered
356  // an inverted cell, we will have placed an infinity
357  // in the aspect_ratio_cell variable, and that value
358  // will survive this max operation.
359  aspect_ratio_cell = std::max(aspect_ratio_cell, ar);
360  }
361  }
362 
363  // fill vector
364  aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
365  }
366  }
367 
368  return aspect_ratio_vector;
369  }
370 
371 
372 
373  template <int dim>
374  double
377  const Quadrature<dim> & quadrature)
378  {
379  Vector<double> aspect_ratio_vector =
380  compute_aspect_ratio_of_cells(mapping, triangulation, quadrature);
381 
382  return VectorTools::compute_global_error(triangulation,
383  aspect_ratio_vector,
385  }
386 
387 
388 
389  template <int dim, int spacedim>
392  {
393  using iterator =
395  const auto predicate = [](const iterator &) { return true; };
396 
397  return compute_bounding_box(
398  tria, std::function<bool(const iterator &)>(predicate));
399  }
400 
401 
402 
403  // Generic functions for appending face data in 2D or 3D. TODO: we can
404  // remove these once we have 'if constexpr'.
405  namespace internal
406  {
407  inline void
408  append_face_data(const CellData<1> &face_data, SubCellData &subcell_data)
409  {
410  subcell_data.boundary_lines.push_back(face_data);
411  }
412 
413 
414 
415  inline void
416  append_face_data(const CellData<2> &face_data, SubCellData &subcell_data)
417  {
418  subcell_data.boundary_quads.push_back(face_data);
419  }
420 
421 
422 
423  // Lexical comparison for sorting CellData objects.
424  template <int structdim>
426  {
427  bool
429  const CellData<structdim> &b) const
430  {
431  // Check vertices:
432  if (std::lexicographical_compare(std::begin(a.vertices),
433  std::end(a.vertices),
434  std::begin(b.vertices),
435  std::end(b.vertices)))
436  return true;
437  // it should never be necessary to check the material or manifold
438  // ids as a 'tiebreaker' (since they must be equal if the vertex
439  // indices are equal). Assert it anyway:
440 #ifdef DEBUG
441  if (std::equal(std::begin(a.vertices),
442  std::end(a.vertices),
443  std::begin(b.vertices)))
444  {
445  Assert(a.material_id == b.material_id &&
446  a.manifold_id == b.manifold_id,
447  ExcMessage(
448  "Two CellData objects with equal vertices must "
449  "have the same material/boundary ids and manifold "
450  "ids."));
451  }
452 #endif
453  return false;
454  }
455  };
456 
457 
467  template <int dim>
469  {
470  public:
474  template <class FaceIteratorType>
475  void
476  insert_face_data(const FaceIteratorType &face)
477  {
478  CellData<dim - 1> face_cell_data;
479  for (unsigned int vertex_n = 0; vertex_n < face->n_vertices();
480  ++vertex_n)
481  face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
482  face_cell_data.boundary_id = face->boundary_id();
483  face_cell_data.manifold_id = face->manifold_id();
484 
485  face_data.insert(face_cell_data);
486  }
487 
492  get()
493  {
494  SubCellData subcell_data;
495 
496  for (const CellData<dim - 1> &face_cell_data : face_data)
497  internal::append_face_data(face_cell_data, subcell_data);
498  return subcell_data;
499  }
500 
501 
502  private:
505  };
506 
507 
508  // Do nothing for dim=1:
509  template <>
510  class FaceDataHelper<1>
511  {
512  public:
513  template <class FaceIteratorType>
514  void
515  insert_face_data(const FaceIteratorType &)
516  {}
517 
519  get()
520  {
521  return SubCellData();
522  }
523  };
524  } // namespace internal
525 
526 
527 
528  template <int dim, int spacedim>
529  std::
530  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
532  {
533  Assert(1 <= tria.n_levels(),
534  ExcMessage("The input triangulation must be non-empty."));
535 
536  std::vector<Point<spacedim>> vertices;
537  std::vector<CellData<dim>> cells;
538 
539  unsigned int max_level_0_vertex_n = 0;
540  for (const auto &cell : tria.cell_iterators_on_level(0))
541  for (const unsigned int cell_vertex_n : cell->vertex_indices())
542  max_level_0_vertex_n =
543  std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
544  vertices.resize(max_level_0_vertex_n + 1);
545 
547  std::set<CellData<1>, internal::CellDataComparator<1>>
548  line_data; // only used in 3D
549 
550  for (const auto &cell : tria.cell_iterators_on_level(0))
551  {
552  // Save cell data
553  CellData<dim> cell_data;
554  for (const unsigned int cell_vertex_n : cell->vertex_indices())
555  {
556  Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
557  ExcInternalError());
558  vertices[cell->vertex_index(cell_vertex_n)] =
559  cell->vertex(cell_vertex_n);
560  cell_data.vertices[cell_vertex_n] =
561  cell->vertex_index(cell_vertex_n);
562  }
563  cell_data.material_id = cell->material_id();
564  cell_data.manifold_id = cell->manifold_id();
565  cells.push_back(cell_data);
566 
567  // Save face data
568  if (dim > 1)
569  {
570  for (const unsigned int face_n : cell->face_indices())
571  face_data.insert_face_data(cell->face(face_n));
572  }
573  // Save line data
574  if (dim == 3)
575  {
576  for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
577  {
578  const auto line = cell->line(line_n);
579  CellData<1> line_cell_data;
580  for (unsigned int vertex_n = 0; vertex_n < line->n_vertices();
581  ++vertex_n)
582  line_cell_data.vertices[vertex_n] =
583  line->vertex_index(vertex_n);
584  line_cell_data.boundary_id = line->boundary_id();
585  line_cell_data.manifold_id = line->manifold_id();
586 
587  line_data.insert(line_cell_data);
588  }
589  }
590  }
591 
592  // Double-check that there are no unused vertices:
593 #ifdef DEBUG
594  {
595  std::vector<bool> used_vertices(vertices.size());
596  for (const CellData<dim> &cell_data : cells)
597  for (const auto v : cell_data.vertices)
598  used_vertices[v] = true;
599  Assert(std::find(used_vertices.begin(), used_vertices.end(), false) ==
600  used_vertices.end(),
601  ExcMessage("The level zero vertices should form a contiguous "
602  "range."));
603  }
604 #endif
605 
606  SubCellData subcell_data = face_data.get();
607 
608  if (dim == 3)
609  for (const CellData<1> &face_line_data : line_data)
610  subcell_data.boundary_lines.push_back(face_line_data);
611 
612  return std::tuple<std::vector<Point<spacedim>>,
613  std::vector<CellData<dim>>,
614  SubCellData>(std::move(vertices),
615  std::move(cells),
616  std::move(subcell_data));
617  }
618 
619 
620 
621  template <int dim, int spacedim>
622  void
624  std::vector<CellData<dim>> & cells,
625  SubCellData & subcelldata)
626  {
627  Assert(
628  subcelldata.check_consistency(dim),
629  ExcMessage(
630  "Invalid SubCellData supplied according to ::check_consistency(). "
631  "This is caused by data containing objects for the wrong dimension."));
632 
633  // first check which vertices are actually used
634  std::vector<bool> vertex_used(vertices.size(), false);
635  for (unsigned int c = 0; c < cells.size(); ++c)
636  for (unsigned int v = 0; v < cells[c].vertices.size(); ++v)
637  {
638  Assert(cells[c].vertices[v] < vertices.size(),
639  ExcMessage("Invalid vertex index encountered! cells[" +
640  Utilities::int_to_string(c) + "].vertices[" +
641  Utilities::int_to_string(v) + "]=" +
642  Utilities::int_to_string(cells[c].vertices[v]) +
643  " is invalid, because only " +
645  " vertices were supplied."));
646  vertex_used[cells[c].vertices[v]] = true;
647  }
648 
649 
650  // then renumber the vertices that are actually used in the same order as
651  // they were beforehand
652  const unsigned int invalid_vertex = numbers::invalid_unsigned_int;
653  std::vector<unsigned int> new_vertex_numbers(vertices.size(),
654  invalid_vertex);
655  unsigned int next_free_number = 0;
656  for (unsigned int i = 0; i < vertices.size(); ++i)
657  if (vertex_used[i] == true)
658  {
659  new_vertex_numbers[i] = next_free_number;
660  ++next_free_number;
661  }
662 
663  // next replace old vertex numbers by the new ones
664  for (unsigned int c = 0; c < cells.size(); ++c)
665  for (auto &v : cells[c].vertices)
666  v = new_vertex_numbers[v];
667 
668  // same for boundary data
669  for (unsigned int c = 0; c < subcelldata.boundary_lines.size(); // NOLINT
670  ++c)
671  for (unsigned int v = 0;
672  v < subcelldata.boundary_lines[c].vertices.size();
673  ++v)
674  {
675  Assert(subcelldata.boundary_lines[c].vertices[v] <
676  new_vertex_numbers.size(),
677  ExcMessage(
678  "Invalid vertex index in subcelldata.boundary_lines. "
679  "subcelldata.boundary_lines[" +
680  Utilities::int_to_string(c) + "].vertices[" +
681  Utilities::int_to_string(v) + "]=" +
683  subcelldata.boundary_lines[c].vertices[v]) +
684  " is invalid, because only " +
685  Utilities::int_to_string(vertices.size()) +
686  " vertices were supplied."));
687  subcelldata.boundary_lines[c].vertices[v] =
688  new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]];
689  }
690 
691  for (unsigned int c = 0; c < subcelldata.boundary_quads.size(); // NOLINT
692  ++c)
693  for (unsigned int v = 0;
694  v < subcelldata.boundary_quads[c].vertices.size();
695  ++v)
696  {
697  Assert(subcelldata.boundary_quads[c].vertices[v] <
698  new_vertex_numbers.size(),
699  ExcMessage(
700  "Invalid vertex index in subcelldata.boundary_quads. "
701  "subcelldata.boundary_quads[" +
702  Utilities::int_to_string(c) + "].vertices[" +
703  Utilities::int_to_string(v) + "]=" +
705  subcelldata.boundary_quads[c].vertices[v]) +
706  " is invalid, because only " +
707  Utilities::int_to_string(vertices.size()) +
708  " vertices were supplied."));
709 
710  subcelldata.boundary_quads[c].vertices[v] =
711  new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]];
712  }
713 
714  // finally copy over the vertices which we really need to a new array and
715  // replace the old one by the new one
716  std::vector<Point<spacedim>> tmp;
717  tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(), true));
718  for (unsigned int v = 0; v < vertices.size(); ++v)
719  if (vertex_used[v] == true)
720  tmp.push_back(vertices[v]);
721  swap(vertices, tmp);
722  }
723 
724 
725 
726  template <int dim, int spacedim>
727  void
729  std::vector<CellData<dim>> & cells,
730  SubCellData & subcelldata,
731  std::vector<unsigned int> & considered_vertices,
732  const double tol)
733  {
734  AssertIndexRange(2, vertices.size());
735  // create a vector of vertex indices. initialize it to the identity, later
736  // on change that if necessary.
737  std::vector<unsigned int> new_vertex_numbers(vertices.size());
738  std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
739 
740  // if the considered_vertices vector is empty, consider all vertices
741  if (considered_vertices.size() == 0)
742  considered_vertices = new_vertex_numbers;
743  Assert(considered_vertices.size() <= vertices.size(), ExcInternalError());
744 
745  // The algorithm below improves upon the naive O(n^2) algorithm by first
746  // sorting vertices by their value in one component and then only
747  // comparing vertices for equality which are nearly equal in that
748  // component. For example, if @p vertices form a cube, then we will only
749  // compare points that have the same x coordinate when we try to find
750  // duplicated vertices.
751 
752  // Start by finding the longest coordinate direction. This minimizes the
753  // number of points that need to be compared against each-other in a
754  // single set for typical geometries.
755  const BoundingBox<spacedim> bbox(vertices);
756  const auto & min = bbox.get_boundary_points().first;
757  const auto & max = bbox.get_boundary_points().second;
758 
759  unsigned int longest_coordinate_direction = 0;
760  double longest_coordinate_length = max[0] - min[0];
761  for (unsigned int d = 1; d < spacedim; ++d)
762  {
763  const double coordinate_length = max[d] - min[d];
764  if (longest_coordinate_length < coordinate_length)
765  {
766  longest_coordinate_length = coordinate_length;
767  longest_coordinate_direction = d;
768  }
769  }
770 
771  // Sort vertices (while preserving their vertex numbers) along that
772  // coordinate direction:
773  std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
774  sorted_vertices.reserve(vertices.size());
775  for (const unsigned int vertex_n : considered_vertices)
776  {
777  AssertIndexRange(vertex_n, vertices.size());
778  sorted_vertices.emplace_back(vertex_n, vertices[vertex_n]);
779  }
780  std::sort(sorted_vertices.begin(),
781  sorted_vertices.end(),
782  [&](const std::pair<unsigned int, Point<spacedim>> &a,
783  const std::pair<unsigned int, Point<spacedim>> &b) {
784  return a.second[longest_coordinate_direction] <
785  b.second[longest_coordinate_direction];
786  });
787 
788  auto within_tolerance = [=](const Point<spacedim> &a,
789  const Point<spacedim> &b) {
790  for (unsigned int d = 0; d < spacedim; ++d)
791  if (std::abs(a[d] - b[d]) > tol)
792  return false;
793  return true;
794  };
795 
796  // Find a range of numbers that have the same component in the longest
797  // coordinate direction:
798  auto range_start = sorted_vertices.begin();
799  while (range_start != sorted_vertices.end())
800  {
801  auto range_end = range_start + 1;
802  while (range_end != sorted_vertices.end() &&
803  std::abs(range_end->second[longest_coordinate_direction] -
804  range_start->second[longest_coordinate_direction]) <
805  tol)
806  ++range_end;
807 
808  // preserve behavior with older versions of this function by replacing
809  // higher vertex numbers by lower vertex numbers
810  std::sort(range_start,
811  range_end,
812  [](const std::pair<unsigned int, Point<spacedim>> &a,
813  const std::pair<unsigned int, Point<spacedim>> &b) {
814  return a.first < b.first;
815  });
816 
817  // Now de-duplicate [range_start, range_end)
818  //
819  // We have identified all points that are within a strip of width 'tol'
820  // in one coordinate direction. Now we need to figure out which of these
821  // are also close in other coordinate directions. If two are close, we
822  // can mark the second one for deletion.
823  for (auto reference = range_start; reference != range_end; ++reference)
824  {
825  if (reference->first != numbers::invalid_unsigned_int)
826  for (auto it = reference + 1; it != range_end; ++it)
827  {
828  if (within_tolerance(reference->second, it->second))
829  {
830  new_vertex_numbers[it->first] = reference->first;
831  // skip the replaced vertex in the future
832  it->first = numbers::invalid_unsigned_int;
833  }
834  }
835  }
836  range_start = range_end;
837  }
838 
839  // now we got a renumbering list. simply renumber all vertices
840  // (non-duplicate vertices get renumbered to themselves, so nothing bad
841  // happens). after that, the duplicate vertices will be unused, so call
842  // delete_unused_vertices() to do that part of the job.
843  for (auto &cell : cells)
844  for (auto &vertex_index : cell.vertices)
845  vertex_index = new_vertex_numbers[vertex_index];
846  for (auto &quad : subcelldata.boundary_quads)
847  for (auto &vertex_index : quad.vertices)
848  vertex_index = new_vertex_numbers[vertex_index];
849  for (auto &line : subcelldata.boundary_lines)
850  for (auto &vertex_index : line.vertices)
851  vertex_index = new_vertex_numbers[vertex_index];
852 
853  delete_unused_vertices(vertices, cells, subcelldata);
854  }
855 
856 
857 
858  // define some transformations
859  namespace internal
860  {
861  template <int spacedim>
862  class Shift
863  {
864  public:
865  explicit Shift(const Tensor<1, spacedim> &shift)
866  : shift(shift)
867  {}
870  {
871  return p + shift;
872  }
873 
874  private:
876  };
877 
878 
879  // Transformation to rotate around one of the cartesian axes.
880  class Rotate3d
881  {
882  public:
883  Rotate3d(const double angle, const unsigned int axis)
884  : angle(angle)
885  , axis(axis)
886  {}
887 
888  Point<3>
889  operator()(const Point<3> &p) const
890  {
891  if (axis == 0)
892  return {p(0),
893  std::cos(angle) * p(1) - std::sin(angle) * p(2),
894  std::sin(angle) * p(1) + std::cos(angle) * p(2)};
895  else if (axis == 1)
896  return {std::cos(angle) * p(0) + std::sin(angle) * p(2),
897  p(1),
898  -std::sin(angle) * p(0) + std::cos(angle) * p(2)};
899  else
900  return {std::cos(angle) * p(0) - std::sin(angle) * p(1),
901  std::sin(angle) * p(0) + std::cos(angle) * p(1),
902  p(2)};
903  }
904 
905  private:
906  const double angle;
907  const unsigned int axis;
908  };
909 
910  template <int spacedim>
911  class Scale
912  {
913  public:
914  explicit Scale(const double factor)
915  : factor(factor)
916  {}
919  {
920  return p * factor;
921  }
922 
923  private:
924  const double factor;
925  };
926  } // namespace internal
927 
928 
929  template <int dim, int spacedim>
930  void
931  shift(const Tensor<1, spacedim> & shift_vector,
933  {
934  transform(internal::Shift<spacedim>(shift_vector), triangulation);
935  }
936 
937 
938  template <int dim>
939  void
940  rotate(const double angle,
941  const unsigned int axis,
943  {
944  Assert(axis < 3, ExcMessage("Invalid axis given!"));
945 
946  transform(internal::Rotate3d(angle, axis), triangulation);
947  }
948 
949  template <int dim, int spacedim>
950  void
951  scale(const double scaling_factor,
953  {
954  Assert(scaling_factor > 0, ExcScalingFactorNotPositive(scaling_factor));
955  transform(internal::Scale<spacedim>(scaling_factor), triangulation);
956  }
957 
958 
959  namespace internal
960  {
966  inline void
968  const AffineConstraints<double> &constraints,
969  Vector<double> & u)
970  {
971  const unsigned int n_dofs = S.n();
972  const auto op = linear_operator(S);
973  const auto SF = constrained_linear_operator(constraints, op);
975  prec.initialize(S, 1.2);
976 
977  SolverControl control(n_dofs, 1.e-10, false, false);
979  SolverCG<Vector<double>> solver(control, mem);
980 
981  Vector<double> f(n_dofs);
982 
983  const auto constrained_rhs =
984  constrained_right_hand_side(constraints, op, f);
985  solver.solve(SF, u, constrained_rhs, prec);
986 
987  constraints.distribute(u);
988  }
989  } // namespace internal
990 
991 
992  // Implementation for dimensions except 1
993  template <int dim>
994  void
995  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
997  const Function<dim> * coefficient,
998  const bool solve_for_absolute_positions)
999  {
1000  if (dim == 1)
1001  Assert(false, ExcNotImplemented());
1002 
1003  // first provide everything that is needed for solving a Laplace
1004  // equation.
1005  FE_Q<dim> q1(1);
1006 
1007  DoFHandler<dim> dof_handler(triangulation);
1008  dof_handler.distribute_dofs(q1);
1009 
1010  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
1011  DoFTools::make_sparsity_pattern(dof_handler, dsp);
1012  dsp.compress();
1013 
1014  SparsityPattern sparsity_pattern;
1015  sparsity_pattern.copy_from(dsp);
1016  sparsity_pattern.compress();
1017 
1018  SparseMatrix<double> S(sparsity_pattern);
1019 
1020  QGauss<dim> quadrature(4);
1021 
1023  StaticMappingQ1<dim>::mapping, dof_handler, quadrature, S, coefficient);
1024 
1025  // set up the boundary values for the laplace problem
1026  std::array<AffineConstraints<double>, dim> constraints;
1027  typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
1028  new_points.end();
1029 
1030  // fill these maps using the data given by new_points
1031  for (const auto &cell : dof_handler.active_cell_iterators())
1032  {
1033  // loop over all vertices of the cell and see if it is listed in the map
1034  // given as first argument of the function
1035  for (const unsigned int vertex_no : cell->vertex_indices())
1036  {
1037  const unsigned int vertex_index = cell->vertex_index(vertex_no);
1038  const Point<dim> & vertex_point = cell->vertex(vertex_no);
1039 
1040  const typename std::map<unsigned int, Point<dim>>::const_iterator
1041  map_iter = new_points.find(vertex_index);
1042 
1043  if (map_iter != map_end)
1044  for (unsigned int i = 0; i < dim; ++i)
1045  {
1046  constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
1047  constraints[i].set_inhomogeneity(
1048  cell->vertex_dof_index(vertex_no, 0),
1049  (solve_for_absolute_positions ?
1050  map_iter->second(i) :
1051  map_iter->second(i) - vertex_point[i]));
1052  }
1053  }
1054  }
1055 
1056  for (unsigned int i = 0; i < dim; ++i)
1057  constraints[i].close();
1058 
1059  // solve the dim problems with different right hand sides.
1060  Vector<double> us[dim];
1061  for (unsigned int i = 0; i < dim; ++i)
1062  us[i].reinit(dof_handler.n_dofs());
1063 
1064  // solve linear systems in parallel
1065  Threads::TaskGroup<> tasks;
1066  for (unsigned int i = 0; i < dim; ++i)
1067  tasks +=
1068  Threads::new_task(&internal::laplace_solve, S, constraints[i], us[i]);
1069  tasks.join_all();
1070 
1071  // change the coordinates of the points of the triangulation
1072  // according to the computed values
1073  std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
1074  for (const auto &cell : dof_handler.active_cell_iterators())
1075  for (const unsigned int vertex_no : cell->vertex_indices())
1076  if (vertex_touched[cell->vertex_index(vertex_no)] == false)
1077  {
1078  Point<dim> &v = cell->vertex(vertex_no);
1079 
1080  const types::global_dof_index dof_index =
1081  cell->vertex_dof_index(vertex_no, 0);
1082  for (unsigned int i = 0; i < dim; ++i)
1083  if (solve_for_absolute_positions)
1084  v(i) = us[i](dof_index);
1085  else
1086  v(i) += us[i](dof_index);
1087 
1088  vertex_touched[cell->vertex_index(vertex_no)] = true;
1089  }
1090  }
1091 
1092  template <int dim, int spacedim>
1093  std::map<unsigned int, Point<spacedim>>
1095  {
1096  std::map<unsigned int, Point<spacedim>> vertex_map;
1098  cell = tria.begin_active(),
1099  endc = tria.end();
1100  for (; cell != endc; ++cell)
1101  {
1102  for (unsigned int i : cell->face_indices())
1103  {
1104  const typename Triangulation<dim, spacedim>::face_iterator &face =
1105  cell->face(i);
1106  if (face->at_boundary())
1107  {
1108  for (unsigned j = 0; j < face->n_vertices(); ++j)
1109  {
1110  const Point<spacedim> &vertex = face->vertex(j);
1111  const unsigned int vertex_index = face->vertex_index(j);
1112  vertex_map[vertex_index] = vertex;
1113  }
1114  }
1115  }
1116  }
1117  return vertex_map;
1118  }
1119 
1124  template <int dim, int spacedim>
1125  void
1126  distort_random(const double factor,
1128  const bool keep_boundary,
1129  const unsigned int seed)
1130  {
1131  // if spacedim>dim we need to make sure that we perturb
1132  // points but keep them on
1133  // the manifold. however, this isn't implemented right now
1134  Assert(spacedim == dim, ExcNotImplemented());
1135 
1136 
1137  // find the smallest length of the
1138  // lines adjacent to the
1139  // vertex. take the initial value
1140  // to be larger than anything that
1141  // might be found: the diameter of
1142  // the triangulation, here
1143  // estimated by adding up the
1144  // diameters of the coarse grid
1145  // cells.
1146  double almost_infinite_length = 0;
1147  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
1148  triangulation.begin(0);
1149  cell != triangulation.end(0);
1150  ++cell)
1151  almost_infinite_length += cell->diameter();
1152 
1153  std::vector<double> minimal_length(triangulation.n_vertices(),
1154  almost_infinite_length);
1155 
1156  // also note if a vertex is at the boundary
1157  std::vector<bool> at_boundary(keep_boundary ? triangulation.n_vertices() :
1158  0,
1159  false);
1160  // for parallel::shared::Triangulation we need to work on all vertices,
1161  // not just the ones related to locally owned cells;
1162  const bool is_parallel_shared =
1164  &triangulation) != nullptr);
1165  for (const auto &cell : triangulation.active_cell_iterators())
1166  if (is_parallel_shared || cell->is_locally_owned())
1167  {
1168  if (dim > 1)
1169  {
1170  for (unsigned int i = 0; i < cell->n_lines(); ++i)
1171  {
1173  line = cell->line(i);
1174 
1175  if (keep_boundary && line->at_boundary())
1176  {
1177  at_boundary[line->vertex_index(0)] = true;
1178  at_boundary[line->vertex_index(1)] = true;
1179  }
1180 
1181  minimal_length[line->vertex_index(0)] =
1182  std::min(line->diameter(),
1183  minimal_length[line->vertex_index(0)]);
1184  minimal_length[line->vertex_index(1)] =
1185  std::min(line->diameter(),
1186  minimal_length[line->vertex_index(1)]);
1187  }
1188  }
1189  else // dim==1
1190  {
1191  if (keep_boundary)
1192  for (unsigned int vertex = 0; vertex < 2; ++vertex)
1193  if (cell->at_boundary(vertex) == true)
1194  at_boundary[cell->vertex_index(vertex)] = true;
1195 
1196  minimal_length[cell->vertex_index(0)] =
1197  std::min(cell->diameter(),
1198  minimal_length[cell->vertex_index(0)]);
1199  minimal_length[cell->vertex_index(1)] =
1200  std::min(cell->diameter(),
1201  minimal_length[cell->vertex_index(1)]);
1202  }
1203  }
1204 
1205  // create a random number generator for the interval [-1,1]
1206  boost::random::mt19937 rng(seed);
1207  boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
1208 
1209  // If the triangulation is distributed, we need to
1210  // exchange the moved vertices across mpi processes
1212  *distributed_triangulation =
1214  &triangulation))
1215  {
1216  const std::vector<bool> locally_owned_vertices =
1217  get_locally_owned_vertices(triangulation);
1218  std::vector<bool> vertex_moved(triangulation.n_vertices(), false);
1219 
1220  // Next move vertices on locally owned cells
1221  for (const auto &cell : triangulation.active_cell_iterators())
1222  if (cell->is_locally_owned())
1223  {
1224  for (const unsigned int vertex_no : cell->vertex_indices())
1225  {
1226  const unsigned global_vertex_no =
1227  cell->vertex_index(vertex_no);
1228 
1229  // ignore this vertex if we shall keep the boundary and
1230  // this vertex *is* at the boundary, if it is already moved
1231  // or if another process moves this vertex
1232  if ((keep_boundary && at_boundary[global_vertex_no]) ||
1233  vertex_moved[global_vertex_no] ||
1234  !locally_owned_vertices[global_vertex_no])
1235  continue;
1236 
1237  // first compute a random shift vector
1238  Point<spacedim> shift_vector;
1239  for (unsigned int d = 0; d < spacedim; ++d)
1240  shift_vector(d) = uniform_distribution(rng);
1241 
1242  shift_vector *= factor * minimal_length[global_vertex_no] /
1243  std::sqrt(shift_vector.square());
1244 
1245  // finally move the vertex
1246  cell->vertex(vertex_no) += shift_vector;
1247  vertex_moved[global_vertex_no] = true;
1248  }
1249  }
1250 
1251 #ifdef DEAL_II_WITH_P4EST
1252  distributed_triangulation->communicate_locally_moved_vertices(
1253  locally_owned_vertices);
1254 #else
1255  (void)distributed_triangulation;
1256  Assert(false, ExcInternalError());
1257 #endif
1258  }
1259  else
1260  // if this is a sequential triangulation, we could in principle
1261  // use the algorithm above, but we'll use an algorithm that we used
1262  // before the parallel::distributed::Triangulation was introduced
1263  // in order to preserve backward compatibility
1264  {
1265  // loop over all vertices and compute their new locations
1266  const unsigned int n_vertices = triangulation.n_vertices();
1267  std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
1268  const std::vector<Point<spacedim>> &old_vertex_locations =
1269  triangulation.get_vertices();
1270 
1271  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1272  {
1273  // ignore this vertex if we will keep the boundary and
1274  // this vertex *is* at the boundary
1275  if (keep_boundary && at_boundary[vertex])
1276  new_vertex_locations[vertex] = old_vertex_locations[vertex];
1277  else
1278  {
1279  // compute a random shift vector
1280  Point<spacedim> shift_vector;
1281  for (unsigned int d = 0; d < spacedim; ++d)
1282  shift_vector(d) = uniform_distribution(rng);
1283 
1284  shift_vector *= factor * minimal_length[vertex] /
1285  std::sqrt(shift_vector.square());
1286 
1287  // record new vertex location
1288  new_vertex_locations[vertex] =
1289  old_vertex_locations[vertex] + shift_vector;
1290  }
1291  }
1292 
1293  // now do the actual move of the vertices
1294  for (const auto &cell : triangulation.active_cell_iterators())
1295  for (const unsigned int vertex_no : cell->vertex_indices())
1296  cell->vertex(vertex_no) =
1297  new_vertex_locations[cell->vertex_index(vertex_no)];
1298  }
1299 
1300  // Correct hanging nodes if necessary
1301  if (dim >= 2)
1302  {
1303  // We do the same as in GridTools::transform
1304  //
1305  // exclude hanging nodes at the boundaries of artificial cells:
1306  // these may belong to ghost cells for which we know the exact
1307  // location of vertices, whereas the artificial cell may or may
1308  // not be further refined, and so we cannot know whether
1309  // the location of the hanging node is correct or not
1311  cell = triangulation.begin_active(),
1312  endc = triangulation.end();
1313  for (; cell != endc; ++cell)
1314  if (!cell->is_artificial())
1315  for (const unsigned int face : cell->face_indices())
1316  if (cell->face(face)->has_children() &&
1317  !cell->face(face)->at_boundary())
1318  {
1319  // this face has hanging nodes
1320  if (dim == 2)
1321  cell->face(face)->child(0)->vertex(1) =
1322  (cell->face(face)->vertex(0) +
1323  cell->face(face)->vertex(1)) /
1324  2;
1325  else if (dim == 3)
1326  {
1327  cell->face(face)->child(0)->vertex(1) =
1328  .5 * (cell->face(face)->vertex(0) +
1329  cell->face(face)->vertex(1));
1330  cell->face(face)->child(0)->vertex(2) =
1331  .5 * (cell->face(face)->vertex(0) +
1332  cell->face(face)->vertex(2));
1333  cell->face(face)->child(1)->vertex(3) =
1334  .5 * (cell->face(face)->vertex(1) +
1335  cell->face(face)->vertex(3));
1336  cell->face(face)->child(2)->vertex(3) =
1337  .5 * (cell->face(face)->vertex(2) +
1338  cell->face(face)->vertex(3));
1339 
1340  // center of the face
1341  cell->face(face)->child(0)->vertex(3) =
1342  .25 * (cell->face(face)->vertex(0) +
1343  cell->face(face)->vertex(1) +
1344  cell->face(face)->vertex(2) +
1345  cell->face(face)->vertex(3));
1346  }
1347  }
1348  }
1349  }
1350 
1351 
1352 
1353  template <int dim, template <int, int> class MeshType, int spacedim>
1354  unsigned int
1355  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1356  const Point<spacedim> & p,
1357  const std::vector<bool> & marked_vertices)
1358  {
1359  // first get the underlying triangulation from the mesh and determine
1360  // vertices and used vertices
1361  const Triangulation<dim, spacedim> &tria = mesh.get_triangulation();
1362 
1363  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
1364 
1365  Assert(tria.get_vertices().size() == marked_vertices.size() ||
1366  marked_vertices.size() == 0,
1367  ExcDimensionMismatch(tria.get_vertices().size(),
1368  marked_vertices.size()));
1369 
1370  // marked_vertices is expected to be a subset of used_vertices. Thus,
1371  // comparing the range marked_vertices.begin() to marked_vertices.end() with
1372  // the range used_vertices.begin() to used_vertices.end() the element in the
1373  // second range must be valid if the element in the first range is valid.
1374  Assert(
1375  marked_vertices.size() == 0 ||
1376  std::equal(marked_vertices.begin(),
1377  marked_vertices.end(),
1378  tria.get_used_vertices().begin(),
1379  [](bool p, bool q) { return !p || q; }),
1380  ExcMessage(
1381  "marked_vertices should be a subset of used vertices in the triangulation "
1382  "but marked_vertices contains one or more vertices that are not used vertices!"));
1383 
1384  // If marked_indices is empty, consider all used_vertices for finding the
1385  // closest vertex to the point. Otherwise, marked_indices is used.
1386  const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
1387  tria.get_used_vertices() :
1388  marked_vertices;
1389 
1390  // At the beginning, the first used vertex is considered to be the closest
1391  // one.
1392  std::vector<bool>::const_iterator first =
1393  std::find(vertices_to_use.begin(), vertices_to_use.end(), true);
1394 
1395  // Assert that at least one vertex is actually used
1396  Assert(first != vertices_to_use.end(), ExcInternalError());
1397 
1398  unsigned int best_vertex = std::distance(vertices_to_use.begin(), first);
1399  double best_dist = (p - vertices[best_vertex]).norm_square();
1400 
1401  // For all remaining vertices, test
1402  // whether they are any closer
1403  for (unsigned int j = best_vertex + 1; j < vertices.size(); j++)
1404  if (vertices_to_use[j])
1405  {
1406  const double dist = (p - vertices[j]).norm_square();
1407  if (dist < best_dist)
1408  {
1409  best_vertex = j;
1410  best_dist = dist;
1411  }
1412  }
1413 
1414  return best_vertex;
1415  }
1416 
1417 
1418 
1419  template <int dim, template <int, int> class MeshType, int spacedim>
1420  unsigned int
1422  const MeshType<dim, spacedim> &mesh,
1423  const Point<spacedim> & p,
1424  const std::vector<bool> & marked_vertices)
1425  {
1426  // Take a shortcut in the simple case.
1427  if (mapping.preserves_vertex_locations() == true)
1428  return find_closest_vertex(mesh, p, marked_vertices);
1429 
1430  // first get the underlying triangulation from the mesh and determine
1431  // vertices and used vertices
1432  const Triangulation<dim, spacedim> &tria = mesh.get_triangulation();
1433 
1434  auto vertices = extract_used_vertices(tria, mapping);
1435 
1436  Assert(tria.get_vertices().size() == marked_vertices.size() ||
1437  marked_vertices.size() == 0,
1438  ExcDimensionMismatch(tria.get_vertices().size(),
1439  marked_vertices.size()));
1440 
1441  // marked_vertices is expected to be a subset of used_vertices. Thus,
1442  // comparing the range marked_vertices.begin() to marked_vertices.end()
1443  // with the range used_vertices.begin() to used_vertices.end() the element
1444  // in the second range must be valid if the element in the first range is
1445  // valid.
1446  Assert(
1447  marked_vertices.size() == 0 ||
1448  std::equal(marked_vertices.begin(),
1449  marked_vertices.end(),
1450  tria.get_used_vertices().begin(),
1451  [](bool p, bool q) { return !p || q; }),
1452  ExcMessage(
1453  "marked_vertices should be a subset of used vertices in the triangulation "
1454  "but marked_vertices contains one or more vertices that are not used vertices!"));
1455 
1456  // Remove from the map unwanted elements.
1457  if (marked_vertices.size() != 0)
1458  for (auto it = vertices.begin(); it != vertices.end();)
1459  {
1460  if (marked_vertices[it->first] == false)
1461  {
1462  it = vertices.erase(it);
1463  }
1464  else
1465  {
1466  ++it;
1467  }
1468  }
1469 
1470  return find_closest_vertex(vertices, p);
1471  }
1472 
1473 
1474 
1475  template <int dim, template <int, int> class MeshType, int spacedim>
1476 #ifndef _MSC_VER
1477  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1478 #else
1479  std::vector<
1480  typename ::internal::
1481  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1482 #endif
1483  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &mesh,
1484  const unsigned int vertex)
1485  {
1486  // make sure that the given vertex is
1487  // an active vertex of the underlying
1488  // triangulation
1489  AssertIndexRange(vertex, mesh.get_triangulation().n_vertices());
1490  Assert(mesh.get_triangulation().get_used_vertices()[vertex],
1491  ExcVertexNotUsed(vertex));
1492 
1493  // use a set instead of a vector
1494  // to ensure that cells are inserted only
1495  // once
1496  std::set<typename ::internal::
1497  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1498  adjacent_cells;
1499 
1500  // go through all active cells and look if the vertex is part of that cell
1501  //
1502  // in 1d, this is all we need to care about. in 2d/3d we also need to worry
1503  // that the vertex might be a hanging node on a face or edge of a cell; in
1504  // this case, we would want to add those cells as well on whose faces the
1505  // vertex is located but for which it is not a vertex itself.
1506  //
1507  // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging
1508  // node can only be in the middle of a face and we can query the neighboring
1509  // cell from the current cell. on the other hand, in 3d a hanging node
1510  // vertex can also be on an edge but there can be many other cells on
1511  // this edge and we can not access them from the cell we are currently
1512  // on.
1513  //
1514  // so, in the 3d case, if we run the algorithm as in 2d, we catch all
1515  // those cells for which the vertex we seek is on a *subface*, but we
1516  // miss the case of cells for which the vertex we seek is on a
1517  // sub-edge for which there is no corresponding sub-face (because the
1518  // immediate neighbor behind this face is not refined), see for example
1519  // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we
1520  // haven't yet found the vertex for the current cell we also need to
1521  // look at the mid-points of edges
1522  //
1523  // as a final note, deciding whether a neighbor is actually coarser is
1524  // simple in the case of isotropic refinement (we just need to look at
1525  // the level of the current and the neighboring cell). however, this
1526  // isn't so simple if we have used anisotropic refinement since then
1527  // the level of a cell is not indicative of whether it is coarser or
1528  // not than the current cell. ultimately, we want to add all cells on
1529  // which the vertex is, independent of whether they are coarser or
1530  // finer and so in the 2d case below we simply add *any* *active* neighbor.
1531  // in the worst case, we add cells multiple times to the adjacent_cells
1532  // list, but std::set throws out those cells already entered
1533  for (const auto &cell : mesh.active_cell_iterators())
1534  {
1535  for (const unsigned int v : cell->vertex_indices())
1536  if (cell->vertex_index(v) == vertex)
1537  {
1538  // OK, we found a cell that contains
1539  // the given vertex. We add it
1540  // to the list.
1541  adjacent_cells.insert(cell);
1542 
1543  // as explained above, in 2+d we need to check whether
1544  // this vertex is on a face behind which there is a
1545  // (possibly) coarser neighbor. if this is the case,
1546  // then we need to also add this neighbor
1547  if (dim >= 2)
1548  for (unsigned int vface = 0; vface < dim; vface++)
1549  {
1550  const unsigned int face =
1551  GeometryInfo<dim>::vertex_to_face[v][vface]; // TODO
1552 
1553  if (!cell->at_boundary(face) &&
1554  cell->neighbor(face)->is_active())
1555  {
1556  // there is a (possibly) coarser cell behind a
1557  // face to which the vertex belongs. the
1558  // vertex we are looking at is then either a
1559  // vertex of that coarser neighbor, or it is a
1560  // hanging node on one of the faces of that
1561  // cell. in either case, it is adjacent to the
1562  // vertex, so add it to the list as well (if
1563  // the cell was already in the list then the
1564  // std::set makes sure that we get it only
1565  // once)
1566  adjacent_cells.insert(cell->neighbor(face));
1567  }
1568  }
1569 
1570  // in any case, we have found a cell, so go to the next cell
1571  goto next_cell;
1572  }
1573 
1574  // in 3d also loop over the edges
1575  if (dim >= 3)
1576  {
1577  for (unsigned int e = 0; e < cell->n_lines(); ++e)
1578  if (cell->line(e)->has_children())
1579  // the only place where this vertex could have been
1580  // hiding is on the mid-edge point of the edge we
1581  // are looking at
1582  if (cell->line(e)->child(0)->vertex_index(1) == vertex)
1583  {
1584  adjacent_cells.insert(cell);
1585 
1586  // jump out of this tangle of nested loops
1587  goto next_cell;
1588  }
1589  }
1590 
1591  // in more than 3d we would probably have to do the same as
1592  // above also for even lower-dimensional objects
1593  Assert(dim <= 3, ExcNotImplemented());
1594 
1595  // move on to the next cell if we have found the
1596  // vertex on the current one
1597  next_cell:;
1598  }
1599 
1600  // if this was an active vertex then there needs to have been
1601  // at least one cell to which it is adjacent!
1602  Assert(adjacent_cells.size() > 0, ExcInternalError());
1603 
1604  // return the result as a vector, rather than the set we built above
1605  return std::vector<
1606  typename ::internal::
1607  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
1608  adjacent_cells.begin(), adjacent_cells.end());
1609  }
1610 
1611 
1612 
1613  template <int dim, int spacedim>
1614  std::vector<std::vector<Tensor<1, spacedim>>>
1616  const Triangulation<dim, spacedim> &mesh,
1617  const std::vector<
1619  &vertex_to_cells)
1620  {
1621  const std::vector<Point<spacedim>> &vertices = mesh.get_vertices();
1622  const unsigned int n_vertices = vertex_to_cells.size();
1623 
1624  AssertDimension(vertices.size(), n_vertices);
1625 
1626 
1627  std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
1628  n_vertices);
1629  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1630  if (mesh.vertex_used(vertex))
1631  {
1632  const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
1633  vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
1634 
1635  typename std::set<typename Triangulation<dim, spacedim>::
1636  active_cell_iterator>::iterator it =
1637  vertex_to_cells[vertex].begin();
1638  for (unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
1639  {
1640  vertex_to_cell_centers[vertex][cell] =
1641  (*it)->center() - vertices[vertex];
1642  vertex_to_cell_centers[vertex][cell] /=
1643  vertex_to_cell_centers[vertex][cell].norm();
1644  }
1645  }
1646  return vertex_to_cell_centers;
1647  }
1648 
1649 
1650  namespace internal
1651  {
1652  template <int spacedim>
1653  bool
1655  const unsigned int a,
1656  const unsigned int b,
1657  const Tensor<1, spacedim> & point_direction,
1658  const std::vector<Tensor<1, spacedim>> &center_directions)
1659  {
1660  const double scalar_product_a = center_directions[a] * point_direction;
1661  const double scalar_product_b = center_directions[b] * point_direction;
1662 
1663  // The function is supposed to return if a is before b. We are looking
1664  // for the alignment of point direction and center direction, therefore
1665  // return if the scalar product of a is larger.
1666  return (scalar_product_a > scalar_product_b);
1667  }
1668  } // namespace internal
1669 
1670  template <int dim, template <int, int> class MeshType, int spacedim>
1671 #ifndef _MSC_VER
1672  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1673 #else
1674  std::pair<typename ::internal::
1675  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1676  Point<dim>>
1677 #endif
1679  const Mapping<dim, spacedim> & mapping,
1680  const MeshType<dim, spacedim> &mesh,
1681  const Point<spacedim> & p,
1682  const std::vector<
1683  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1684  & vertex_to_cells,
1685  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1686  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
1687  const std::vector<bool> & marked_vertices,
1688  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree,
1689  const double tolerance)
1690  {
1691  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1692  Point<dim>>
1693  cell_and_position;
1694  // To handle points at the border we keep track of points which are close to
1695  // the unit cell:
1696  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1697  Point<dim>>
1698  cell_and_position_approx;
1699 
1700  bool found_cell = false;
1701  bool approx_cell = false;
1702 
1703  unsigned int closest_vertex_index = 0;
1704  Tensor<1, spacedim> vertex_to_point;
1705  auto current_cell = cell_hint;
1706 
1707  while (found_cell == false)
1708  {
1709  // First look at the vertices of the cell cell_hint. If it's an
1710  // invalid cell, then query for the closest global vertex
1711  if (current_cell.state() == IteratorState::valid)
1712  {
1713  const auto cell_vertices = mapping.get_vertices(current_cell);
1714  const unsigned int closest_vertex =
1715  find_closest_vertex_of_cell<dim, spacedim>(current_cell,
1716  p,
1717  mapping);
1718  vertex_to_point = p - cell_vertices[closest_vertex];
1719  closest_vertex_index = current_cell->vertex_index(closest_vertex);
1720  }
1721  else
1722  {
1723  if (!used_vertices_rtree.empty())
1724  {
1725  // If we have an rtree at our disposal, use it.
1726  using ValueType = std::pair<Point<spacedim>, unsigned int>;
1727  std::function<bool(const ValueType &)> marked;
1728  if (marked_vertices.size() == mesh.n_vertices())
1729  marked = [&marked_vertices](const ValueType &value) -> bool {
1730  return marked_vertices[value.second];
1731  };
1732  else
1733  marked = [](const ValueType &) -> bool { return true; };
1734 
1735  std::vector<std::pair<Point<spacedim>, unsigned int>> res;
1736  used_vertices_rtree.query(
1737  boost::geometry::index::nearest(p, 1) &&
1738  boost::geometry::index::satisfies(marked),
1739  std::back_inserter(res));
1740 
1741  // We should have one and only one result
1742  AssertDimension(res.size(), 1);
1743  closest_vertex_index = res[0].second;
1744  }
1745  else
1746  {
1747  closest_vertex_index = GridTools::find_closest_vertex(
1748  mapping, mesh, p, marked_vertices);
1749  }
1750  vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
1751  }
1752 
1753  const double vertex_point_norm = vertex_to_point.norm();
1754  if (vertex_point_norm > 0)
1755  vertex_to_point /= vertex_point_norm;
1756 
1757  const unsigned int n_neighbor_cells =
1758  vertex_to_cells[closest_vertex_index].size();
1759 
1760  // Create a corresponding map of vectors from vertex to cell center
1761  std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
1762 
1763  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
1764  neighbor_permutation[i] = i;
1765 
1766  auto comp = [&](const unsigned int a, const unsigned int b) -> bool {
1767  return internal::compare_point_association<spacedim>(
1768  a,
1769  b,
1770  vertex_to_point,
1771  vertex_to_cell_centers[closest_vertex_index]);
1772  };
1773 
1774  std::sort(neighbor_permutation.begin(),
1775  neighbor_permutation.end(),
1776  comp);
1777  // It is possible the vertex is close
1778  // to an edge, thus we add a tolerance
1779  // to keep also the "best" cell
1780  double best_distance = tolerance;
1781 
1782  // Search all of the cells adjacent to the closest vertex of the cell
1783  // hint Most likely we will find the point in them.
1784  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
1785  {
1786  try
1787  {
1788  auto cell = vertex_to_cells[closest_vertex_index].begin();
1789  std::advance(cell, neighbor_permutation[i]);
1790 
1791  if (!(*cell)->is_artificial())
1792  {
1793  const Point<dim> p_unit =
1794  mapping.transform_real_to_unit_cell(*cell, p);
1796  tolerance))
1797  {
1798  cell_and_position.first = *cell;
1799  cell_and_position.second = p_unit;
1800  found_cell = true;
1801  approx_cell = false;
1802  break;
1803  }
1804  // The point is not inside this cell: checking how far
1805  // outside it is and whether we want to use this cell as a
1806  // backup if we can't find a cell within which the point
1807  // lies.
1808  const double dist =
1810  if (dist < best_distance)
1811  {
1812  best_distance = dist;
1813  cell_and_position_approx.first = *cell;
1814  cell_and_position_approx.second = p_unit;
1815  approx_cell = true;
1816  }
1817  }
1818  }
1819  catch (typename Mapping<dim>::ExcTransformationFailed &)
1820  {}
1821  }
1822 
1823  if (found_cell == true)
1824  return cell_and_position;
1825  else if (approx_cell == true)
1826  return cell_and_position_approx;
1827 
1828  // The first time around, we check for vertices in the hint_cell. If
1829  // that does not work, we set the cell iterator to an invalid one, and
1830  // look for a global vertex close to the point. If that does not work,
1831  // we are in trouble, and just throw an exception.
1832  //
1833  // If we got here, then we did not find the point. If the
1834  // current_cell.state() here is not IteratorState::valid, it means that
1835  // the user did not provide a hint_cell, and at the beginning of the
1836  // while loop we performed an actual global search on the mesh
1837  // vertices. Not finding the point then means the point is outside the
1838  // domain, or that we've had problems with the algorithm above. Try as a
1839  // last resort the other (simpler) algorithm.
1840  if (current_cell.state() != IteratorState::valid)
1842  mapping, mesh, p, marked_vertices, tolerance);
1843 
1844  current_cell = typename MeshType<dim, spacedim>::active_cell_iterator();
1845  }
1846  return cell_and_position;
1847  }
1848 
1849 
1850 
1851  template <int dim, int spacedim>
1852  unsigned int
1855  const Point<spacedim> & position,
1856  const Mapping<dim, spacedim> & mapping)
1857  {
1858  const auto vertices = mapping.get_vertices(cell);
1859  double minimum_distance = position.distance_square(vertices[0]);
1860  unsigned int closest_vertex = 0;
1861 
1862  for (unsigned int v = 1; v < cell->n_vertices(); ++v)
1863  {
1864  const double vertex_distance = position.distance_square(vertices[v]);
1865  if (vertex_distance < minimum_distance)
1866  {
1867  closest_vertex = v;
1868  minimum_distance = vertex_distance;
1869  }
1870  }
1871  return closest_vertex;
1872  }
1873 
1874 
1875 
1876  namespace internal
1877  {
1878  namespace BoundingBoxPredicate
1879  {
1880  template <class MeshType>
1881  std::tuple<BoundingBox<MeshType::space_dimension>, bool>
1883  const typename MeshType::cell_iterator &parent_cell,
1884  const std::function<
1885  bool(const typename MeshType::active_cell_iterator &)> &predicate)
1886  {
1887  bool has_predicate =
1888  false; // Start assuming there's no cells with predicate inside
1889  std::vector<typename MeshType::active_cell_iterator> active_cells;
1890  if (parent_cell->is_active())
1891  active_cells = {parent_cell};
1892  else
1893  // Finding all active cells descendants of the current one (or the
1894  // current one if it is active)
1895  active_cells = get_active_child_cells<MeshType>(parent_cell);
1896 
1897  const unsigned int spacedim = MeshType::space_dimension;
1898 
1899  // Looking for the first active cell which has the property predicate
1900  unsigned int i = 0;
1901  while (i < active_cells.size() && !predicate(active_cells[i]))
1902  ++i;
1903 
1904  // No active cells or no active cells with property
1905  if (active_cells.size() == 0 || i == active_cells.size())
1906  {
1907  BoundingBox<spacedim> bbox;
1908  return std::make_tuple(bbox, has_predicate);
1909  }
1910 
1911  // The two boundary points defining the boundary box
1912  Point<spacedim> maxp = active_cells[i]->vertex(0);
1913  Point<spacedim> minp = active_cells[i]->vertex(0);
1914 
1915  for (; i < active_cells.size(); ++i)
1916  if (predicate(active_cells[i]))
1917  for (const unsigned int v : active_cells[i]->vertex_indices())
1918  for (unsigned int d = 0; d < spacedim; ++d)
1919  {
1920  minp[d] = std::min(minp[d], active_cells[i]->vertex(v)[d]);
1921  maxp[d] = std::max(maxp[d], active_cells[i]->vertex(v)[d]);
1922  }
1923 
1924  has_predicate = true;
1925  BoundingBox<spacedim> bbox(std::make_pair(minp, maxp));
1926  return std::make_tuple(bbox, has_predicate);
1927  }
1928  } // namespace BoundingBoxPredicate
1929  } // namespace internal
1930 
1931 
1932 
1933  template <class MeshType>
1934  std::vector<BoundingBox<MeshType::space_dimension>>
1936  const MeshType &mesh,
1937  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1938  & predicate,
1939  const unsigned int refinement_level,
1940  const bool allow_merge,
1941  const unsigned int max_boxes)
1942  {
1943  // Algorithm brief description: begin with creating bounding boxes of all
1944  // cells at refinement_level (and coarser levels if there are active cells)
1945  // which have the predicate property. These are then merged
1946 
1947  Assert(
1948  refinement_level <= mesh.n_levels(),
1949  ExcMessage(
1950  "Error: refinement level is higher then total levels in the triangulation!"));
1951 
1952  const unsigned int spacedim = MeshType::space_dimension;
1953  std::vector<BoundingBox<spacedim>> bounding_boxes;
1954 
1955  // Creating a bounding box for all active cell on coarser level
1956 
1957  for (unsigned int i = 0; i < refinement_level; ++i)
1958  for (const typename MeshType::cell_iterator &cell :
1959  mesh.active_cell_iterators_on_level(i))
1960  {
1961  bool has_predicate = false;
1962  BoundingBox<spacedim> bbox;
1963  std::tie(bbox, has_predicate) =
1965  MeshType>(cell, predicate);
1966  if (has_predicate)
1967  bounding_boxes.push_back(bbox);
1968  }
1969 
1970  // Creating a Bounding Box for all cells on the chosen refinement_level
1971  for (const typename MeshType::cell_iterator &cell :
1972  mesh.cell_iterators_on_level(refinement_level))
1973  {
1974  bool has_predicate = false;
1975  BoundingBox<spacedim> bbox;
1976  std::tie(bbox, has_predicate) =
1978  MeshType>(cell, predicate);
1979  if (has_predicate)
1980  bounding_boxes.push_back(bbox);
1981  }
1982 
1983  if (!allow_merge)
1984  // If merging is not requested return the created bounding_boxes
1985  return bounding_boxes;
1986  else
1987  {
1988  // Merging part of the algorithm
1989  // Part 1: merging neighbors
1990  // This array stores the indices of arrays we have already merged
1991  std::vector<unsigned int> merged_boxes_idx;
1992  bool found_neighbors = true;
1993 
1994  // We merge only neighbors which can be expressed by a single bounding
1995  // box e.g. in 1d [0,1] and [1,2] can be described with [0,2] without
1996  // losing anything
1997  while (found_neighbors)
1998  {
1999  found_neighbors = false;
2000  for (unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
2001  {
2002  if (std::find(merged_boxes_idx.begin(),
2003  merged_boxes_idx.end(),
2004  i) == merged_boxes_idx.end())
2005  for (unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
2006  if (std::find(merged_boxes_idx.begin(),
2007  merged_boxes_idx.end(),
2008  j) == merged_boxes_idx.end() &&
2009  bounding_boxes[i].get_neighbor_type(
2010  bounding_boxes[j]) ==
2012  {
2013  bounding_boxes[i].merge_with(bounding_boxes[j]);
2014  merged_boxes_idx.push_back(j);
2015  found_neighbors = true;
2016  }
2017  }
2018  }
2019 
2020  // Copying the merged boxes into merged_b_boxes
2021  std::vector<BoundingBox<spacedim>> merged_b_boxes;
2022  for (unsigned int i = 0; i < bounding_boxes.size(); ++i)
2023  if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
2024  merged_boxes_idx.end())
2025  merged_b_boxes.push_back(bounding_boxes[i]);
2026 
2027  // Part 2: if there are too many bounding boxes, merging smaller boxes
2028  // This has sense only in dimension 2 or greater, since in dimension 1,
2029  // neighboring intervals can always be merged without problems
2030  if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
2031  {
2032  std::vector<double> volumes;
2033  for (unsigned int i = 0; i < merged_b_boxes.size(); ++i)
2034  volumes.push_back(merged_b_boxes[i].volume());
2035 
2036  while (merged_b_boxes.size() > max_boxes)
2037  {
2038  unsigned int min_idx =
2039  std::min_element(volumes.begin(), volumes.end()) -
2040  volumes.begin();
2041  volumes.erase(volumes.begin() + min_idx);
2042  // Finding a neighbor
2043  bool not_removed = true;
2044  for (unsigned int i = 0;
2045  i < merged_b_boxes.size() && not_removed;
2046  ++i)
2047  // We merge boxes if we have "attached" or "mergeable"
2048  // neighbors, even though mergeable should be dealt with in
2049  // Part 1
2050  if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
2051  merged_b_boxes[min_idx]) ==
2053  merged_b_boxes[i].get_neighbor_type(
2054  merged_b_boxes[min_idx]) ==
2056  {
2057  merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
2058  merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
2059  not_removed = false;
2060  }
2061  Assert(!not_removed,
2062  ExcMessage("Error: couldn't merge bounding boxes!"));
2063  }
2064  }
2065  Assert(merged_b_boxes.size() <= max_boxes,
2066  ExcMessage(
2067  "Error: couldn't reach target number of bounding boxes!"));
2068  return merged_b_boxes;
2069  }
2070  }
2071 
2072 
2073 
2074  template <int spacedim>
2075 #ifndef DOXYGEN
2076  std::tuple<std::vector<std::vector<unsigned int>>,
2077  std::map<unsigned int, unsigned int>,
2078  std::map<unsigned int, std::vector<unsigned int>>>
2079 #else
2080  return_type
2081 #endif
2083  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
2084  const std::vector<Point<spacedim>> & points)
2085  {
2086  unsigned int n_procs = global_bboxes.size();
2087  std::vector<std::vector<unsigned int>> point_owners(n_procs);
2088  std::map<unsigned int, unsigned int> map_owners_found;
2089  std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
2090 
2091  unsigned int n_points = points.size();
2092  for (unsigned int pt = 0; pt < n_points; ++pt)
2093  {
2094  // Keep track of how many processes we guess to own the point
2095  std::vector<unsigned int> owners_found;
2096  // Check in which other processes the point might be
2097  for (unsigned int rk = 0; rk < n_procs; ++rk)
2098  {
2099  for (const BoundingBox<spacedim> &bbox : global_bboxes[rk])
2100  if (bbox.point_inside(points[pt]))
2101  {
2102  point_owners[rk].emplace_back(pt);
2103  owners_found.emplace_back(rk);
2104  break; // We can check now the next process
2105  }
2106  }
2107  Assert(owners_found.size() > 0,
2108  ExcMessage("No owners found for the point " +
2109  std::to_string(pt)));
2110  if (owners_found.size() == 1)
2111  map_owners_found[pt] = owners_found[0];
2112  else
2113  // Multiple owners
2114  map_owners_guessed[pt] = owners_found;
2115  }
2116 
2117  return std::make_tuple(std::move(point_owners),
2118  std::move(map_owners_found),
2119  std::move(map_owners_guessed));
2120  }
2121 
2122  template <int spacedim>
2123 #ifndef DOXYGEN
2124  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2125  std::map<unsigned int, unsigned int>,
2126  std::map<unsigned int, std::vector<unsigned int>>>
2127 #else
2128  return_type
2129 #endif
2131  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
2132  const std::vector<Point<spacedim>> & points)
2133  {
2134  std::map<unsigned int, std::vector<unsigned int>> point_owners;
2135  std::map<unsigned int, unsigned int> map_owners_found;
2136  std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
2137  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> search_result;
2138 
2139  unsigned int n_points = points.size();
2140  for (unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
2141  {
2142  search_result.clear(); // clearing last output
2143 
2144  // Running tree search
2145  covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
2146  std::back_inserter(search_result));
2147 
2148  // Keep track of how many processes we guess to own the point
2149  std::set<unsigned int> owners_found;
2150  // Check in which other processes the point might be
2151  for (const auto &rank_bbox : search_result)
2152  {
2153  // Try to add the owner to the owners found,
2154  // and check if it was already present
2155  const bool pt_inserted = owners_found.insert(pt_n).second;
2156  if (pt_inserted)
2157  point_owners[rank_bbox.second].emplace_back(pt_n);
2158  }
2159  Assert(owners_found.size() > 0,
2160  ExcMessage("No owners found for the point " +
2161  std::to_string(pt_n)));
2162  if (owners_found.size() == 1)
2163  map_owners_found[pt_n] = *owners_found.begin();
2164  else
2165  // Multiple owners
2166  std::copy(owners_found.begin(),
2167  owners_found.end(),
2168  std::back_inserter(map_owners_guessed[pt_n]));
2169  }
2170 
2171  return std::make_tuple(std::move(point_owners),
2172  std::move(map_owners_found),
2173  std::move(map_owners_guessed));
2174  }
2175 
2176 
2177  template <int dim, int spacedim>
2178  std::vector<
2179  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2181  {
2182  std::vector<
2183  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2184  vertex_to_cell_map(triangulation.n_vertices());
2186  cell = triangulation.begin_active(),
2187  endc = triangulation.end();
2188  for (; cell != endc; ++cell)
2189  for (const unsigned int i : cell->vertex_indices())
2190  vertex_to_cell_map[cell->vertex_index(i)].insert(cell);
2191 
2192  // Take care of hanging nodes
2193  cell = triangulation.begin_active();
2194  for (; cell != endc; ++cell)
2195  {
2196  for (unsigned int i : cell->face_indices())
2197  {
2198  if ((cell->at_boundary(i) == false) &&
2199  (cell->neighbor(i)->is_active()))
2200  {
2202  adjacent_cell = cell->neighbor(i);
2203  for (unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
2204  vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert(
2205  adjacent_cell);
2206  }
2207  }
2208 
2209  // in 3d also loop over the edges
2210  if (dim == 3)
2211  {
2212  for (unsigned int i = 0; i < cell->n_lines(); ++i)
2213  if (cell->line(i)->has_children())
2214  // the only place where this vertex could have been
2215  // hiding is on the mid-edge point of the edge we
2216  // are looking at
2217  vertex_to_cell_map[cell->line(i)->child(0)->vertex_index(1)]
2218  .insert(cell);
2219  }
2220  }
2221 
2222  return vertex_to_cell_map;
2223  }
2224 
2225 
2226 
2227  template <int dim, int spacedim>
2228  std::map<unsigned int, types::global_vertex_index>
2231  {
2232  std::map<unsigned int, types::global_vertex_index>
2233  local_to_global_vertex_index;
2234 
2235 #ifndef DEAL_II_WITH_MPI
2236 
2237  // without MPI, this function doesn't make sense because on cannot
2238  // use parallel::distributed::Triangulation in any meaningful
2239  // way
2240  (void)triangulation;
2241  Assert(false,
2242  ExcMessage("This function does not make any sense "
2243  "for parallel::distributed::Triangulation "
2244  "objects if you do not have MPI enabled."));
2245 
2246 #else
2247 
2248  using active_cell_iterator =
2250  const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
2251  vertex_to_cell_map(triangulation);
2252 
2253  // Create a local index for the locally "owned" vertices
2254  types::global_vertex_index next_index = 0;
2255  unsigned int max_cellid_size = 0;
2256  std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
2257  vertices_added;
2258  std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
2259  std::map<types::subdomain_id,
2260  std::vector<std::tuple<types::global_vertex_index,
2262  std::string>>>
2263  vertices_to_send;
2264  active_cell_iterator cell = triangulation.begin_active(),
2265  endc = triangulation.end();
2266  std::set<active_cell_iterator> missing_vert_cells;
2267  std::set<unsigned int> used_vertex_index;
2268  for (; cell != endc; ++cell)
2269  {
2270  if (cell->is_locally_owned())
2271  {
2272  for (const unsigned int i : cell->vertex_indices())
2273  {
2274  types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
2275  typename std::set<active_cell_iterator>::iterator
2276  adjacent_cell = vertex_to_cell[cell->vertex_index(i)].begin(),
2277  end_adj_cell = vertex_to_cell[cell->vertex_index(i)].end();
2278  for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
2279  lowest_subdomain_id =
2280  std::min(lowest_subdomain_id,
2281  (*adjacent_cell)->subdomain_id());
2282 
2283  // See if I "own" this vertex
2284  if (lowest_subdomain_id == cell->subdomain_id())
2285  {
2286  // Check that the vertex we are working on a vertex that has
2287  // not be dealt with yet
2288  if (used_vertex_index.find(cell->vertex_index(i)) ==
2289  used_vertex_index.end())
2290  {
2291  // Set the local index
2292  local_to_global_vertex_index[cell->vertex_index(i)] =
2293  next_index++;
2294 
2295  // Store the information that will be sent to the
2296  // adjacent cells on other subdomains
2297  adjacent_cell =
2298  vertex_to_cell[cell->vertex_index(i)].begin();
2299  for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
2300  if ((*adjacent_cell)->subdomain_id() !=
2301  cell->subdomain_id())
2302  {
2303  std::pair<types::subdomain_id,
2304  types::global_vertex_index>
2305  tmp((*adjacent_cell)->subdomain_id(),
2306  cell->vertex_index(i));
2307  if (vertices_added.find(tmp) ==
2308  vertices_added.end())
2309  {
2310  vertices_to_send[(*adjacent_cell)
2311  ->subdomain_id()]
2312  .emplace_back(i,
2313  cell->vertex_index(i),
2314  cell->id().to_string());
2315  if (cell->id().to_string().size() >
2316  max_cellid_size)
2317  max_cellid_size =
2318  cell->id().to_string().size();
2319  vertices_added.insert(tmp);
2320  }
2321  }
2322  used_vertex_index.insert(cell->vertex_index(i));
2323  }
2324  }
2325  else
2326  {
2327  // We don't own the vertex so we will receive its global
2328  // index
2329  vertices_to_recv[lowest_subdomain_id].insert(
2330  cell->vertex_index(i));
2331  missing_vert_cells.insert(cell);
2332  }
2333  }
2334  }
2335 
2336  // Some hanging nodes are vertices of ghost cells. They need to be
2337  // received.
2338  if (cell->is_ghost())
2339  {
2340  for (unsigned int i : cell->face_indices())
2341  {
2342  if (cell->at_boundary(i) == false)
2343  {
2344  if (cell->neighbor(i)->is_active())
2345  {
2346  typename Triangulation<dim,
2347  spacedim>::active_cell_iterator
2348  adjacent_cell = cell->neighbor(i);
2349  if ((adjacent_cell->is_locally_owned()))
2350  {
2351  types::subdomain_id adj_subdomain_id =
2352  adjacent_cell->subdomain_id();
2353  if (cell->subdomain_id() < adj_subdomain_id)
2354  for (unsigned int j = 0;
2355  j < cell->face(i)->n_vertices();
2356  ++j)
2357  {
2358  vertices_to_recv[cell->subdomain_id()].insert(
2359  cell->face(i)->vertex_index(j));
2360  missing_vert_cells.insert(cell);
2361  }
2362  }
2363  }
2364  }
2365  }
2366  }
2367  }
2368 
2369  // Get the size of the largest CellID string
2370  max_cellid_size =
2371  Utilities::MPI::max(max_cellid_size, triangulation.get_communicator());
2372 
2373  // Make indices global by getting the number of vertices owned by each
2374  // processors and shifting the indices accordingly
2375  types::global_vertex_index shift = 0;
2376  int ierr = MPI_Exscan(&next_index,
2377  &shift,
2378  1,
2380  MPI_SUM,
2381  triangulation.get_communicator());
2382  AssertThrowMPI(ierr);
2383 
2384  std::map<unsigned int, types::global_vertex_index>::iterator
2385  global_index_it = local_to_global_vertex_index.begin(),
2386  global_index_end = local_to_global_vertex_index.end();
2387  for (; global_index_it != global_index_end; ++global_index_it)
2388  global_index_it->second += shift;
2389 
2390 
2391  const int mpi_tag = Utilities::MPI::internal::Tags::
2393  const int mpi_tag2 = Utilities::MPI::internal::Tags::
2395 
2396 
2397  // In a first message, send the global ID of the vertices and the local
2398  // positions in the cells. In a second messages, send the cell ID as a
2399  // resize string. This is done in two messages so that types are not mixed
2400 
2401  // Send the first message
2402  std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
2403  vertices_to_send.size());
2404  std::vector<MPI_Request> first_requests(vertices_to_send.size());
2405  typename std::map<types::subdomain_id,
2406  std::vector<std::tuple<types::global_vertex_index,
2408  std::string>>>::iterator
2409  vert_to_send_it = vertices_to_send.begin(),
2410  vert_to_send_end = vertices_to_send.end();
2411  for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
2412  ++vert_to_send_it, ++i)
2413  {
2414  int destination = vert_to_send_it->first;
2415  const unsigned int n_vertices = vert_to_send_it->second.size();
2416  const int buffer_size = 2 * n_vertices;
2417  vertices_send_buffers[i].resize(buffer_size);
2418 
2419  // fill the buffer
2420  for (unsigned int j = 0; j < n_vertices; ++j)
2421  {
2422  vertices_send_buffers[i][2 * j] =
2423  std::get<0>(vert_to_send_it->second[j]);
2424  vertices_send_buffers[i][2 * j + 1] =
2425  local_to_global_vertex_index[std::get<1>(
2426  vert_to_send_it->second[j])];
2427  }
2428 
2429  // Send the message
2430  ierr = MPI_Isend(vertices_send_buffers[i].data(),
2431  buffer_size,
2433  destination,
2434  mpi_tag,
2435  triangulation.get_communicator(),
2436  &first_requests[i]);
2437  AssertThrowMPI(ierr);
2438  }
2439 
2440  // Receive the first message
2441  std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
2442  vertices_to_recv.size());
2443  typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
2444  vert_to_recv_it = vertices_to_recv.begin(),
2445  vert_to_recv_end = vertices_to_recv.end();
2446  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2447  ++vert_to_recv_it, ++i)
2448  {
2449  int source = vert_to_recv_it->first;
2450  const unsigned int n_vertices = vert_to_recv_it->second.size();
2451  const int buffer_size = 2 * n_vertices;
2452  vertices_recv_buffers[i].resize(buffer_size);
2453 
2454  // Receive the message
2455  ierr = MPI_Recv(vertices_recv_buffers[i].data(),
2456  buffer_size,
2458  source,
2459  mpi_tag,
2460  triangulation.get_communicator(),
2461  MPI_STATUS_IGNORE);
2462  AssertThrowMPI(ierr);
2463  }
2464 
2465 
2466  // Send second message
2467  std::vector<std::vector<char>> cellids_send_buffers(
2468  vertices_to_send.size());
2469  std::vector<MPI_Request> second_requests(vertices_to_send.size());
2470  vert_to_send_it = vertices_to_send.begin();
2471  for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
2472  ++vert_to_send_it, ++i)
2473  {
2474  int destination = vert_to_send_it->first;
2475  const unsigned int n_vertices = vert_to_send_it->second.size();
2476  const int buffer_size = max_cellid_size * n_vertices;
2477  cellids_send_buffers[i].resize(buffer_size);
2478 
2479  // fill the buffer
2480  unsigned int pos = 0;
2481  for (unsigned int j = 0; j < n_vertices; ++j)
2482  {
2483  std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
2484  for (unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
2485  {
2486  if (k < cell_id.size())
2487  cellids_send_buffers[i][pos] = cell_id[k];
2488  // if necessary fill up the reserved part of the buffer with an
2489  // invalid value
2490  else
2491  cellids_send_buffers[i][pos] = '-';
2492  }
2493  }
2494 
2495  // Send the message
2496  ierr = MPI_Isend(cellids_send_buffers[i].data(),
2497  buffer_size,
2498  MPI_CHAR,
2499  destination,
2500  mpi_tag2,
2501  triangulation.get_communicator(),
2502  &second_requests[i]);
2503  AssertThrowMPI(ierr);
2504  }
2505 
2506  // Receive the second message
2507  std::vector<std::vector<char>> cellids_recv_buffers(
2508  vertices_to_recv.size());
2509  vert_to_recv_it = vertices_to_recv.begin();
2510  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2511  ++vert_to_recv_it, ++i)
2512  {
2513  int source = vert_to_recv_it->first;
2514  const unsigned int n_vertices = vert_to_recv_it->second.size();
2515  const int buffer_size = max_cellid_size * n_vertices;
2516  cellids_recv_buffers[i].resize(buffer_size);
2517 
2518  // Receive the message
2519  ierr = MPI_Recv(cellids_recv_buffers[i].data(),
2520  buffer_size,
2521  MPI_CHAR,
2522  source,
2523  mpi_tag2,
2524  triangulation.get_communicator(),
2525  MPI_STATUS_IGNORE);
2526  AssertThrowMPI(ierr);
2527  }
2528 
2529 
2530  // Match the data received with the required vertices
2531  vert_to_recv_it = vertices_to_recv.begin();
2532  for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2533  ++i, ++vert_to_recv_it)
2534  {
2535  for (unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
2536  {
2537  const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
2538  const types::global_vertex_index global_id_recv =
2539  vertices_recv_buffers[i][2 * j + 1];
2540  const std::string cellid_recv(
2541  &cellids_recv_buffers[i][max_cellid_size * j],
2542  &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
2543  bool found = false;
2544  typename std::set<active_cell_iterator>::iterator
2545  cell_set_it = missing_vert_cells.begin(),
2546  end_cell_set = missing_vert_cells.end();
2547  for (; (found == false) && (cell_set_it != end_cell_set);
2548  ++cell_set_it)
2549  {
2550  typename std::set<active_cell_iterator>::iterator
2551  candidate_cell =
2552  vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
2553  end_cell =
2554  vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
2555  for (; candidate_cell != end_cell; ++candidate_cell)
2556  {
2557  std::string current_cellid =
2558  (*candidate_cell)->id().to_string();
2559  current_cellid.resize(max_cellid_size, '-');
2560  if (current_cellid.compare(cellid_recv) == 0)
2561  {
2562  local_to_global_vertex_index
2563  [(*candidate_cell)->vertex_index(local_pos_recv)] =
2564  global_id_recv;
2565  found = true;
2566 
2567  break;
2568  }
2569  }
2570  }
2571  }
2572  }
2573 #endif
2574 
2575  return local_to_global_vertex_index;
2576  }
2577 
2578 
2579 
2580  template <int dim, int spacedim>
2581  void
2584  DynamicSparsityPattern & cell_connectivity)
2585  {
2586  cell_connectivity.reinit(triangulation.n_active_cells(),
2587  triangulation.n_active_cells());
2588 
2589  // loop over all cells and their neighbors to build the sparsity
2590  // pattern. note that it's a bit hard to enter all the connections when a
2591  // neighbor has children since we would need to find out which of its
2592  // children is adjacent to the current cell. this problem can be omitted
2593  // if we only do something if the neighbor has no children -- in that case
2594  // it is either on the same or a coarser level than we are. in return, we
2595  // have to add entries in both directions for both cells
2596  for (const auto &cell : triangulation.active_cell_iterators())
2597  {
2598  const unsigned int index = cell->active_cell_index();
2599  cell_connectivity.add(index, index);
2600  for (auto f : cell->face_indices())
2601  if ((cell->at_boundary(f) == false) &&
2602  (cell->neighbor(f)->has_children() == false))
2603  {
2604  const unsigned int other_index =
2605  cell->neighbor(f)->active_cell_index();
2606  cell_connectivity.add(index, other_index);
2607  cell_connectivity.add(other_index, index);
2608  }
2609  }
2610  }
2611 
2612 
2613 
2614  template <int dim, int spacedim>
2615  void
2618  DynamicSparsityPattern & cell_connectivity)
2619  {
2620  std::vector<std::vector<unsigned int>> vertex_to_cell(
2621  triangulation.n_vertices());
2622  for (const auto &cell : triangulation.active_cell_iterators())
2623  {
2624  for (const unsigned int v : cell->vertex_indices())
2625  vertex_to_cell[cell->vertex_index(v)].push_back(
2626  cell->active_cell_index());
2627  }
2628 
2629  cell_connectivity.reinit(triangulation.n_active_cells(),
2630  triangulation.n_active_cells());
2631  for (const auto &cell : triangulation.active_cell_iterators())
2632  {
2633  for (const unsigned int v : cell->vertex_indices())
2634  for (unsigned int n = 0;
2635  n < vertex_to_cell[cell->vertex_index(v)].size();
2636  ++n)
2637  cell_connectivity.add(cell->active_cell_index(),
2638  vertex_to_cell[cell->vertex_index(v)][n]);
2639  }
2640  }
2641 
2642 
2643  template <int dim, int spacedim>
2644  void
2647  const unsigned int level,
2648  DynamicSparsityPattern & cell_connectivity)
2649  {
2650  std::vector<std::vector<unsigned int>> vertex_to_cell(
2651  triangulation.n_vertices());
2652  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
2653  triangulation.begin(level);
2654  cell != triangulation.end(level);
2655  ++cell)
2656  {
2657  for (const unsigned int v : cell->vertex_indices())
2658  vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
2659  }
2660 
2661  cell_connectivity.reinit(triangulation.n_cells(level),
2662  triangulation.n_cells(level));
2663  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
2664  triangulation.begin(level);
2665  cell != triangulation.end(level);
2666  ++cell)
2667  {
2668  for (const unsigned int v : cell->vertex_indices())
2669  for (unsigned int n = 0;
2670  n < vertex_to_cell[cell->vertex_index(v)].size();
2671  ++n)
2672  cell_connectivity.add(cell->index(),
2673  vertex_to_cell[cell->vertex_index(v)][n]);
2674  }
2675  }
2676 
2677 
2678 
2679  template <int dim, int spacedim>
2680  void
2681  partition_triangulation(const unsigned int n_partitions,
2683  const SparsityTools::Partitioner partitioner)
2684  {
2686  &triangulation) == nullptr),
2687  ExcMessage("Objects of type parallel::distributed::Triangulation "
2688  "are already partitioned implicitly and can not be "
2689  "partitioned again explicitly."));
2690 
2691  std::vector<unsigned int> cell_weights;
2692 
2693  // Get cell weighting if a signal has been attached to the triangulation
2694  if (!triangulation.signals.cell_weight.empty())
2695  {
2696  cell_weights.resize(triangulation.n_active_cells(), 0U);
2697 
2698  // In a first step, obtain the weights of the locally owned
2699  // cells. For all others, the weight remains at the zero the
2700  // vector was initialized with above.
2701  for (const auto &cell : triangulation.active_cell_iterators())
2702  if (cell->is_locally_owned())
2703  cell_weights[cell->active_cell_index()] =
2704  triangulation.signals.cell_weight(
2706 
2707  // If this is a parallel triangulation, we then need to also
2708  // get the weights for all other cells. We have asserted above
2709  // that this function can't be used for
2710  // parallel::distribute::Triangulation objects, so the only
2711  // ones we have to worry about here are
2712  // parallel::shared::Triangulation
2713  if (const auto shared_tria =
2715  &triangulation))
2716  Utilities::MPI::sum(cell_weights,
2717  shared_tria->get_communicator(),
2718  cell_weights);
2719  }
2720 
2721  // Call the other more general function
2722  partition_triangulation(n_partitions,
2723  cell_weights,
2724  triangulation,
2725  partitioner);
2726  }
2727 
2728 
2729 
2730  template <int dim, int spacedim>
2731  void
2732  partition_triangulation(const unsigned int n_partitions,
2733  const std::vector<unsigned int> &cell_weights,
2735  const SparsityTools::Partitioner partitioner)
2736  {
2738  &triangulation) == nullptr),
2739  ExcMessage("Objects of type parallel::distributed::Triangulation "
2740  "are already partitioned implicitly and can not be "
2741  "partitioned again explicitly."));
2742  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
2743 
2744  // check for an easy return
2745  if (n_partitions == 1)
2746  {
2747  for (const auto &cell : triangulation.active_cell_iterators())
2748  cell->set_subdomain_id(0);
2749  return;
2750  }
2751 
2752  // we decompose the domain by first
2753  // generating the connection graph of all
2754  // cells with their neighbors, and then
2755  // passing this graph off to METIS.
2756  // finally defer to the other function for
2757  // partitioning and assigning subdomain ids
2758  DynamicSparsityPattern cell_connectivity;
2759  get_face_connectivity_of_cells(triangulation, cell_connectivity);
2760 
2761  SparsityPattern sp_cell_connectivity;
2762  sp_cell_connectivity.copy_from(cell_connectivity);
2763  partition_triangulation(n_partitions,
2764  cell_weights,
2765  sp_cell_connectivity,
2766  triangulation,
2767  partitioner);
2768  }
2769 
2770 
2771 
2772  template <int dim, int spacedim>
2773  void
2774  partition_triangulation(const unsigned int n_partitions,
2775  const SparsityPattern & cell_connection_graph,
2777  const SparsityTools::Partitioner partitioner)
2778  {
2780  &triangulation) == nullptr),
2781  ExcMessage("Objects of type parallel::distributed::Triangulation "
2782  "are already partitioned implicitly and can not be "
2783  "partitioned again explicitly."));
2784 
2785  std::vector<unsigned int> cell_weights;
2786 
2787  // Get cell weighting if a signal has been attached to the triangulation
2788  if (!triangulation.signals.cell_weight.empty())
2789  {
2790  cell_weights.resize(triangulation.n_active_cells(), 0U);
2791 
2792  // In a first step, obtain the weights of the locally owned
2793  // cells. For all others, the weight remains at the zero the
2794  // vector was initialized with above.
2795  for (const auto &cell : triangulation.active_cell_iterators())
2796  if (cell->is_locally_owned())
2797  cell_weights[cell->active_cell_index()] =
2798  triangulation.signals.cell_weight(
2800 
2801  // If this is a parallel triangulation, we then need to also
2802  // get the weights for all other cells. We have asserted above
2803  // that this function can't be used for
2804  // parallel::distribute::Triangulation objects, so the only
2805  // ones we have to worry about here are
2806  // parallel::shared::Triangulation
2807  if (const auto shared_tria =
2809  &triangulation))
2810  Utilities::MPI::sum(cell_weights,
2811  shared_tria->get_communicator(),
2812  cell_weights);
2813  }
2814 
2815  // Call the other more general function
2816  partition_triangulation(n_partitions,
2817  cell_weights,
2818  cell_connection_graph,
2819  triangulation,
2820  partitioner);
2821  }
2822 
2823 
2824 
2825  template <int dim, int spacedim>
2826  void
2827  partition_triangulation(const unsigned int n_partitions,
2828  const std::vector<unsigned int> &cell_weights,
2829  const SparsityPattern & cell_connection_graph,
2831  const SparsityTools::Partitioner partitioner)
2832  {
2834  &triangulation) == nullptr),
2835  ExcMessage("Objects of type parallel::distributed::Triangulation "
2836  "are already partitioned implicitly and can not be "
2837  "partitioned again explicitly."));
2838  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
2839  Assert(cell_connection_graph.n_rows() == triangulation.n_active_cells(),
2840  ExcMessage("Connectivity graph has wrong size"));
2841  Assert(cell_connection_graph.n_cols() == triangulation.n_active_cells(),
2842  ExcMessage("Connectivity graph has wrong size"));
2843 
2844  // signal that partitioning is going to happen
2845  triangulation.signals.pre_partition();
2846 
2847  // check for an easy return
2848  if (n_partitions == 1)
2849  {
2850  for (const auto &cell : triangulation.active_cell_iterators())
2851  cell->set_subdomain_id(0);
2852  return;
2853  }
2854 
2855  // partition this connection graph and get
2856  // back a vector of indices, one per degree
2857  // of freedom (which is associated with a
2858  // cell)
2859  std::vector<unsigned int> partition_indices(triangulation.n_active_cells());
2860  SparsityTools::partition(cell_connection_graph,
2861  cell_weights,
2862  n_partitions,
2863  partition_indices,
2864  partitioner);
2865 
2866  // finally loop over all cells and set the subdomain ids
2867  for (const auto &cell : triangulation.active_cell_iterators())
2868  cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
2869  }
2870 
2871 
2872  namespace internal
2873  {
2877  template <class IT>
2878  void
2880  unsigned int & current_proc_idx,
2881  unsigned int & current_cell_idx,
2882  const unsigned int n_active_cells,
2883  const unsigned int n_partitions)
2884  {
2885  if (cell->is_active())
2886  {
2887  while (current_cell_idx >=
2888  std::floor(static_cast<uint_least64_t>(n_active_cells) *
2889  (current_proc_idx + 1) / n_partitions))
2890  ++current_proc_idx;
2891  cell->set_subdomain_id(current_proc_idx);
2892  ++current_cell_idx;
2893  }
2894  else
2895  {
2896  for (unsigned int n = 0; n < cell->n_children(); ++n)
2898  current_proc_idx,
2899  current_cell_idx,
2901  n_partitions);
2902  }
2903  }
2904  } // namespace internal
2905 
2906  template <int dim, int spacedim>
2907  void
2908  partition_triangulation_zorder(const unsigned int n_partitions,
2910  const bool group_siblings)
2911  {
2913  &triangulation) == nullptr),
2914  ExcMessage("Objects of type parallel::distributed::Triangulation "
2915  "are already partitioned implicitly and can not be "
2916  "partitioned again explicitly."));
2917  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
2918 
2919  // signal that partitioning is going to happen
2920  triangulation.signals.pre_partition();
2921 
2922  // check for an easy return
2923  if (n_partitions == 1)
2924  {
2925  for (const auto &cell : triangulation.active_cell_iterators())
2926  cell->set_subdomain_id(0);
2927  return;
2928  }
2929 
2930  // Duplicate the coarse cell reordoring
2931  // as done in p4est
2932  std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
2933  std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
2934 
2935  DynamicSparsityPattern cell_connectivity;
2937  0,
2938  cell_connectivity);
2939  coarse_cell_to_p4est_tree_permutation.resize(triangulation.n_cells(0));
2940  SparsityTools::reorder_hierarchical(cell_connectivity,
2941  coarse_cell_to_p4est_tree_permutation);
2942 
2943  p4est_tree_to_coarse_cell_permutation =
2944  Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
2945 
2946  unsigned int current_proc_idx = 0;
2947  unsigned int current_cell_idx = 0;
2948  const unsigned int n_active_cells = triangulation.n_active_cells();
2949 
2950  // set subdomain id for active cell descendants
2951  // of each coarse cell in permuted order
2952  for (unsigned int idx = 0; idx < triangulation.n_cells(0); ++idx)
2953  {
2954  const unsigned int coarse_cell_idx =
2955  p4est_tree_to_coarse_cell_permutation[idx];
2956  typename Triangulation<dim, spacedim>::cell_iterator coarse_cell(
2957  &triangulation, 0, coarse_cell_idx);
2958 
2960  current_proc_idx,
2961  current_cell_idx,
2962  n_active_cells,
2963  n_partitions);
2964  }
2965 
2966  // if all children of a cell are active (e.g. we
2967  // have a cell that is refined once and no part
2968  // is refined further), p4est places all of them
2969  // on the same processor. The new owner will be
2970  // the processor with the largest number of children
2971  // (ties are broken by picking the lower rank).
2972  // Duplicate this logic here.
2973  if (group_siblings)
2974  {
2976  cell = triangulation.begin(),
2977  endc = triangulation.end();
2978  for (; cell != endc; ++cell)
2979  {
2980  if (cell->is_active())
2981  continue;
2982  bool all_children_active = true;
2983  std::map<unsigned int, unsigned int> map_cpu_n_cells;
2984  for (unsigned int n = 0; n < cell->n_children(); ++n)
2985  if (!cell->child(n)->is_active())
2986  {
2987  all_children_active = false;
2988  break;
2989  }
2990  else
2991  ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
2992 
2993  if (!all_children_active)
2994  continue;
2995 
2996  unsigned int new_owner = cell->child(0)->subdomain_id();
2997  for (std::map<unsigned int, unsigned int>::iterator it =
2998  map_cpu_n_cells.begin();
2999  it != map_cpu_n_cells.end();
3000  ++it)
3001  if (it->second > map_cpu_n_cells[new_owner])
3002  new_owner = it->first;
3003 
3004  for (unsigned int n = 0; n < cell->n_children(); ++n)
3005  cell->child(n)->set_subdomain_id(new_owner);
3006  }
3007  }
3008  }
3009 
3010 
3011  template <int dim, int spacedim>
3012  void
3014  {
3015  unsigned int n_levels = triangulation.n_levels();
3016  for (int lvl = n_levels - 1; lvl >= 0; --lvl)
3017  {
3019  cell = triangulation.begin(lvl),
3020  endc = triangulation.end(lvl);
3021  for (; cell != endc; ++cell)
3022  {
3023  if (cell->is_active())
3024  cell->set_level_subdomain_id(cell->subdomain_id());
3025  else
3026  {
3027  Assert(cell->child(0)->level_subdomain_id() !=
3029  ExcInternalError());
3030  cell->set_level_subdomain_id(
3031  cell->child(0)->level_subdomain_id());
3032  }
3033  }
3034  }
3035  }
3036 
3037 
3038  template <int dim, int spacedim>
3039  void
3041  std::vector<types::subdomain_id> & subdomain)
3042  {
3043  Assert(subdomain.size() == triangulation.n_active_cells(),
3044  ExcDimensionMismatch(subdomain.size(),
3045  triangulation.n_active_cells()));
3046  for (const auto &cell : triangulation.active_cell_iterators())
3047  subdomain[cell->active_cell_index()] = cell->subdomain_id();
3048  }
3049 
3050 
3051 
3052  template <int dim, int spacedim>
3053  unsigned int
3056  const types::subdomain_id subdomain)
3057  {
3058  unsigned int count = 0;
3059  for (const auto &cell : triangulation.active_cell_iterators())
3060  if (cell->subdomain_id() == subdomain)
3061  ++count;
3062 
3063  return count;
3064  }
3065 
3066 
3067 
3068  template <int dim, int spacedim>
3069  std::vector<bool>
3071  {
3072  // start with all vertices
3073  std::vector<bool> locally_owned_vertices =
3074  triangulation.get_used_vertices();
3075 
3076  // if the triangulation is distributed, eliminate those that
3077  // are owned by other processors -- either because the vertex is
3078  // on an artificial cell, or because it is on a ghost cell with
3079  // a smaller subdomain
3082  *>(&triangulation))
3083  for (const auto &cell : triangulation.active_cell_iterators())
3084  if (cell->is_artificial() ||
3085  (cell->is_ghost() &&
3086  (cell->subdomain_id() < tr->locally_owned_subdomain())))
3087  for (const unsigned int v : cell->vertex_indices())
3088  locally_owned_vertices[cell->vertex_index(v)] = false;
3089 
3090  return locally_owned_vertices;
3091  }
3092 
3093 
3094 
3095  namespace internal
3096  {
3097  template <int dim, int spacedim>
3098  double
3100  const Mapping<dim, spacedim> &mapping)
3101  {
3102  // see also TriaAccessor::diameter()
3103 
3104  const auto vertices = mapping.get_vertices(cell);
3105  switch (cell->reference_cell_type())
3106  {
3108  return (vertices[1] - vertices[0]).norm();
3110  return std::max(std::max((vertices[1] - vertices[0]).norm(),
3111  (vertices[2] - vertices[1]).norm()),
3112  (vertices[2] - vertices[0]).norm());
3114  return std::max((vertices[3] - vertices[0]).norm(),
3115  (vertices[2] - vertices[1]).norm());
3117  return std::max(
3118  std::max(std::max((vertices[1] - vertices[0]).norm(),
3119  (vertices[2] - vertices[0]).norm()),
3120  std::max((vertices[2] - vertices[1]).norm(),
3121  (vertices[3] - vertices[0]).norm())),
3122  std::max((vertices[3] - vertices[1]).norm(),
3123  (vertices[3] - vertices[2]).norm()));
3125  return std::max(std::max((vertices[7] - vertices[0]).norm(),
3126  (vertices[6] - vertices[1]).norm()),
3127  std::max((vertices[2] - vertices[5]).norm(),
3128  (vertices[3] - vertices[4]).norm()));
3129  default:
3130  Assert(false, ExcNotImplemented());
3131  return -1e10;
3132  }
3133  }
3134  } // namespace internal
3135 
3136 
3137  template <int dim, int spacedim>
3138  double
3140  const Mapping<dim, spacedim> & mapping)
3141  {
3142  double min_diameter = std::numeric_limits<double>::max();
3143  for (const auto &cell : triangulation.active_cell_iterators())
3144  if (!cell->is_artificial())
3145  min_diameter =
3146  std::min(min_diameter,
3147  internal::diameter<dim, spacedim>(cell, mapping));
3148 
3149  double global_min_diameter = 0;
3150 
3151 #ifdef DEAL_II_WITH_MPI
3152  if (const parallel::TriangulationBase<dim, spacedim> *p_tria =
3153  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3154  &triangulation))
3155  global_min_diameter =
3156  Utilities::MPI::min(min_diameter, p_tria->get_communicator());
3157  else
3158 #endif
3159  global_min_diameter = min_diameter;
3160 
3161  return global_min_diameter;
3162  }
3163 
3164 
3165 
3166  template <int dim, int spacedim>
3167  double
3169  const Mapping<dim, spacedim> & mapping)
3170  {
3171  double max_diameter = 0.;
3172  for (const auto &cell : triangulation.active_cell_iterators())
3173  if (!cell->is_artificial())
3174  max_diameter =
3175  std::max(max_diameter, internal::diameter(cell, mapping));
3176 
3177  double global_max_diameter = 0;
3178 
3179 #ifdef DEAL_II_WITH_MPI
3180  if (const parallel::TriangulationBase<dim, spacedim> *p_tria =
3181  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3182  &triangulation))
3183  global_max_diameter =
3184  Utilities::MPI::max(max_diameter, p_tria->get_communicator());
3185  else
3186 #endif
3187  global_max_diameter = max_diameter;
3188 
3189  return global_max_diameter;
3190  }
3191 
3192 
3193 
3194  namespace internal
3195  {
3196  namespace FixUpDistortedChildCells
3197  {
3198  // compute the mean square
3199  // deviation of the alternating
3200  // forms of the children of the
3201  // given object from that of
3202  // the object itself. for
3203  // objects with
3204  // structdim==spacedim, the
3205  // alternating form is the
3206  // determinant of the jacobian,
3207  // whereas for faces with
3208  // structdim==spacedim-1, the
3209  // alternating form is the
3210  // (signed and scaled) normal
3211  // vector
3212  //
3213  // this average square
3214  // deviation is computed for an
3215  // object where the center node
3216  // has been replaced by the
3217  // second argument to this
3218  // function
3219  template <typename Iterator, int spacedim>
3220  double
3221  objective_function(const Iterator & object,
3222  const Point<spacedim> &object_mid_point)
3223  {
3224  const unsigned int structdim =
3225  Iterator::AccessorType::structure_dimension;
3226  Assert(spacedim == Iterator::AccessorType::dimension,
3227  ExcInternalError());
3228 
3229  // everything below is wrong
3230  // if not for the following
3231  // condition
3232  Assert(object->refinement_case() ==
3234  ExcNotImplemented());
3235  // first calculate the
3236  // average alternating form
3237  // for the parent cell/face
3240  Tensor<spacedim - structdim, spacedim>
3241  parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
3242 
3243  for (const unsigned int i : object->vertex_indices())
3244  parent_vertices[i] = object->vertex(i);
3245 
3247  parent_vertices, parent_alternating_forms);
3248 
3249  const Tensor<spacedim - structdim, spacedim>
3250  average_parent_alternating_form =
3251  std::accumulate(parent_alternating_forms,
3252  parent_alternating_forms +
3255 
3256  // now do the same
3257  // computation for the
3258  // children where we use the
3259  // given location for the
3260  // object mid point instead of
3261  // the one the triangulation
3262  // currently reports
3266  Tensor<spacedim - structdim, spacedim> child_alternating_forms
3269 
3270  for (unsigned int c = 0; c < object->n_children(); ++c)
3271  for (const unsigned int i : object->child(c)->vertex_indices())
3272  child_vertices[c][i] = object->child(c)->vertex(i);
3273 
3274  // replace mid-object
3275  // vertex. note that for
3276  // child i, the mid-object
3277  // vertex happens to have the
3278  // number
3279  // max_children_per_cell-i
3280  for (unsigned int c = 0; c < object->n_children(); ++c)
3281  child_vertices[c][GeometryInfo<structdim>::max_children_per_cell - c -
3282  1] = object_mid_point;
3283 
3284  for (unsigned int c = 0; c < object->n_children(); ++c)
3286  child_vertices[c], child_alternating_forms[c]);
3287 
3288  // on a uniformly refined
3289  // hypercube object, the child
3290  // alternating forms should
3291  // all be smaller by a factor
3292  // of 2^structdim than the
3293  // ones of the parent. as a
3294  // consequence, we'll use the
3295  // squared deviation from
3296  // this ideal value as an
3297  // objective function
3298  double objective = 0;
3299  for (unsigned int c = 0; c < object->n_children(); ++c)
3300  for (const unsigned int i : object->child(c)->vertex_indices())
3301  objective +=
3302  (child_alternating_forms[c][i] -
3303  average_parent_alternating_form / std::pow(2., 1. * structdim))
3304  .norm_square();
3305 
3306  return objective;
3307  }
3308 
3309 
3315  template <typename Iterator>
3317  get_face_midpoint(const Iterator & object,
3318  const unsigned int f,
3319  std::integral_constant<int, 1>)
3320  {
3321  return object->vertex(f);
3322  }
3323 
3324 
3325 
3331  template <typename Iterator>
3333  get_face_midpoint(const Iterator & object,
3334  const unsigned int f,
3335  std::integral_constant<int, 2>)
3336  {
3337  return object->line(f)->center();
3338  }
3339 
3340 
3341 
3347  template <typename Iterator>
3349  get_face_midpoint(const Iterator & object,
3350  const unsigned int f,
3351  std::integral_constant<int, 3>)
3352  {
3353  return object->face(f)->center();
3354  }
3355 
3356 
3357 
3380  template <typename Iterator>
3381  double
3382  minimal_diameter(const Iterator &object)
3383  {
3384  const unsigned int structdim =
3385  Iterator::AccessorType::structure_dimension;
3386 
3387  double diameter = object->diameter();
3388  for (const unsigned int f : object->face_indices())
3389  for (unsigned int e = f + 1; e < object->n_faces(); ++e)
3390  diameter = std::min(
3391  diameter,
3392  get_face_midpoint(object,
3393  f,
3394  std::integral_constant<int, structdim>())
3395  .distance(get_face_midpoint(
3396  object, e, std::integral_constant<int, structdim>())));
3397 
3398  return diameter;
3399  }
3400 
3401 
3402 
3407  template <typename Iterator>
3408  bool
3409  fix_up_object(const Iterator &object)
3410  {
3411  const unsigned int structdim =
3412  Iterator::AccessorType::structure_dimension;
3413  const unsigned int spacedim = Iterator::AccessorType::space_dimension;
3414 
3415  // right now we can only deal with cells that have been refined
3416  // isotropically because that is the only case where we have a cell
3417  // mid-point that can be moved around without having to consider
3418  // boundary information
3419  Assert(object->has_children(), ExcInternalError());
3420  Assert(object->refinement_case() ==
3422  ExcNotImplemented());
3423 
3424  // get the current location of the object mid-vertex:
3425  Point<spacedim> object_mid_point = object->child(0)->vertex(
3427 
3428  // now do a few steepest descent steps to reduce the objective
3429  // function. compute the diameter in the helper function above
3430  unsigned int iteration = 0;
3431  const double diameter = minimal_diameter(object);
3432 
3433  // current value of objective function and initial delta
3434  double current_value = objective_function(object, object_mid_point);
3435  double initial_delta = 0;
3436 
3437  do
3438  {
3439  // choose a step length that is initially 1/4 of the child
3440  // objects' diameter, and a sequence whose sum does not converge
3441  // (to avoid premature termination of the iteration)
3442  const double step_length = diameter / 4 / (iteration + 1);
3443 
3444  // compute the objective function's derivative using a two-sided
3445  // difference formula with eps=step_length/10
3446  Tensor<1, spacedim> gradient;
3447  for (unsigned int d = 0; d < spacedim; ++d)
3448  {
3449  const double eps = step_length / 10;
3450 
3452  h[d] = eps / 2;
3453 
3454  gradient[d] =
3456  object, project_to_object(object, object_mid_point + h)) -
3458  object, project_to_object(object, object_mid_point - h))) /
3459  eps;
3460  }
3461 
3462  // there is nowhere to go
3463  if (gradient.norm() == 0)
3464  break;
3465 
3466  // We need to go in direction -gradient. the optimal value of the
3467  // objective function is zero, so assuming that the model is
3468  // quadratic we would have to go -2*val/||gradient|| in this
3469  // direction, make sure we go at most step_length into this
3470  // direction
3471  object_mid_point -=
3472  std::min(2 * current_value / (gradient * gradient),
3473  step_length / gradient.norm()) *
3474  gradient;
3475  object_mid_point = project_to_object(object, object_mid_point);
3476 
3477  // compute current value of the objective function
3478  const double previous_value = current_value;
3479  current_value = objective_function(object, object_mid_point);
3480 
3481  if (iteration == 0)
3482  initial_delta = (previous_value - current_value);
3483 
3484  // stop if we aren't moving much any more
3485  if ((iteration >= 1) &&
3486  ((previous_value - current_value < 0) ||
3487  (std::fabs(previous_value - current_value) <
3488  0.001 * initial_delta)))
3489  break;
3490 
3491  ++iteration;
3492  }
3493  while (iteration < 20);
3494 
3495  // verify that the new
3496  // location is indeed better
3497  // than the one before. check
3498  // this by comparing whether
3499  // the minimum value of the
3500  // products of parent and
3501  // child alternating forms is
3502  // positive. for cells this
3503  // means that the
3504  // determinants have the same
3505  // sign, for faces that the
3506  // face normals of parent and
3507  // children point in the same
3508  // general direction
3509  double old_min_product, new_min_product;
3510 
3513  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3514  parent_vertices[i] = object->vertex(i);
3515 
3516  Tensor<spacedim - structdim, spacedim>
3517  parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
3519  parent_vertices, parent_alternating_forms);
3520 
3524 
3525  for (unsigned int c = 0; c < object->n_children(); ++c)
3526  for (const unsigned int i : object->child(c)->vertex_indices())
3527  child_vertices[c][i] = object->child(c)->vertex(i);
3528 
3529  Tensor<spacedim - structdim, spacedim> child_alternating_forms
3532 
3533  for (unsigned int c = 0; c < object->n_children(); ++c)
3535  child_vertices[c], child_alternating_forms[c]);
3536 
3537  old_min_product =
3538  child_alternating_forms[0][0] * parent_alternating_forms[0];
3539  for (unsigned int c = 0; c < object->n_children(); ++c)
3540  for (const unsigned int i : object->child(c)->vertex_indices())
3541  for (const unsigned int j : object->vertex_indices())
3542  old_min_product = std::min<double>(old_min_product,
3543  child_alternating_forms[c][i] *
3544  parent_alternating_forms[j]);
3545 
3546  // for the new minimum value,
3547  // replace mid-object
3548  // vertex. note that for child
3549  // i, the mid-object vertex
3550  // happens to have the number
3551  // max_children_per_cell-i
3552  for (unsigned int c = 0; c < object->n_children(); ++c)
3553  child_vertices[c][GeometryInfo<structdim>::max_children_per_cell - c -
3554  1] = object_mid_point;
3555 
3556  for (unsigned int c = 0; c < object->n_children(); ++c)
3558  child_vertices[c], child_alternating_forms[c]);
3559 
3560  new_min_product =
3561  child_alternating_forms[0][0] * parent_alternating_forms[0];
3562  for (unsigned int c = 0; c < object->n_children(); ++c)
3563  for (const unsigned int i : object->child(c)->vertex_indices())
3564  for (const unsigned int j : object->vertex_indices())
3565  new_min_product = std::min<double>(new_min_product,
3566  child_alternating_forms[c][i] *
3567  parent_alternating_forms[j]);
3568 
3569  // if new minimum value is
3570  // better than before, then set the
3571  // new mid point. otherwise
3572  // return this object as one of
3573  // those that can't apparently
3574  // be fixed
3575  if (new_min_product >= old_min_product)
3576  object->child(0)->vertex(
3578  object_mid_point;
3579 
3580  // return whether after this
3581  // operation we have an object that
3582  // is well oriented
3583  return (std::max(new_min_product, old_min_product) > 0);
3584  }
3585 
3586 
3587 
3588  // possibly fix up the faces of a cell by moving around its mid-points
3589  template <int dim, int spacedim>
3590  void
3592  const typename ::Triangulation<dim, spacedim>::cell_iterator
3593  &cell,
3594  std::integral_constant<int, dim>,
3595  std::integral_constant<int, spacedim>)
3596  {
3597  // see if we first can fix up some of the faces of this object. We can
3598  // mess with faces if and only if the neighboring cell is not even
3599  // more refined than we are (since in that case the sub-faces have
3600  // themselves children that we can't move around any more). however,
3601  // the latter case shouldn't happen anyway: if the current face is
3602  // distorted but the neighbor is even more refined, then the face had
3603  // been deformed before already, and had been ignored at the time; we
3604  // should then also be able to ignore it this time as well
3605  for (auto f : cell->face_indices())
3606  {
3607  Assert(cell->face(f)->has_children(), ExcInternalError());
3608  Assert(cell->face(f)->refinement_case() ==
3609  RefinementCase<dim - 1>::isotropic_refinement,
3610  ExcInternalError());
3611 
3612  bool subface_is_more_refined = false;
3613  for (unsigned int g = 0;
3614  g < GeometryInfo<dim>::max_children_per_face;
3615  ++g)
3616  if (cell->face(f)->child(g)->has_children())
3617  {
3618  subface_is_more_refined = true;
3619  break;
3620  }
3621 
3622  if (subface_is_more_refined == true)
3623  continue;
3624 
3625  // we finally know that we can do something about this face
3626  fix_up_object(cell->face(f));
3627  }
3628  }
3629  } /* namespace FixUpDistortedChildCells */
3630  } /* namespace internal */
3631 
3632 
3633  template <int dim, int spacedim>
3637  &distorted_cells,
3638  Triangulation<dim, spacedim> & /*triangulation*/)
3639  {
3640  static_assert(
3641  dim != 1 && spacedim != 1,
3642  "This function is only valid when dim != 1 or spacedim != 1.");
3643  typename Triangulation<dim, spacedim>::DistortedCellList unfixable_subset;
3644 
3645  // loop over all cells that we have to fix up
3646  for (typename std::list<
3647  typename Triangulation<dim, spacedim>::cell_iterator>::const_iterator
3648  cell_ptr = distorted_cells.distorted_cells.begin();
3649  cell_ptr != distorted_cells.distorted_cells.end();
3650  ++cell_ptr)
3651  {
3652  const typename Triangulation<dim, spacedim>::cell_iterator cell =
3653  *cell_ptr;
3654 
3655  Assert(!cell->is_active(),
3656  ExcMessage(
3657  "This function is only valid for a list of cells that "
3658  "have children (i.e., no cell in the list may be active)."));
3659 
3661  cell,
3662  std::integral_constant<int, dim>(),
3663  std::integral_constant<int, spacedim>());
3664 
3665  // If possible, fix up the object.
3667  unfixable_subset.distorted_cells.push_back(cell);
3668  }
3669 
3670  return unfixable_subset;
3671  }
3672 
3673 
3674 
3675  template <int dim, int spacedim>
3676  void
3678  const bool reset_boundary_ids)
3679  {
3680  const auto src_boundary_ids = tria.get_boundary_ids();
3681  std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
3682  auto m_it = dst_manifold_ids.begin();
3683  for (const auto b : src_boundary_ids)
3684  {
3685  *m_it = static_cast<types::manifold_id>(b);
3686  ++m_it;
3687  }
3688  const std::vector<types::boundary_id> reset_boundary_id =
3689  reset_boundary_ids ?
3690  std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
3691  src_boundary_ids;
3692  map_boundary_to_manifold_ids(src_boundary_ids,
3693  dst_manifold_ids,
3694  tria,
3695  reset_boundary_id);
3696  }
3697 
3698 
3699 
3700  template <int dim, int spacedim>
3701  void
3703  const std::vector<types::boundary_id> &src_boundary_ids,
3704  const std::vector<types::manifold_id> &dst_manifold_ids,
3706  const std::vector<types::boundary_id> &reset_boundary_ids_)
3707  {
3708  AssertDimension(src_boundary_ids.size(), dst_manifold_ids.size());
3709  const auto reset_boundary_ids =
3710  reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
3711  AssertDimension(reset_boundary_ids.size(), src_boundary_ids.size());
3712 
3713  // in 3d, we not only have to copy boundary ids of faces, but also of edges
3714  // because we see them twice (once from each adjacent boundary face),
3715  // we cannot immediately reset their boundary ids. thus, copy first
3716  // and reset later
3717  if (dim >= 3)
3718  for (const auto &cell : tria.active_cell_iterators())
3719  for (auto f : cell->face_indices())
3720  if (cell->face(f)->at_boundary())
3721  for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
3722  {
3723  const auto bid = cell->face(f)->line(e)->boundary_id();
3724  const unsigned int ind = std::find(src_boundary_ids.begin(),
3725  src_boundary_ids.end(),
3726  bid) -
3727  src_boundary_ids.begin();
3728  if (ind < src_boundary_ids.size())
3729  cell->face(f)->line(e)->set_manifold_id(
3730  dst_manifold_ids[ind]);
3731  }
3732 
3733  // now do cells
3734  for (const auto &cell : tria.active_cell_iterators())
3735  for (auto f : cell->face_indices())
3736  if (cell->face(f)->at_boundary())
3737  {
3738  const auto bid = cell->face(f)->boundary_id();
3739  const unsigned int ind =
3740  std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
3741  src_boundary_ids.begin();
3742 
3743  if (ind < src_boundary_ids.size())
3744  {
3745  // assign the manifold id
3746  cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
3747  // then reset boundary id
3748  cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
3749  }
3750 
3751  if (dim >= 3)
3752  for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
3753  {
3754  const auto bid = cell->face(f)->line(e)->boundary_id();
3755  const unsigned int ind = std::find(src_boundary_ids.begin(),
3756  src_boundary_ids.end(),
3757  bid) -
3758  src_boundary_ids.begin();
3759  if (ind < src_boundary_ids.size())
3760  cell->face(f)->line(e)->set_boundary_id(
3761  reset_boundary_ids[ind]);
3762  }
3763  }
3764  }
3765 
3766 
3767  template <int dim, int spacedim>
3768  void
3770  const bool compute_face_ids)
3771  {
3773  cell = tria.begin_active(),
3774  endc = tria.end();
3775 
3776  for (; cell != endc; ++cell)
3777  {
3778  cell->set_manifold_id(cell->material_id());
3779  if (compute_face_ids == true)
3780  {
3781  for (auto f : cell->face_indices())
3782  {
3783  if (cell->at_boundary(f) == false)
3784  cell->face(f)->set_manifold_id(
3785  std::min(cell->material_id(),
3786  cell->neighbor(f)->material_id()));
3787  else
3788  cell->face(f)->set_manifold_id(cell->material_id());
3789  }
3790  }
3791  }
3792  }
3793 
3794 
3795  template <int dim, int spacedim>
3796  void
3799  const std::function<types::manifold_id(
3800  const std::set<types::manifold_id> &)> &disambiguation_function,
3801  bool overwrite_only_flat_manifold_ids)
3802  {
3803  // Easy case first:
3804  if (dim == 1)
3805  return;
3806  const unsigned int n_subobjects =
3807  dim == 2 ? tria.n_lines() : tria.n_lines() + tria.n_quads();
3808 
3809  // If user index is zero, then it has not been set.
3810  std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
3811  std::vector<unsigned int> backup;
3812  tria.save_user_indices(backup);
3813  tria.clear_user_data();
3814 
3815  unsigned next_index = 1;
3816  for (auto &cell : tria.active_cell_iterators())
3817  {
3818  if (dim > 1)
3819  for (unsigned int l = 0; l < cell->n_lines(); ++l)
3820  {
3821  if (cell->line(l)->user_index() == 0)
3822  {
3823  AssertIndexRange(next_index, n_subobjects + 1);
3824  manifold_ids[next_index].insert(cell->manifold_id());
3825  cell->line(l)->set_user_index(next_index++);
3826  }
3827  else
3828  manifold_ids[cell->line(l)->user_index()].insert(
3829  cell->manifold_id());
3830  }
3831  if (dim > 2)
3832  for (unsigned int l = 0; l < cell->n_faces(); ++l)
3833  {
3834  if (cell->quad(l)->user_index() == 0)
3835  {
3836  AssertIndexRange(next_index, n_subobjects + 1);
3837  manifold_ids[next_index].insert(cell->manifold_id());
3838  cell->quad(l)->set_user_index(next_index++);
3839  }
3840  else
3841  manifold_ids[cell->quad(l)->user_index()].insert(
3842  cell->manifold_id());
3843  }
3844  }
3845  for (auto &cell : tria.active_cell_iterators())
3846  {
3847  if (dim > 1)
3848  for (unsigned int l = 0; l < cell->n_lines(); ++l)
3849  {
3850  const auto id = cell->line(l)->user_index();
3851  // Make sure we change the manifold indicator only once
3852  if (id != 0)
3853  {
3854  if (cell->line(l)->manifold_id() ==
3856  overwrite_only_flat_manifold_ids == false)
3857  cell->line(l)->set_manifold_id(
3858  disambiguation_function(manifold_ids[id]));
3859  cell->line(l)->set_user_index(0);
3860  }
3861  }
3862  if (dim > 2)
3863  for (unsigned int l = 0; l < cell->n_faces(); ++l)
3864  {
3865  const auto id = cell->quad(l)->user_index();
3866  // Make sure we change the manifold indicator only once
3867  if (id != 0)
3868  {
3869  if (cell->quad(l)->manifold_id() ==
3871  overwrite_only_flat_manifold_ids == false)
3872  cell->quad(l)->set_manifold_id(
3873  disambiguation_function(manifold_ids[id]));
3874  cell->quad(l)->set_user_index(0);
3875  }
3876  }
3877  }
3878  tria.load_user_indices(backup);
3879  }
3880 
3881 
3882 
3883  template <int dim, int spacedim>
3884  std::pair<unsigned int, double>
3887  {
3888  double max_ratio = 1;
3889  unsigned int index = 0;
3890 
3891  for (unsigned int i = 0; i < dim; ++i)
3892  for (unsigned int j = i + 1; j < dim; ++j)
3893  {
3894  unsigned int ax = i % dim;
3895  unsigned int next_ax = j % dim;
3896 
3897  double ratio =
3898  cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
3899 
3900  if (ratio > max_ratio)
3901  {
3902  max_ratio = ratio;
3903  index = ax;
3904  }
3905  else if (1.0 / ratio > max_ratio)
3906  {
3907  max_ratio = 1.0 / ratio;
3908  index = next_ax;
3909  }
3910  }
3911  return std::make_pair(index, max_ratio);
3912  }
3913 
3914 
3915  template <int dim, int spacedim>
3916  void
3918  const bool isotropic,
3919  const unsigned int max_iterations)
3920  {
3921  unsigned int iter = 0;
3922  bool continue_refinement = true;
3923 
3924  while (continue_refinement && (iter < max_iterations))
3925  {
3926  if (max_iterations != numbers::invalid_unsigned_int)
3927  iter++;
3928  continue_refinement = false;
3929 
3930  for (const auto &cell : tria.active_cell_iterators())
3931  for (const unsigned int j : cell->face_indices())
3932  if (cell->at_boundary(j) == false &&
3933  cell->neighbor(j)->has_children())
3934  {
3935  if (isotropic)
3936  {
3937  cell->set_refine_flag();
3938  continue_refinement = true;
3939  }
3940  else
3941  continue_refinement |= cell->flag_for_face_refinement(j);
3942  }
3943 
3945  }
3946  }
3947 
3948  template <int dim, int spacedim>
3949  void
3951  const double max_ratio,
3952  const unsigned int max_iterations)
3953  {
3954  unsigned int iter = 0;
3955  bool continue_refinement = true;
3956 
3957  while (continue_refinement && (iter < max_iterations))
3958  {
3959  iter++;
3960  continue_refinement = false;
3961  for (const auto &cell : tria.active_cell_iterators())
3962  {
3963  std::pair<unsigned int, double> info =
3964  GridTools::get_longest_direction<dim, spacedim>(cell);
3965  if (info.second > max_ratio)
3966  {
3967  cell->set_refine_flag(
3968  RefinementCase<dim>::cut_axis(info.first));
3969  continue_refinement = true;
3970  }
3971  }
3973  }
3974  }
3975 
3976 
3977  template <int dim, int spacedim>
3978  void
3980  const double limit_angle_fraction)
3981  {
3982  if (dim == 1)
3983  return; // Nothing to do
3984 
3985  // Check that we don't have hanging nodes
3987  ExcMessage("The input Triangulation cannot "
3988  "have hanging nodes."));
3989 
3990 
3991  bool has_cells_with_more_than_dim_faces_on_boundary = true;
3992  bool has_cells_with_dim_faces_on_boundary = false;
3993 
3994  unsigned int refinement_cycles = 0;
3995 
3996  while (has_cells_with_more_than_dim_faces_on_boundary)
3997  {
3998  has_cells_with_more_than_dim_faces_on_boundary = false;
3999 
4000  for (const auto &cell : tria.active_cell_iterators())
4001  {
4002  unsigned int boundary_face_counter = 0;
4003  for (auto f : cell->face_indices())
4004  if (cell->face(f)->at_boundary())
4005  boundary_face_counter++;
4006  if (boundary_face_counter > dim)
4007  {
4008  has_cells_with_more_than_dim_faces_on_boundary = true;
4009  break;
4010  }
4011  else if (boundary_face_counter == dim)
4012  has_cells_with_dim_faces_on_boundary = true;
4013  }
4014  if (has_cells_with_more_than_dim_faces_on_boundary)
4015  {
4016  tria.refine_global(1);
4017  refinement_cycles++;
4018  }
4019  }
4020 
4021  if (has_cells_with_dim_faces_on_boundary)
4022  {
4023  tria.refine_global(1);
4024  refinement_cycles++;
4025  }
4026  else
4027  {
4028  while (refinement_cycles > 0)
4029  {
4030  for (const auto &cell : tria.active_cell_iterators())
4031  cell->set_coarsen_flag();
4033  refinement_cycles--;
4034  }
4035  return;
4036  }
4037 
4038  std::vector<bool> cells_to_remove(tria.n_active_cells(), false);
4039  std::vector<Point<spacedim>> vertices = tria.get_vertices();
4040 
4041  std::vector<bool> faces_to_remove(tria.n_raw_faces(), false);
4042 
4043  std::vector<CellData<dim>> cells_to_add;
4044  SubCellData subcelldata_to_add;
4045 
4046  // Trick compiler for dimension independent things
4047  const unsigned int v0 = 0, v1 = 1, v2 = (dim > 1 ? 2 : 0),
4048  v3 = (dim > 1 ? 3 : 0);
4049 
4050  for (const auto &cell : tria.active_cell_iterators())
4051  {
4052  double angle_fraction = 0;
4053  unsigned int vertex_at_corner = numbers::invalid_unsigned_int;
4054 
4055  if (dim == 2)
4056  {
4058  p0[spacedim > 1 ? 1 : 0] = 1;
4060  p1[0] = 1;
4061 
4062  if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
4063  {
4064  p0 = cell->vertex(v0) - cell->vertex(v2);
4065  p1 = cell->vertex(v3) - cell->vertex(v2);
4066  vertex_at_corner = v2;
4067  }
4068  else if (cell->face(v3)->at_boundary() &&
4069  cell->face(v1)->at_boundary())
4070  {
4071  p0 = cell->vertex(v2) - cell->vertex(v3);
4072  p1 = cell->vertex(v1) - cell->vertex(v3);
4073  vertex_at_corner = v3;
4074  }
4075  else if (cell->face(1)->at_boundary() &&
4076  cell->face(2)->at_boundary())
4077  {
4078  p0 = cell->vertex(v0) - cell->vertex(v1);
4079  p1 = cell->vertex(v3) - cell->vertex(v1);
4080  vertex_at_corner = v1;
4081  }
4082  else if (cell->face(2)->at_boundary() &&
4083  cell->face(0)->at_boundary())
4084  {
4085  p0 = cell->vertex(v2) - cell->vertex(v0);
4086  p1 = cell->vertex(v1) - cell->vertex(v0);
4087  vertex_at_corner = v0;
4088  }
4089  p0 /= p0.norm();
4090  p1 /= p1.norm();
4091  angle_fraction = std::acos(p0 * p1) / numbers::PI;
4092  }
4093  else
4094  {
4095  Assert(false, ExcNotImplemented());
4096  }
4097 
4098  if (angle_fraction > limit_angle_fraction)
4099  {
4100  auto flags_removal = [&](unsigned int f1,
4101  unsigned int f2,
4102  unsigned int n1,
4103  unsigned int n2) -> void {
4104  cells_to_remove[cell->active_cell_index()] = true;
4105  cells_to_remove[cell->neighbor(n1)->active_cell_index()] = true;
4106  cells_to_remove[cell->neighbor(n2)->active_cell_index()] = true;
4107 
4108  faces_to_remove[cell->face(f1)->index()] = true;
4109  faces_to_remove[cell->face(f2)->index()] = true;
4110 
4111  faces_to_remove[cell->neighbor(n1)->face(f1)->index()] = true;
4112  faces_to_remove[cell->neighbor(n2)->face(f2)->index()] = true;
4113  };
4114 
4115  auto cell_creation = [&](const unsigned int vv0,
4116  const unsigned int vv1,
4117  const unsigned int f0,
4118  const unsigned int f1,
4119 
4120  const unsigned int n0,
4121  const unsigned int v0n0,
4122  const unsigned int v1n0,
4123 
4124  const unsigned int n1,
4125  const unsigned int v0n1,
4126  const unsigned int v1n1) {
4127  CellData<dim> c1, c2;
4128  CellData<1> l1, l2;
4129 
4130  c1.vertices[v0] = cell->vertex_index(vv0);
4131  c1.vertices[v1] = cell->vertex_index(vv1);
4132  c1.vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
4133  c1.vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
4134 
4135  c1.manifold_id = cell->manifold_id();
4136  c1.material_id = cell->material_id();
4137 
4138  c2.vertices[v0] = cell->vertex_index(vv0);
4139  c2.vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
4140  c2.vertices[v2] = cell->vertex_index(vv1);
4141  c2.vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
4142 
4143  c2.manifold_id = cell->manifold_id();
4144  c2.material_id = cell->material_id();
4145 
4146  l1.vertices[0] = cell->vertex_index(vv0);
4147  l1.vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
4148 
4149  l1.boundary_id = cell->line(f0)->boundary_id();
4150  l1.manifold_id = cell->line(f0)->manifold_id();
4151  subcelldata_to_add.boundary_lines.push_back(l1);
4152 
4153  l2.vertices[0] = cell->vertex_index(vv0);
4154  l2.vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
4155 
4156  l2.boundary_id = cell->line(f1)->boundary_id();
4157  l2.manifold_id = cell->line(f1)->manifold_id();
4158  subcelldata_to_add.boundary_lines.push_back(l2);
4159 
4160  cells_to_add.push_back(c1);
4161  cells_to_add.push_back(c2);
4162  };
4163 
4164  if (dim == 2)
4165  {
4166  switch (vertex_at_corner)
4167  {
4168  case 0:
4169  flags_removal(0, 2, 3, 1);
4170  cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
4171  break;
4172  case 1:
4173  flags_removal(1, 2, 3, 0);
4174  cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
4175  break;
4176  case 2:
4177  flags_removal(3, 0, 1, 2);
4178  cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
4179  break;
4180  case 3:
4181  flags_removal(3, 1, 0, 2);
4182  cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
4183  break;
4184  }
4185  }
4186  else
4187  {
4188  Assert(false, ExcNotImplemented());
4189  }
4190  }
4191  }
4192 
4193  // if no cells need to be added, then no regularization is necessary.
4194  // Restore things as they were before this function was called.
4195  if (cells_to_add.size() == 0)
4196  {
4197  while (refinement_cycles > 0)
4198  {
4199  for (const auto &cell : tria.active_cell_iterators())
4200  cell->set_coarsen_flag();
4202  refinement_cycles--;
4203  }
4204  return;
4205  }
4206 
4207  // add the cells that were not marked as skipped
4208  for (const auto &cell : tria.active_cell_iterators())
4209  {
4210  if (cells_to_remove[cell->active_cell_index()] == false)
4211  {
4212  CellData<dim> c;
4213  for (const unsigned int v : cell->vertex_indices())
4214  c.vertices[v] = cell->vertex_index(v);
4215  c.manifold_id = cell->manifold_id();
4216  c.material_id = cell->material_id();
4217  cells_to_add.push_back(c);
4218  }
4219  }
4220 
4221  // Face counter for both dim == 2 and dim == 3
4223  face = tria.begin_active_face(),
4224  endf = tria.end_face();
4225  for (; face != endf; ++face)
4226  if ((face->at_boundary() ||
4227  face->manifold_id() != numbers::flat_manifold_id) &&
4228  faces_to_remove[face->index()] == false)
4229  {
4230  for (unsigned int l = 0; l < face->n_lines(); ++l)
4231  {
4232  CellData<1> line;
4233  if (dim == 2)
4234  {
4235  for (const unsigned int v : face->vertex_indices())
4236  line.vertices[v] = face->vertex_index(v);
4237  line.boundary_id = face->boundary_id();
4238  line.manifold_id = face->manifold_id();
4239  }
4240  else
4241  {
4242  for (const unsigned int v : face->line(l)->vertex_indices())
4243  line.vertices[v] = face->line(l)->vertex_index(v);
4244  line.boundary_id = face->line(l)->boundary_id();
4245  line.manifold_id = face->line(l)->manifold_id();
4246  }
4247  subcelldata_to_add.boundary_lines.push_back(line);
4248  }
4249  if (dim == 3)
4250  {
4251  CellData<2> quad;
4252  for (const unsigned int v : face->vertex_indices())
4253  quad.vertices[v] = face->vertex_index(v);
4254  quad.boundary_id = face->boundary_id();
4255  quad.manifold_id = face->manifold_id();
4256  subcelldata_to_add.boundary_quads.push_back(quad);
4257  }
4258  }
4260  cells_to_add,
4261  subcelldata_to_add);
4263 
4264  // Save manifolds
4265  auto manifold_ids = tria.get_manifold_ids();
4266  std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
4267  manifolds;
4268  // Set manifolds in new Triangulation
4269  for (const auto manifold_id : manifold_ids)
4271  manifolds[manifold_id] = tria.get_manifold(manifold_id).clone();
4272 
4273  tria.clear();
4274 
4275  tria.create_triangulation(vertices, cells_to_add, subcelldata_to_add);
4276 
4277  // Restore manifolds
4278  for (const auto manifold_id : manifold_ids)
4280  tria.set_manifold(manifold_id, *manifolds[manifold_id]);
4281  }
4282 
4283 
4284 
4285  template <int dim, int spacedim>
4286 #ifndef DOXYGEN
4287  std::tuple<
4288  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
4289  std::vector<std::vector<Point<dim>>>,
4290  std::vector<std::vector<unsigned int>>>
4291 #else
4292  return_type
4293 #endif
4295  const Cache<dim, spacedim> & cache,
4296  const std::vector<Point<spacedim>> &points,
4298  &cell_hint)
4299  {
4300  const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
4301  // Splitting the tuple's components
4302  auto &cells = std::get<0>(cqmp);
4303  auto &qpoints = std::get<1>(cqmp);
4304  auto &maps = std::get<2>(cqmp);
4305  auto &missing_points = std::get<3>(cqmp);
4306  // If a point was not found, throwing an error, as the old
4307  // implementation of compute_point_locations would have done
4308  AssertThrow(std::get<3>(cqmp).size() == 0,
4309  ExcPointNotFound<spacedim>(points[missing_points[0]]));
4310 
4311  (void)missing_points;
4312 
4313  return std::make_tuple(std::move(cells),
4314  std::move(qpoints),
4315  std::move(maps));
4316  }
4317 
4318 
4319 
4320  template <int dim, int spacedim>
4321 #ifndef DOXYGEN
4322  std::tuple<
4323  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
4324  std::vector<std::vector<Point<dim>>>,
4325  std::vector<std::vector<unsigned int>>,
4326  std::vector<unsigned int>>
4327 #else
4328  return_type
4329 #endif
4331  const Cache<dim, spacedim> & cache,
4332  const std::vector<Point<spacedim>> &points,
4334  &cell_hint)
4335  {
4336  // How many points are here?
4337  const unsigned int np = points.size();
4338 
4339  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
4340  cells_out;
4341  std::vector<std::vector<Point<dim>>> qpoints_out;
4342  std::vector<std::vector<unsigned int>> maps_out;
4343  std::vector<unsigned int> missing_points_out;
4344 
4345  // Now the easy case.
4346  if (np == 0)
4347  return std::make_tuple(std::move(cells_out),
4348  std::move(qpoints_out),
4349  std::move(maps_out),
4350  std::move(missing_points_out));
4351 
4352  // For the search we shall use the following tree
4353  const auto &b_tree = cache.get_cell_bounding_boxes_rtree();
4354 
4355  // We begin by finding the cell/transform of the first point
4356  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4357  Point<dim>>
4358  my_pair;
4359 
4360  bool found = false;
4361  unsigned int points_checked = 0;
4362 
4363  // If a hint cell was given, use it
4364  if (cell_hint.state() == IteratorState::valid)
4365  {
4366  try
4367  {
4369  points[0],
4370  cell_hint);
4371  found = true;
4372  }
4373  catch (const GridTools::ExcPointNotFound<dim> &)
4374  {
4375  missing_points_out.emplace_back(0);
4376  }
4377  ++points_checked;
4378  }
4379 
4380  // The tree search returns
4381  // - a bounding box covering the cell
4382  // - the active cell iterator
4383  std::vector<
4384  std::pair<BoundingBox<spacedim>,
4386  box_cell;
4387 
4388  // This is used as an index for box_cell
4389  int cell_candidate_idx = -1;
4390  // If any of the cells in box_cell is a ghost cell,
4391  // an artificial cell or at the boundary,
4392  // we want to use try/catch
4393  bool use_try = false;
4394 
4395  while (!found && points_checked < np)
4396  {
4397  box_cell.clear();
4398  b_tree.query(boost::geometry::index::intersects(points[points_checked]),
4399  std::back_inserter(box_cell));
4400 
4401  // Checking box_cell result for a suitable candidate
4402  cell_candidate_idx = -1;
4403  for (unsigned int i = 0; i < box_cell.size(); ++i)
4404  {
4405  // As a candidate we don't want artificial cells
4406  if (!box_cell[i].second->is_artificial())
4407  cell_candidate_idx = i;
4408 
4409  // If the cell is not locally owned or at boundary
4410  // we check for exceptions
4411  if (cell_candidate_idx != -1 &&
4412  (!box_cell[i].second->is_locally_owned() ||
4413  box_cell[i].second->at_boundary()))
4414  use_try = true;
4415 
4416 
4417  if (cell_candidate_idx != -1)
4418  break;
4419  }
4420 
4421  // If a suitable cell was found, use it as hint
4422  if (cell_candidate_idx != -1)
4423  {
4424  if (use_try)
4425  {
4426  try
4427  {
4429  cache,
4430  points[points_checked],
4431  box_cell[cell_candidate_idx].second);
4432  found = true;
4433  }
4434  catch (const GridTools::ExcPointNotFound<dim> &)
4435  {
4436  missing_points_out.emplace_back(points_checked);
4437  }
4438  }
4439  else
4440  {
4442  cache,
4443  points[points_checked],
4444  box_cell[cell_candidate_idx].second);
4445  found = true;
4446  }
4447  }
4448  else
4449  {
4450  try
4451  {
4453  cache, points[points_checked]);
4454  // If we arrive here the cell was not among
4455  // the candidates returned by the tree, so we're adding it
4456  // by hand
4457  found = true;
4458  cell_candidate_idx = box_cell.size();
4459  box_cell.push_back(
4460  std::make_pair(my_pair.first->bounding_box(), my_pair.first));
4461  }
4462  catch (const GridTools::ExcPointNotFound<dim> &)
4463  {
4464  missing_points_out.emplace_back(points_checked);
4465  }
4466  }
4467 
4468  // Updating the position of the analyzed points
4469  ++points_checked;
4470  }
4471 
4472  // If the point has been found in a cell, adding it
4473  if (found)
4474  {
4475  cells_out.emplace_back(my_pair.first);
4476  qpoints_out.emplace_back(1, my_pair.second);
4477  maps_out.emplace_back(1, points_checked - 1);
4478  }
4479 
4480  // Now the second easy case.
4481  if (np == qpoints_out.size())
4482  return std::make_tuple(std::move(cells_out),
4483  std::move(qpoints_out),
4484  std::move(maps_out),
4485  std::move(missing_points_out));
4486 
4487  // Cycle over all points left
4488  for (unsigned int p = points_checked; p < np; ++p)
4489  {
4490  // We assume the last used cell contains the point: checking it
4491  if (cell_candidate_idx != -1)
4492  if (!box_cell[cell_candidate_idx].first.point_inside(points[p]))
4493  // Point outside candidate cell: we have no candidate
4494  cell_candidate_idx = -1;
4495 
4496  // If there's no candidate, run a tree search
4497  if (cell_candidate_idx == -1)
4498  {
4499  // Using the b_tree to find new candidates
4500  box_cell.clear();
4501  b_tree.query(boost::geometry::index::intersects(points[p]),
4502  std::back_inserter(box_cell));
4503  // Checking the returned bounding boxes/cells
4504  use_try = false;
4505  cell_candidate_idx = -1;
4506  for (unsigned int i = 0; i < box_cell.size(); ++i)
4507  {
4508  // As a candidate we don't want artificial cells
4509  if (!box_cell[i].second->is_artificial())
4510  cell_candidate_idx = i;
4511 
4512  // If the cell is not locally owned or at boundary
4513  // we check for exceptions
4514  if (cell_candidate_idx != -1 &&
4515  (!box_cell[i].second->is_locally_owned() ||
4516  box_cell[i].second->at_boundary()))
4517  use_try = true;
4518 
4519  // If a cell candidate was found we can stop
4520  if (cell_candidate_idx != -1)
4521  break;
4522  }
4523  }
4524 
4525  if (cell_candidate_idx == -1)
4526  {
4527  // No candidate cell, but the cell might
4528  // still be inside the mesh, this is our final check:
4529  try
4530  {
4531  my_pair =
4532  GridTools::find_active_cell_around_point(cache, points[p]);
4533  // If we arrive here the cell was not among
4534  // the candidates returned by the tree, so we're adding it
4535  // by hand
4536  cell_candidate_idx = box_cell.size();
4537  box_cell.push_back(
4538  std::make_pair(my_pair.first->bounding_box(), my_pair.first));
4539  }
4540  catch (const GridTools::ExcPointNotFound<dim> &)
4541  {
4542  missing_points_out.emplace_back(p);
4543  continue;
4544  }
4545  }
4546  else
4547  {
4548  // We have a candidate cell
4549  if (use_try)
4550  {
4551  try
4552  {
4554  cache, points[p], box_cell[cell_candidate_idx].second);
4555  }
4556  catch (const GridTools::ExcPointNotFound<dim> &)
4557  {
4558  missing_points_out.push_back(p);
4559  continue;
4560  }
4561  }
4562  else
4563  {
4565  cache, points[p], box_cell[cell_candidate_idx].second);
4566  }
4567 
4568  // If the point was found in another cell,
4569  // updating cell_candidate_idx
4570  if (my_pair.first != box_cell[cell_candidate_idx].second)
4571  {
4572  for (unsigned int i = 0; i < box_cell.size(); ++i)
4573  {
4574  if (my_pair.first == box_cell[i].second)
4575  {
4576  cell_candidate_idx = i;
4577  break;
4578  }
4579  }
4580 
4581  if (my_pair.first != box_cell[cell_candidate_idx].second)
4582  {
4583  // The cell was not among the candidates returned by the
4584  // tree
4585  cell_candidate_idx = box_cell.size();
4586  box_cell.push_back(
4587  std::make_pair(my_pair.first->bounding_box(),
4588  my_pair.first));
4589  }
4590  }
4591  }
4592 
4593 
4594  // Assuming the point is more likely to be in the last
4595  // used cell
4596  if (my_pair.first == cells_out.back())
4597  {
4598  // Found in the last cell: adding the data
4599  qpoints_out.back().emplace_back(my_pair.second);
4600  maps_out.back().emplace_back(p);
4601  }
4602  else
4603  {
4604  // Check if it is in another cell already found
4605  typename std::vector<typename Triangulation<dim, spacedim>::
4606  active_cell_iterator>::iterator cells_it =
4607  std::find(cells_out.begin(), cells_out.end() - 1, my_pair.first);
4608 
4609  if (cells_it == cells_out.end() - 1)
4610  {
4611  // Cell not found: adding a new cell
4612  cells_out.emplace_back(my_pair.first);
4613  qpoints_out.emplace_back(1, my_pair.second);
4614  maps_out.emplace_back(1, p);
4615  }
4616  else
4617  {
4618  // Cell found: just adding the point index and qpoint to the
4619  // list
4620  unsigned int current_cell = cells_it - cells_out.begin();
4621  qpoints_out[current_cell].emplace_back(my_pair.second);
4622  maps_out[current_cell].emplace_back(p);
4623  }
4624  }
4625  }
4626 
4627  // Debug Checking
4628  Assert(cells_out.size() == maps_out.size(),
4629  ExcDimensionMismatch(cells_out.size(), maps_out.size()));
4630 
4631  Assert(cells_out.size() == qpoints_out.size(),
4632  ExcDimensionMismatch(cells_out.size(), qpoints_out.size()));
4633 
4634 #ifdef DEBUG
4635  unsigned int c = cells_out.size();
4636  unsigned int qps = 0;
4637  // The number of points in all
4638  // the cells must be the same as
4639  // the number of points we
4640  // started off from,
4641  // plus the points which were ignored
4642  for (unsigned int n = 0; n < c; ++n)
4643  {
4644  Assert(qpoints_out[n].size() == maps_out[n].size(),
4645  ExcDimensionMismatch(qpoints_out[n].size(), maps_out[n].size()));
4646  qps += qpoints_out[n].size();
4647  }
4648 
4649  Assert(qps + missing_points_out.size() == np,
4650  ExcDimensionMismatch(qps + missing_points_out.size(), np));
4651 #endif
4652 
4653  return std::make_tuple(std::move(cells_out),
4654  std::move(qpoints_out),
4655  std::move(maps_out),
4656  std::move(missing_points_out));
4657  }
4658 
4659 
4660 
4661  namespace internal
4662  {
4663  // Functions used for distributed compute point locations
4664  namespace DistributedComputePointLocations
4665  {
4666  // Hash function for cells; needed for unordered maps/multimaps
4667  template <int dim, int spacedim>
4668  struct cell_hash
4669  {
4670  std::size_t
4673  const
4674  {
4675  // Return active cell index, which is faster than CellId to compute
4676  return k->active_cell_index();
4677  }
4678  };
4679 
4680 
4681 
4682  // Compute point locations; internal version which returns an unordered
4683  // map. The algorithm is the same as for
4684  // GridTools::compute_point_locations.
4685  template <int dim, int spacedim>
4686  std::unordered_map<
4688  std::pair<std::vector<Point<dim>>, std::vector<unsigned int>>,
4691  const std::vector<Point<spacedim>> & points)
4692  {
4693  const unsigned int n_points = points.size();
4694  // Creating the output tuple
4695  std::unordered_map<
4696  typename Triangulation<dim, spacedim>::active_cell_iterator,
4697  std::pair<std::vector<Point<dim>>, std::vector<unsigned int>>,
4699  cell_qpoint_map;
4700 
4701  // Now the easy case.
4702  if (n_points == 0)
4703  return cell_qpoint_map;
4704 
4705  // We begin by finding the cell/transform of the first point
4706  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4707  Point<dim>>
4708  point_and_reference_location;
4709 
4710  unsigned int counter = 0;
4711 
4712  while (counter < n_points)
4713  try
4714  {
4715  unsigned int i = counter;
4716  ++counter;
4717 
4718  point_and_reference_location =
4719  GridTools::find_active_cell_around_point(cache, points[i]);
4720  break;
4721  }
4722  catch (...)
4723  {
4724  if (counter == n_points)
4725  return cell_qpoint_map;
4726  }
4727 
4728  auto last_cell = cell_qpoint_map.emplace(std::make_pair(
4729  point_and_reference_location.first,
4730  std::make_pair(
4731  std::vector<Point<dim>>{point_and_reference_location.second},
4732  std::vector<unsigned int>{counter - 1})));
4733 
4734  // Now the second easy case.
4735  if (n_points == 1)
4736  return cell_qpoint_map;
4737 
4738  Point<spacedim> cell_center =
4739  point_and_reference_location.first->center();
4740  double cell_diameter = point_and_reference_location.first->diameter() *
4742 
4743  // Cycle over all points left
4744  for (unsigned int p = counter; p < n_points; ++p)
4745  {
4746  // Checking if the point is close to the cell center, in which
4747  // case calling find active cell with a cell hint
4748  if (cell_center.distance(points[p]) < cell_diameter)
4749  try
4750  {
4751  point_and_reference_location =
4753  cache, points[p], last_cell.first->first);
4754  }
4755  catch (...)
4756  {
4757  continue;
4758  }
4759  else
4760  try
4761  {
4762  point_and_reference_location =
4763  GridTools::find_active_cell_around_point(cache, points[p]);
4764  }
4765  catch (...)
4766  {
4767  continue;
4768  }
4769 
4770  if (last_cell.first->first == point_and_reference_location.first)
4771  {
4772  last_cell.first->second.first.emplace_back(
4773  point_and_reference_location.second);
4774  last_cell.first->second.second.emplace_back(p);
4775  }
4776  else
4777  {
4778  // Check if it is in another cell already found
4779  last_cell = cell_qpoint_map.emplace(
4780  std::make_pair(point_and_reference_location.first,
4781  std::make_pair(
4782  std::vector<Point<dim>>{
4783  point_and_reference_location.second},
4784  std::vector<unsigned int>{p})));
4785 
4786  if (last_cell.second == false)
4787  {
4788  // Cell already present: adding the new point
4789  last_cell.first->second.first.emplace_back(
4790  point_and_reference_location.second);
4791  last_cell.first->second.second.emplace_back(p);
4792  }
4793  else
4794  {
4795  // New cell was added, updating center and diameter
4796  cell_center = point_and_reference_location.first->center();
4797  cell_diameter =
4798  point_and_reference_location.first->diameter() *
4800  }
4801  }
4802  }
4803 
4804 #ifdef DEBUG
4805  unsigned int inserted_points = 0;
4806  // The number of points in all
4807  // the cells must be the same as
4808  // the number of points we
4809  // started off from.
4810  for (const auto &map_entry : cell_qpoint_map)
4811  {
4812  Assert(map_entry.second.second.size() ==
4813  map_entry.second.first.size(),
4814  ExcDimensionMismatch(map_entry.second.second.size(),
4815  map_entry.second.first.size()));
4816  inserted_points += map_entry.second.second.size();
4817  }
4818 #endif
4819  return cell_qpoint_map;
4820  }
4821 
4822 
4823 
4824  // Merge the input data to the existing map point_locations. If the cell
4825  // is already present in the map add information about the new points.
4826  // If the cell is not present add the cell with all information.
4827  //
4828  // Notice we call "information" the data associated with a point of the
4829  // sort: containing cell, coordinates on reference cell, index,
4830  // rank of the owner etc.
4831  template <int dim, int spacedim>
4832  void
4834  const std::vector<
4835  typename Triangulation<dim, spacedim>::active_cell_iterator> &cells,
4836  const std::vector<std::vector<Point<dim>>> & qpoints,
4837  const std::vector<std::vector<unsigned int>> & maps,
4838  const std::vector<std::vector<Point<spacedim>>> & points,
4839  const unsigned int rank,
4840  std::unordered_map<
4841  typename Triangulation<dim, spacedim>::active_cell_iterator,
4842  std::tuple<std::vector<Point<dim>>,
4843  std::vector<unsigned int>,
4844  std::vector<Point<spacedim>>,
4845  std::vector<unsigned int>>,
4846  cell_hash<dim, spacedim>> &point_locations)
4847  {
4848  // Adding cells
4849  for (unsigned int c = 0; c < cells.size(); ++c)
4850  {
4851  // Attempt to add a new cell with its relative data
4852  auto current_c = point_locations.emplace(
4853  std::make_pair(cells[c],
4854  std::make_tuple(qpoints[c],
4855  maps[c],
4856  points[c],
4857  std::vector<unsigned int>(
4858  points[c].size(), rank))));
4859 
4860  // If the flag is false the cell already existed
4861  if (current_c.second == false)
4862  {
4863  // Add the information to the cell at current_c.first:
4864  auto &cell_qpts = std::get<0>(current_c.first->second);
4865  auto &cell_maps = std::get<1>(current_c.first->second);
4866  auto &cell_pts = std::get<2>(current_c.first->second);
4867  auto &cell_ranks = std::get<3>(current_c.first->second);
4868 
4869  cell_qpts.insert(cell_qpts.end(),
4870  qpoints[c].begin(),
4871  qpoints[c].end());
4872  cell_maps.insert(cell_maps.end(),
4873  maps[c].begin(),
4874  maps[c].end());
4875  cell_pts.insert(cell_pts.end(),
4876  points[c].begin(),
4877  points[c].end());
4878  std::vector<unsigned int> ranks_tmp(points[c].size(), rank);
4879  cell_ranks.insert(cell_ranks.end(),
4880  ranks_tmp.begin(),
4881  ranks_tmp.end());
4882  }
4883  }
4884  }
4885 
4886 
4887 
4888  // This function calls compute point locations for all local_points
4889  // and sorts them in those which are probably locally owned, this which
4890  // are probably in ghost cells, and dismisses those in artificial cells
4891  // Output quantities are:
4892  // - locally_owned_locations: points, with relative information, inside
4893  // locally owned
4894  // cells
4895  // - ghost_cell_locations: points, with relative information, inside ghost
4896  // cells
4897  // - classified pts: indices of all points returned in
4898  // locally_owned_locations and
4899  // ghost_cell_locations (dropping those that were not found)
4900  template <int dim, int spacedim>
4901  void
4903  const GridTools::Cache<dim, spacedim> &cache,
4904  const std::vector<Point<spacedim>> & local_points,
4905  const std::vector<unsigned int> & local_points_idx,
4906  std::unordered_map<
4907  typename Triangulation<dim, spacedim>::active_cell_iterator,
4908  std::tuple<std::vector<Point<dim>>,
4909  std::vector<unsigned int>,
4910  std::vector<Point<spacedim>>,
4911  std::vector<unsigned int>>,
4912  cell_hash<dim, spacedim>> &locally_owned_locations,
4913  std::map<unsigned int,
4914  std::tuple<std::vector<CellId>,
4915  std::vector<std::vector<Point<dim>>>,
4916  std::vector<std::vector<unsigned int>>,
4917  std::vector<std::vector<Point<spacedim>>>>>
4918  & ghost_cell_locations,
4919  std::vector<unsigned int> &found_location_indices)
4920  {
4921  auto point_location_data =
4923  cache, local_points);
4924 
4925  // Sort output into locally owned cells, ghost cells, and artificial
4926  // cells.
4927  for (const auto &cell_tuples : point_location_data)
4928  {
4929  auto &cell = cell_tuples.first;
4930  auto &q_loc = std::get<0>(cell_tuples.second);
4931  auto &indices_loc = std::get<1>(cell_tuples.second);
4932 
4933  // Store the data for points in locally owned cells
4934  if (cell->is_locally_owned())
4935  {
4936  std::vector<Point<spacedim>> cell_points(indices_loc.size());
4937  std::vector<unsigned int> cell_points_idx(indices_loc.size());
4938  for (unsigned int i = 0; i < indices_loc.size(); ++i)
4939  {
4940  // Adding the point to the cell points
4941  cell_points[i] = local_points[indices_loc[i]];
4942 
4943  // Storing the index: notice indices loc refer to the local
4944  // points vector, but we need to return the index with
4945  // respect of the points owned by the current process
4946  cell_points_idx[i] = local_points_idx[indices_loc[i]];
4947  found_location_indices.emplace_back(
4948  local_points_idx[indices_loc[i]]);
4949  }
4950  locally_owned_locations.emplace(
4951  std::make_pair(cell,
4952  std::make_tuple(q_loc,
4953  cell_points_idx,
4954  cell_points,
4955  std::vector<unsigned int>(
4956  indices_loc.size(),
4957  cell->subdomain_id()))));
4958  }
4959  // Store the data for points in ghost cells and prepare transfer
4960  else if (cell->is_ghost())
4961  {
4962  std::vector<Point<spacedim>> cell_points(indices_loc.size());
4963  std::vector<unsigned int> cell_points_idx(indices_loc.size());
4964  for (unsigned int i = 0; i < indices_loc.size(); ++i)
4965  {
4966  cell_points[i] = local_points[indices_loc[i]];
4967  cell_points_idx[i] = local_points_idx[indices_loc[i]];
4968  found_location_indices.emplace_back(
4969  local_points_idx[indices_loc[i]]);
4970  }
4971  // Each key of the following map represents a process,
4972  // each mapped value is a tuple containing the information to be
4973  // sent: preparing the output for the owner, which has rank
4974  // subdomain id
4975  auto &map_tuple_owner =
4976  ghost_cell_locations[cell->subdomain_id()];
4977  // To identify the cell on the other process we use the cell id
4978  std::get<0>(map_tuple_owner).emplace_back(cell->id());
4979  std::get<1>(map_tuple_owner).emplace_back(q_loc);
4980  std::get<2>(map_tuple_owner).emplace_back(cell_points_idx);
4981  std::get<3>(map_tuple_owner).emplace_back(cell_points);
4982  }
4983  // else: the cell is artificial, nothing to do
4984  }
4985  }
4986 
4987 
4988 
4989  // Given the map received_point_locations obtained from a communication,
4990  // where the key is rank and the mapped value is a pair of
4991  // (points,indices), calls compute_point_locations; its output is then
4992  // merged with output tuple. If check_owned is set to true only points
4993  // lying inside locally owned cells are merged, otherwise all points are
4994  // merged into point_locations.
4995  template <int dim, int spacedim>
4996  void
4998  const GridTools::Cache<dim, spacedim> &cache,
4999  const std::map<
5000  unsigned int,
5001  std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>>
5002  &received_point_locations,
5003  std::unordered_map<
5004  typename Triangulation<dim, spacedim>::active_cell_iterator,
5005  std::tuple<std::vector<Point<dim>>,
5006  std::vector<unsigned int>,
5007  std::vector<Point<spacedim>>,
5008  std::vector<unsigned int>>,
5009  cell_hash<dim, spacedim>> &point_locations,
5010  const bool check_owned)
5011  {
5012  // rank and points is a pair: first rank, then a pair of vectors
5013  // (points, indices)
5014  for (const auto &rank_and_points : received_point_locations)
5015  {
5016  // Rewriting the contents of the map in human readable format
5017  const auto &received_process = rank_and_points.first;
5018  const auto &received_points = rank_and_points.second.first;
5019  const auto &received_map = rank_and_points.second.second;
5020 
5021  // Initializing the vectors needed to store the result of compute
5022  // point location
5023  std::vector<
5024  typename Triangulation<dim, spacedim>::active_cell_iterator>
5025  in_cell;
5026  std::vector<std::vector<Point<dim>>> in_qpoints;
5027  std::vector<std::vector<unsigned int>> in_maps;
5028  std::vector<std::vector<Point<spacedim>>> in_points;
5029 
5030  const auto computed_point_locations =
5032  compute_point_locations(cache, rank_and_points.second.first);
5033  for (const auto &map_c_pt_idx : computed_point_locations)
5034  {
5035  // Human-readable variables:
5036  const auto &proc_cell = map_c_pt_idx.first;
5037  const auto &proc_qpoints = map_c_pt_idx.second.first;
5038  const auto &proc_maps = map_c_pt_idx.second.second;
5039 
5040  // store either if we're not checking if the cell is
5041  // owned or if the cell is locally owned
5042  if (check_owned == false || proc_cell->is_locally_owned())
5043  {
5044  in_cell.emplace_back(proc_cell);
5045  in_qpoints.emplace_back(proc_qpoints);
5046  // The other two vectors need to be built
5047  unsigned int loc_size = proc_qpoints.size();
5048  std::vector<unsigned int> cell_maps(loc_size);
5049  std::vector<Point<spacedim>> cell_points(loc_size);
5050  for (unsigned int pt = 0; pt < loc_size; ++pt)
5051  {
5052  cell_maps[pt] = received_map[proc_maps[pt]];
5053  cell_points[pt] = received_points[proc_maps[pt]];
5054  }
5055  in_maps.emplace_back(cell_maps);
5056  in_points.emplace_back(cell_points);
5057  }
5058  }
5059 
5060  // Merge everything from the current process
5063  in_qpoints,
5064  in_maps,
5065  in_points,
5066  received_process,
5067  point_locations);
5068  }
5069  }
5070  } // namespace DistributedComputePointLocations
5071  } // namespace internal
5072 
5073 
5074 
5075  template <int dim, int spacedim>
5076 #ifndef DOXYGEN
5077  std::tuple<
5078  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5079  std::vector<std::vector<Point<dim>>>,
5080  std::vector<std::vector<unsigned int>>,
5081  std::vector<std::vector<Point<spacedim>>>,
5082  std::vector<std::vector<unsigned int>>>
5083 #else
5084  return_type
5085 #endif
5087  const GridTools::Cache<dim, spacedim> & cache,
5088  const std::vector<Point<spacedim>> & local_points,
5089  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes)
5090  {
5091 #ifndef DEAL_II_WITH_MPI
5092  (void)cache;
5093  (void)local_points;
5094  (void)global_bboxes;
5095  Assert(false,
5096  ExcMessage(
5097  "GridTools::distributed_compute_point_locations() requires MPI."));
5098  std::tuple<
5099  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5100  std::vector<std::vector<Point<dim>>>,
5101  std::vector<std::vector<unsigned int>>,
5102  std::vector<std::vector<Point<spacedim>>>,
5103  std::vector<std::vector<unsigned int>>>
5104  tup;
5105  return tup;
5106 #else
5107  // Recovering the mpi communicator used to create the triangulation
5108  const auto &tria_mpi =
5109  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
5110  &cache.get_triangulation());
5111  // If the dynamic cast failed we can't recover the mpi communicator:
5112  // throwing an assertion error
5113  Assert(
5114  tria_mpi,
5115  ExcMessage(
5116  "GridTools::distributed_compute_point_locations() requires a parallel triangulation."));
5117  auto mpi_communicator = tria_mpi->get_communicator();
5118  // Preparing the output tuple
5119  std::tuple<
5120  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5121  std::vector<std::vector<Point<dim>>>,
5122  std::vector<std::vector<unsigned int>>,
5123  std::vector<std::vector<Point<spacedim>>>,
5124  std::vector<std::vector<unsigned int>>>
5125  output_tuple;
5126 
5127  // Preparing the map that will be filled with found points
5128  std::unordered_map<
5130  std::tuple<std::vector<Point<dim>>,
5131  std::vector<unsigned int>,
5132  std::vector<Point<spacedim>>,
5133  std::vector<unsigned int>>,
5135  found_points;
5136 
5137  // Step 1 (part 1): Using the bounding boxes to guess the owner of each
5138  // point in local_points
5139  const unsigned int my_rank =
5140  Utilities::MPI::this_mpi_process(mpi_communicator);
5141 
5142  // Using global bounding boxes to guess/find owner/s of each point
5143  std::tuple<std::vector<std::vector<unsigned int>>,
5144  std::map<unsigned int, unsigned int>,
5145  std::map<unsigned int, std::vector<unsigned int>>>
5146  guessed_points;
5147  guessed_points = GridTools::guess_point_owner(global_bboxes, local_points);
5148 
5149  // Preparing to call compute_point_locations on points which may be local
5150  const auto &guess_loc_idx = std::get<0>(guessed_points)[my_rank];
5151  const unsigned int n_local_guess = guess_loc_idx.size();
5152 
5153  // Vector containing points which are probably local
5154  std::vector<Point<spacedim>> guess_local_points(n_local_guess);
5155  for (unsigned int i = 0; i < n_local_guess; ++i)
5156  guess_local_points[i] = local_points[guess_loc_idx[i]];
5157 
5158  // Preparing the map with data on points lying on ghost cells
5159  std::map<unsigned int,
5160  std::tuple<std::vector<CellId>,
5161  std::vector<std::vector<Point<dim>>>,
5162  std::vector<std::vector<unsigned int>>,
5163  std::vector<std::vector<Point<spacedim>>>>>
5164  found_ghost_points;
5165 
5166  // Vector containing indices of points lying either on locally owned
5167  // cells or ghost cells, to avoid computing them more than once
5168  std::vector<unsigned int> found_point_indices;
5169 
5170  // Thread used to call compute point locations on guess local pts
5171  Threads::Task<void> compute_locations_task =
5172  Threads::new_task(&internal::DistributedComputePointLocations::
5173  compute_and_classify_points<dim, spacedim>,
5174  cache,
5175  guess_local_points,
5176  guess_loc_idx,
5177  found_points,
5178  found_ghost_points,
5179  found_point_indices);
5180 
5181  // Step 1 (part 2): communicate points which are owned by a certain process
5182  // Preparing the map with points whose owner is known with certainty:
5183  const auto &not_locally_owned_idx = std::get<1>(guessed_points);
5184  std::map<unsigned int,
5185  std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>>
5186  not_locally_owned_points;
5187 
5188  for (const auto &indices : not_locally_owned_idx)
5189  if (indices.second != my_rank)
5190  {
5191  // Finding the list of points to be sent to this rank
5192  auto &points_to_send = not_locally_owned_points[indices.second];
5193  // Indices.first is the index of the considered point in local points
5194  points_to_send.first.emplace_back(local_points[indices.first]);
5195  points_to_send.second.emplace_back(indices.first);
5196  }
5197 
5198  // Communicating the points whose owner is sure
5199  auto received_points =
5200  Utilities::MPI::some_to_some(mpi_communicator, not_locally_owned_points);
5201  // Waiting for part 1 to finish to avoid concurrency problems
5202  compute_locations_task.join();
5203 
5204  // Step 2 (part 1): merge received points which are owned by us
5205  Threads::Task<void> merge_locally_owned_points_task =
5206  Threads::new_task(&internal::DistributedComputePointLocations::
5207  merge_received_point_locations<dim, spacedim>,
5208  cache,
5209  received_points,
5210  found_points,
5211  false);
5212 
5213  // Step 2 (part 2): communicate info on points lying on ghost cells
5214  auto received_ghost_points =
5215  Utilities::MPI::some_to_some(mpi_communicator, found_ghost_points);
5216 
5217  // Step 3: construct vectors containing points with uncertain owner i.e.
5218  // those which have multiple guesses. The map goes from rank of the probable
5219  // owner to a pair of vectors: the first containing the points, the second
5220  // containing the ranks in the current process
5221  std::map<unsigned int,
5222  std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>>
5223  uncertain_points;
5224 
5225  // This map goes from the point index to a vector of
5226  // ranks of probable owners
5227  const std::map<unsigned int, std::vector<unsigned int>>
5228  &points_to_probable_owners = std::get<2>(guessed_points);
5229 
5230  // Points in found_point_indices need not to be communicated;
5231  // sorting the array classified pts in order to use
5232  // binary search when checking if the points needs to be
5233  // communicated
5234  // Note that found_point_indices is a vector of integer indexes
5235  std::sort(found_point_indices.begin(), found_point_indices.end());
5236 
5237  for (const auto &probable_owners : points_to_probable_owners)
5238  {
5239  const auto &point_idx = probable_owners.first;
5240  const auto &probable_owner_ranks = probable_owners.second;
5241  if (!std::binary_search(found_point_indices.begin(),
5242  found_point_indices.end(),
5243  point_idx))
5244  // The point wasn't found in ghost or locally owned cells: send it
5245  for (const unsigned int probable_owner_rank : probable_owner_ranks)
5246  if (probable_owner_rank != my_rank)
5247  {
5248  // add to the data for probable_owner_rank
5249  auto &points_to_send = uncertain_points[probable_owner_rank];
5250  points_to_send.first.emplace_back(local_points[point_idx]);
5251  points_to_send.second.emplace_back(point_idx);
5252  }
5253  }
5254 
5255  // Step 4: send around uncertain points
5256  const auto received_uncertain_points =
5257  Utilities::MPI::some_to_some(mpi_communicator, uncertain_points);
5258  // Before proceeding, merging threads to avoid concurrency problems
5259  merge_locally_owned_points_task.join();
5260 
5261  // Step 5: add the received ghost cell data to output
5262  for (const auto &received_ghost_point : received_ghost_points)
5263  {
5264  // Transforming CellsIds into Tria iterators
5265  const auto &cell_ids = std::get<0>(received_ghost_point.second);
5266  const unsigned int n_cells = cell_ids.size();
5267  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5268  cell_iter(n_cells);
5269  for (unsigned int c = 0; c < n_cells; ++c)
5270  cell_iter[c] = cell_ids[c].to_cell(cache.get_triangulation());
5271 
5273  cell_iter,
5274  std::get<1>(received_ghost_point.second),
5275  std::get<2>(received_ghost_point.second),
5276  std::get<3>(received_ghost_point.second),
5277  received_ghost_point.first,
5278  found_points);
5279  }
5280 
5281  // Step 6: use compute point locations on the uncertain points and
5282  // merge output
5284  cache, received_uncertain_points, found_points, true);
5285 
5286  // Copying data from the unordered map to the tuple
5287  // and returning output
5288  const unsigned int size_output = found_points.size();
5289  auto &out_cells = std::get<0>(output_tuple);
5290  auto &out_qpoints = std::get<1>(output_tuple);
5291  auto &out_maps = std::get<2>(output_tuple);
5292  auto &out_points = std::get<3>(output_tuple);
5293  auto &out_ranks = std::get<4>(output_tuple);
5294 
5295  out_cells.resize(size_output);
5296  out_qpoints.resize(size_output);
5297  out_maps.resize(size_output);
5298  out_points.resize(size_output);
5299  out_ranks.resize(size_output);
5300 
5301  unsigned int c = 0;
5302  for (const auto &cell_and_data : found_points)
5303  {
5304  out_cells[c] = cell_and_data.first;
5305  out_qpoints[c] = std::get<0>(cell_and_data.second);
5306  out_maps[c] = std::get<1>(cell_and_data.second);
5307  out_points[c] = std::get<2>(cell_and_data.second);
5308  out_ranks[c] = std::get<3>(cell_and_data.second);
5309  ++c;
5310  }
5311 
5312  return output_tuple;
5313 #endif
5314  }
5315 
5316 
5317  template <int dim, int spacedim>
5318  std::map<unsigned int, Point<spacedim>>
5320  const Mapping<dim, spacedim> & mapping)
5321  {
5322  std::map<unsigned int, Point<spacedim>> result;
5323  for (const auto &cell : container.active_cell_iterators())
5324  {
5325  if (!cell->is_artificial())
5326  {
5327  const auto vs = mapping.get_vertices(cell);
5328  for (unsigned int i = 0; i < vs.size(); ++i)
5329  result[cell->vertex_index(i)] = vs[i];
5330  }
5331  }
5332  return result;
5333  }
5334 
5335 
5336  template <int spacedim>
5337  unsigned int
5338  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
5339  const Point<spacedim> & p)
5340  {
5341  auto id_and_v = std::min_element(
5342  vertices.begin(),
5343  vertices.end(),
5344  [&](const std::pair<const unsigned int, Point<spacedim>> &p1,
5345  const std::pair<const unsigned int, Point<spacedim>> &p2) -> bool {
5346  return p1.second.distance(p) < p2.second.distance(p);
5347  });
5348  return id_and_v->first;
5349  }
5350 
5351 
5352  template <int dim, int spacedim>
5353  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5354  Point<dim>>
5356  const Cache<dim, spacedim> &cache,
5357  const Point<spacedim> & p,
5359  & cell_hint,
5360  const std::vector<bool> &marked_vertices,
5361  const double tolerance)
5362  {
5363  const auto &mesh = cache.get_triangulation();
5364  const auto &mapping = cache.get_mapping();
5365  const auto &vertex_to_cells = cache.get_vertex_to_cell_map();
5366  const auto &vertex_to_cell_centers =
5368  const auto &used_vertices_rtree = cache.get_used_vertices_rtree();
5369 
5370  return find_active_cell_around_point(mapping,
5371  mesh,
5372  p,
5373  vertex_to_cells,
5374  vertex_to_cell_centers,
5375  cell_hint,
5376  marked_vertices,
5377  used_vertices_rtree,
5378  tolerance);
5379  }
5380 
5381  template <int spacedim>
5382  std::vector<std::vector<BoundingBox<spacedim>>>
5384  const std::vector<BoundingBox<spacedim>> &local_bboxes,
5385  const MPI_Comm & mpi_communicator)
5386  {
5387 #ifndef DEAL_II_WITH_MPI
5388  (void)local_bboxes;
5389  (void)mpi_communicator;
5390  Assert(false,
5391  ExcMessage(
5392  "GridTools::exchange_local_bounding_boxes() requires MPI."));
5393  return {};
5394 #else
5395  // Step 1: preparing data to be sent
5396  unsigned int n_bboxes = local_bboxes.size();
5397  // Dimension of the array to be exchanged (number of double)
5398  int n_local_data = 2 * spacedim * n_bboxes;
5399  // data array stores each entry of each point describing the bounding boxes
5400  std::vector<double> loc_data_array(n_local_data);
5401  for (unsigned int i = 0; i < n_bboxes; ++i)
5402  for (unsigned int d = 0; d < spacedim; ++d)
5403  {
5404  // Extracting the coordinates of each boundary point
5405  loc_data_array[2 * i * spacedim + d] =
5406  local_bboxes[i].get_boundary_points().first[d];
5407  loc_data_array[2 * i * spacedim + spacedim + d] =
5408  local_bboxes[i].get_boundary_points().second[d];
5409  }
5410 
5411  // Step 2: exchanging the size of local data
5412  unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_communicator);
5413 
5414  // Vector to store the size of loc_data_array for every process
5415  std::vector<int> size_all_data(n_procs);
5416 
5417  // Exchanging the number of bboxes
5418  int ierr = MPI_Allgather(&n_local_data,
5419  1,
5420  MPI_INT,
5421  size_all_data.data(),
5422  1,
5423  MPI_INT,
5424  mpi_communicator);
5425  AssertThrowMPI(ierr);
5426 
5427  // Now computing the the displacement, relative to recvbuf,
5428  // at which to store the incoming data
5429  std::vector<int> rdispls(n_procs);
5430  rdispls[0] = 0;
5431  for (unsigned int i = 1; i < n_procs; ++i)
5432  rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
5433 
5434  // Step 3: exchange the data and bounding boxes:
5435  // Allocating a vector to contain all the received data
5436  std::vector<double> data_array(rdispls.back() + size_all_data.back());
5437 
5438  ierr = MPI_Allgatherv(loc_data_array.data(),
5439  n_local_data,
5440  MPI_DOUBLE,
5441  data_array.data(),
5442  size_all_data.data(),
5443  rdispls.data(),
5444  MPI_DOUBLE,
5445  mpi_communicator);
5446  AssertThrowMPI(ierr);
5447 
5448  // Step 4: create the array of bboxes for output
5449  std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
5450  unsigned int begin_idx = 0;
5451  for (unsigned int i = 0; i < n_procs; ++i)
5452  {
5453  // Number of local bounding boxes
5454  unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
5455  global_bboxes[i].resize(n_bbox_i);
5456  for (unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
5457  {
5458  Point<spacedim> p1, p2; // boundary points for bbox
5459  for (unsigned int d = 0; d < spacedim; ++d)
5460  {
5461  p1[d] = data_array[begin_idx + 2 * bbox * spacedim + d];
5462  p2[d] =
5463  data_array[begin_idx + 2 * bbox * spacedim + spacedim + d];
5464  }
5465  BoundingBox<spacedim> loc_bbox(std::make_pair(p1, p2));
5466  global_bboxes[i][bbox] = loc_bbox;
5467  }
5468  // Shifting the first index to the start of the next vector
5469  begin_idx += size_all_data[i];
5470  }
5471  return global_bboxes;
5472 #endif // DEAL_II_WITH_MPI
5473  }
5474 
5475 
5476 
5477  template <int spacedim>
5480  const std::vector<BoundingBox<spacedim>> &local_description,
5481  const MPI_Comm & mpi_communicator)
5482  {
5483 #ifndef DEAL_II_WITH_MPI
5484  (void)mpi_communicator;
5485  // Building a tree with the only boxes available without MPI
5486  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> boxes_index(
5487  local_description.size());
5488  // Adding to each box the rank of the process owning it
5489  for (unsigned int i = 0; i < local_description.size(); ++i)
5490  boxes_index[i] = std::make_pair(local_description[i], 0u);
5491  return pack_rtree(boxes_index);
5492 #else
5493  // Exchanging local bounding boxes
5494  const std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes =
5495  Utilities::MPI::all_gather(mpi_communicator, local_description);
5496 
5497  // Preparing to flatten the vector
5498  const unsigned int n_procs =
5499  Utilities::MPI::n_mpi_processes(mpi_communicator);
5500  // The i'th element of the following vector contains the index of the first
5501  // local bounding box from the process of rank i
5502  std::vector<unsigned int> bboxes_position(n_procs);
5503 
5504  unsigned int tot_bboxes = 0;
5505  for (const auto &process_bboxes : global_bboxes)
5506  tot_bboxes += process_bboxes.size();
5507 
5508  // Now flattening the vector
5509  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
5510  flat_global_bboxes;
5511  flat_global_bboxes.reserve(tot_bboxes);
5512  unsigned int process_index = 0;
5513  for (const auto &process_bboxes : global_bboxes)
5514  {
5515  // Initialize a vector containing bounding boxes and rank of a process
5516  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
5517  boxes_and_indices(process_bboxes.size());
5518 
5519  // Adding to each box the rank of the process owning it
5520  for (unsigned int i = 0; i < process_bboxes.size(); ++i)
5521  boxes_and_indices[i] =
5522  std::make_pair(process_bboxes[i], process_index);
5523 
5524  flat_global_bboxes.insert(flat_global_bboxes.end(),
5525  boxes_and_indices.begin(),
5526  boxes_and_indices.end());
5527 
5528  ++process_index;
5529  }
5530 
5531  // Build a tree out of the bounding boxes. We avoid using the
5532  // insert method so that boost uses the packing algorithm
5533  return RTree<std::pair<BoundingBox<spacedim>, unsigned int>>(
5534  flat_global_bboxes.begin(), flat_global_bboxes.end());
5535 #endif // DEAL_II_WITH_MPI
5536  }
5537 
5538 
5539 
5540  template <int dim, int spacedim>
5541  void
5543  const Triangulation<dim, spacedim> & tria,
5544  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
5545  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group)
5546  {
5547  // 1) determine for each vertex a vertex it concides with and
5548  // put it into a map
5549  {
5550  static const int lookup_table_2d[2][2] =
5551  // flip:
5552  {
5553  {0, 1}, // false
5554  {1, 0} // true
5555  };
5556 
5557  static const int lookup_table_3d[2][2][2][4] =
5558  // orientation flip rotation
5559  {{{
5560  {0, 2, 1, 3}, // false false false
5561  {2, 3, 0, 1} // false false true
5562  },
5563  {
5564  {3, 1, 2, 0}, // false true false
5565  {1, 0, 3, 2} // false true true
5566  }},
5567  {{
5568  {0, 1, 2, 3}, // true false false
5569  {1, 3, 0, 2} // true false true
5570  },
5571  {
5572  {3, 2, 1, 0}, // true true false
5573  {2, 0, 3, 1} // true true true
5574  }}};
5575 
5576  // loop over all periodic face pairs
5577  for (const auto &pair : tria.get_periodic_face_map())
5578  {
5579  if (pair.first.first->level() != pair.second.first.first->level())
5580  continue;
5581 
5582  const auto face_a = pair.first.first->face(pair.first.second);
5583  const auto face_b =
5584  pair.second.first.first->face(pair.second.first.second);
5585  const auto mask = pair.second.second;
5586 
5587  AssertDimension(face_a->n_vertices(), face_b->n_vertices());
5588 
5589  // loop over all vertices on face
5590  for (unsigned int i = 0; i < face_a->n_vertices(); ++i)
5591  {
5592  const bool face_orientation = mask[0];
5593  const bool face_flip = mask[1];
5594  const bool face_rotation = mask[2];
5595 
5596  // find the right local vertex index for the second face
5597  unsigned int j = 0;
5598  switch (dim)
5599  {
5600  case 1:
5601  j = i;
5602  break;
5603  case 2:
5604  j = lookup_table_2d[face_flip][i];
5605  break;
5606  case 3:
5607  j = lookup_table_3d[face_orientation][face_flip]
5608  [face_rotation][i];
5609  break;
5610  default:
5611  AssertThrow(false, ExcNotImplemented());
5612  }
5613 
5614  // get vertex indices and store in map
5615  const auto vertex_a = face_a->vertex_index(i);
5616  const auto vertex_b = face_b->vertex_index(j);
5617  unsigned int temp = std::min(vertex_a, vertex_b);
5618 
5619  auto it_a = vertex_to_coinciding_vertex_group.find(vertex_a);
5620  if (it_a != vertex_to_coinciding_vertex_group.end())
5621  temp = std::min(temp, it_a->second);
5622 
5623  auto it_b = vertex_to_coinciding_vertex_group.find(vertex_b);
5624  if (it_b != vertex_to_coinciding_vertex_group.end())
5625  temp = std::min(temp, it_b->second);
5626 
5627  if (it_a != vertex_to_coinciding_vertex_group.end())
5628  it_a->second = temp;
5629  else
5630  vertex_to_coinciding_vertex_group[vertex_a] = temp;
5631 
5632  if (it_b != vertex_to_coinciding_vertex_group.end())
5633  it_b->second = temp;
5634  else
5635  vertex_to_coinciding_vertex_group[vertex_b] = temp;
5636  }
5637  }
5638 
5639  // 2) compress map: let vertices point to the coinciding vertex with
5640  // the smallest index
5641  for (auto &p : vertex_to_coinciding_vertex_group)
5642  {
5643  if (p.first == p.second)
5644  continue;
5645  unsigned int temp = p.second;
5646  while (temp != vertex_to_coinciding_vertex_group[temp])
5647  temp = vertex_to_coinciding_vertex_group[temp];
5648  p.second = temp;
5649  }
5650 
5651  // 3) create a map: smallest index of coinciding index -> all
5652  // coinciding indices
5653  for (auto p : vertex_to_coinciding_vertex_group)
5654  coinciding_vertex_groups[p.second] = {};
5655 
5656  for (auto p : vertex_to_coinciding_vertex_group)
5657  coinciding_vertex_groups[p.second].push_back(p.first);
5658  }
5659  }
5660 
5661 
5662 
5663  template <int dim, int spacedim>
5664  std::map<unsigned int, std::set<::types::subdomain_id>>
5666  const Triangulation<dim, spacedim> &tria)
5667  {
5668  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
5669  &tria) == nullptr) // nothing to do for a serial triangulation
5670  return {};
5671 
5672  // 1) collect for each vertex on periodic faces all vertices it coincides
5673  // with
5674  std::map<unsigned int, std::vector<unsigned int>> coinciding_vertex_groups;
5675  std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
5676 
5678  coinciding_vertex_groups,
5679  vertex_to_coinciding_vertex_group);
5680 
5681  // 2) collect vertices belonging to local cells
5682  std::vector<bool> vertex_of_own_cell(tria.n_vertices(), false);
5683  for (const auto &cell : tria.active_cell_iterators())
5684  if (cell->is_locally_owned())
5685  for (const unsigned int v : cell->vertex_indices())
5686  vertex_of_own_cell[cell->vertex_index(v)] = true;
5687 
5688  // 3) for each vertex belonging to a locally owned cell all ghost
5689  // neighbors (including the periodic own)
5690  std::map<unsigned int, std::set<types::subdomain_id>> result;
5691 
5692  // loop over all active ghost cells
5693  for (const auto &cell : tria.active_cell_iterators())
5694  if (cell->is_ghost())
5695  {
5696  const types::subdomain_id owner = cell->subdomain_id();
5697 
5698  // loop over all its vertices
5699  for (const unsigned int v : cell->vertex_indices())
5700  {
5701  // set owner if vertex belongs to a local cell
5702  if (vertex_of_own_cell[cell->vertex_index(v)])
5703  result[cell->vertex_index(v)].insert(owner);
5704 
5705  // mark also nodes coinciding due to periodicity
5706  auto coinciding_vertex_group =
5707  vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
5708  if (coinciding_vertex_group !=
5709  vertex_to_coinciding_vertex_group.end())
5710  for (auto coinciding_vertex :
5711  coinciding_vertex_groups[coinciding_vertex_group->second])
5712  if (vertex_of_own_cell[coinciding_vertex])
5713  result[coinciding_vertex].insert(owner);
5714  }
5715  }
5716 
5717  return result;
5718  }
5719 
5720 } /* namespace GridTools */
5721 
5722 
5723 // explicit instantiations
5724 #include "grid_tools.inst"
5725 
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:3917
void map_boundary_to_manifold_ids(const std::vector< types::boundary_id > &src_boundary_ids, const std::vector< types::manifold_id > &dst_manifold_ids, Triangulation< dim, spacedim > &tria, const std::vector< types::boundary_id > &reset_boundary_ids={})
Definition: grid_tools.cc:3702
std::vector< CellData< 1 > > boundary_lines
Transformed quadrature weights.
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, const MPI_Comm &mpi_communicator)
Definition: grid_tools.cc:5479
void laplace_transform(const std::map< unsigned int, Point< dim >> &new_points, Triangulation< dim > &tria, const Function< dim, double > *coefficient=nullptr, const bool solve_for_absolute_positions=false)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
unsigned int n_active_cells() const
Definition: tria.cc:12036
void insert_face_data(const FaceIteratorType &)
Definition: grid_tools.cc:515
const Triangulation< dim, spacedim > & get_triangulation() const
unsigned int n_vertices() const
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
static void reorder_cells(std::vector< CellData< dim >> &original_cells, const bool use_new_style_ordering=false)
const types::manifold_id flat_manifold_id
Definition: types.h:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
unsigned int manifold_id
Definition: types.h:141
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:1094
double objective_function(const Iterator &object, const Point< spacedim > &object_mid_point)
Definition: grid_tools.cc:3221
double diameter(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Mapping< dim, spacedim > &mapping)
Definition: grid_tools.cc:3099
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:3677
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:2082
active_face_iterator begin_active_face() const
Definition: tria.cc:11603
void create_laplace_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, SparseMatrix< double > &matrix, const Function< spacedim > *const a=nullptr, const AffineConstraints< double > &constraints=AffineConstraints< double >())
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:78
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1444
Rotate3d(const double angle, const unsigned int axis)
Definition: grid_tools.cc:883
GridTools::compute_local_to_global_vertex_index_map.
Definition: mpi_tags.h:105
virtual bool has_hanging_nodes() const
Definition: tria.cc:12168
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:309
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: grid_tools.cc:5319
Point< 3 > operator()(const Point< 3 > &p) const
Definition: grid_tools.cc:889
unsigned int n_cells() const
Definition: tria.cc:12028
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:3885
const Mapping< dim, spacedim > & get_mapping() const
static const double KA[GeometryInfo< dim >::vertices_per_cell][dim]
Definition: grid_tools.cc:210
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:391
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Task< RT > new_task(const std::function< RT()> &function)
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:3979
void add(const size_type i, const size_type j)
Volume element.
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:951
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:134
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:11546
Point< spacedim > operator()(const Point< spacedim > p) const
Definition: grid_tools.cc:869
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2180
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
std::vector< unsigned int > vertices
void join() const
Shift(const Tensor< 1, spacedim > &shift)
Definition: grid_tools.cc:865
bool compare_point_association(const unsigned int a, const unsigned int b, const Tensor< 1, spacedim > &point_direction, const std::vector< Tensor< 1, spacedim >> &center_directions)
Definition: grid_tools.cc:1654
return_type distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &local_points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes)
Definition: grid_tools.cc:5086
LinearOperator< Range, Domain, Payload > constrained_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11374
static const char U
std::map< unsigned int, types::global_vertex_index > compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2229
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
Point< 2 > second
Definition: grid_out.cc:4341
std::size_t operator()(const typename Triangulation< dim, spacedim >::active_cell_iterator &k) const
Definition: grid_tools.cc:4671
void fix_up_faces(const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, std::integral_constant< int, dim >, std::integral_constant< int, spacedim >)
Definition: grid_tools.cc:3591
types::boundary_id boundary_id
std::tuple< BoundingBox< MeshType::space_dimension >, bool > compute_cell_predicate_bounding_box(const typename MeshType::cell_iterator &parent_cell, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
Definition: grid_tools.cc:1882
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
const DerivativeForm< 1, dim, spacedim > & jacobian(const unsigned int quadrature_point) const
static const double Kb[GeometryInfo< dim >::vertices_per_cell]
Definition: grid_tools.cc:211
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11354
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3168
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_cell_bounding_boxes_rtree() const
void insert_face_data(const FaceIteratorType &face)
Definition: grid_tools.cc:476
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3013
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3139
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition: rtree.h:145
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
unsigned int n_levels() const
void merge_into_point_locations(const std::vector< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells, const std::vector< std::vector< Point< dim >>> &qpoints, const std::vector< std::vector< unsigned int >> &maps, const std::vector< std::vector< Point< spacedim >>> &points, const unsigned int rank, std::unordered_map< typename Triangulation< dim, spacedim >::active_cell_iterator, std::tuple< std::vector< Point< dim >>, std::vector< unsigned int >, std::vector< Point< spacedim >>, std::vector< unsigned int >>, cell_hash< dim, spacedim >> &point_locations)
Definition: grid_tools.cc:4833
const double angle
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:2681
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
Definition: tria.cc:9580
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:385
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, const MPI_Comm &mpi_communicator)
Definition: grid_tools.cc:5383
static double distance_to_unit_cell(const Point< dim > &p)
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:623
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2616
cell_iterator end() const
Definition: tria.cc:11440
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:531
size_type n() const
bool operator()(const CellData< structdim > &a, const CellData< structdim > &b) const
Definition: grid_tools.cc:428
virtual void execute_coarsening_and_refinement()
Definition: tria.cc:12703
void set_subdomain_id_in_zorder_recursively(IT cell, unsigned int &current_proc_idx, unsigned int &current_cell_idx, const unsigned int n_active_cells, const unsigned int n_partitions)
Definition: grid_tools.cc:2879
RTree< typename LeafTypeIterator::value_type, IndexType, IndexableGetter > pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
Definition: tria.cc:11557
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
bool check_consistency(const unsigned int dim) const
Definition: fe_q.h:548
unsigned int subdomain_id
Definition: types.h:43
Scale(const double factor)
Definition: grid_tools.cc:914
T sum(const T &t, const MPI_Comm &mpi_communicator)
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2645
Expression acos(const Expression &x)
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
Definition: tria.cc:9849
#define Assert(cond, exc)
Definition: exceptions.h:1466
Signals signals
Definition: tria.h:2272
IteratorRange< active_cell_iterator > active_cell_iterators() const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
types::global_dof_index n_dofs() const
Abstract base class for mapping classes.
Definition: mapping.h:301
std::list< typename Triangulation< dim, spacedim >::cell_iterator > distorted_cells
Definition: tria.h:1538
unsigned int n_quads() const
Definition: tria.cc:12393
bool fix_up_object(const Iterator &object)
Definition: grid_tools.cc:3409
const Tensor< 1, spacedim > shift
Definition: grid_tools.cc:875
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim - dim, spacedim >(&forms)[vertices_per_cell])
std::vector< BoundingBox< MeshType::space_dimension > > compute_mesh_predicate_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const unsigned int refinement_level=0, const bool allow_merge=false, const unsigned int max_boxes=numbers::invalid_unsigned_int)
Definition: grid_tools.cc:1935
void save_user_indices(std::vector< unsigned int > &v) const
Definition: tria.cc:10907
types::material_id material_id
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
void load_user_indices(const std::vector< unsigned int > &v)
Definition: tria.cc:10939
unsigned int level
Definition: grid_out.cc:4343
unsigned int n_lines() const
Definition: tria.cc:12180
const RTree< std::pair< Point< spacedim >, unsigned int > > & get_used_vertices_rtree() const
VectorType::value_type * end(VectorType &V)
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={}, const double tolerance=1.e-10)
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:3950
std::string to_string(const T &t)
Definition: patterns.h:2342
Point< 3 > vertices[4]
double minimal_diameter(const Iterator &object)
Definition: grid_tools.cc:3382
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1419
void collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int >> &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
Definition: grid_tools.cc:5542
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim >> &vertices)
Definition: grid_tools.cc:285
uint64_t global_vertex_index
Definition: types.h:48
Expression fabs(const Expression &x)
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:3769
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:11986
Triangulation< dim, spacedim >::DistortedCellList fix_up_distorted_child_cells(const typename Triangulation< dim, spacedim >::DistortedCellList &distorted_cells, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3635
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
void rotate(const double angle, Triangulation< dim > &triangulation)
PackagedOperation< Range > constrained_right_hand_side(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop, const Range &right_hand_side)
void append_face_data(const CellData< 1 > &face_data, SubCellData &subcell_data)
Definition: grid_tools.cc:408
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:11979
unsigned int size() const
virtual const MPI_Comm & get_communicator() const
Definition: tria_base.cc:139
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Point< 2 > first
Definition: grid_out.cc:4340
number singular_value(const size_type i) const
types::manifold_id manifold_id
const std::vector< std::vector< Tensor< 1, spacedim > > > & get_vertex_to_cell_centers_directions() const
unsigned int n_raw_faces() const
Definition: tria.cc:12071
std::vector< std::vector< Tensor< 1, spacedim > > > vertex_to_cell_centers_directions(const Triangulation< dim, spacedim > &mesh, const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator >> &vertex_to_cells)
Definition: grid_tools.cc:1615
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
static const char A
void reorder_hierarchical(const DynamicSparsityPattern &sparsity, std::vector< DynamicSparsityPattern::size_type > &new_indices)
Point< Iterator::AccessorType::space_dimension > get_face_midpoint(const Iterator &object, const unsigned int f, std::integral_constant< int, 3 >)
Definition: grid_tools.cc:3349
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
__global__ void set(Number *val, const Number s, const size_type N)
std::set< CellData< dim - 1 >, internal::CellDataComparator< dim - 1 > > face_data
Definition: grid_tools.cc:504
return_type compute_point_locations_try_all(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:4330
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
GridTools::compute_local_to_global_vertex_index_map second tag.
Definition: mpi_tags.h:107
const unsigned int axis
Definition: grid_tools.cc:907
void advance(std::tuple< I1, I2 > &t, const unsigned int n)
void distribute(VectorType &vec) const
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1747
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
static constexpr double PI
Definition: numbers.h:231
const std::vector< bool > & get_used_vertices() const
Definition: tria.cc:12580
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:422
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
VectorType::value_type * begin(VectorType &V)
void laplace_solve(const SparseMatrix< double > &S, const AffineConstraints< double > &constraints, Vector< double > &u)
Definition: grid_tools.cc:967
T min(const T &t, const MPI_Comm &mpi_communicator)
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
Definition: grid_tools.cc:1126
std::vector< typename MeshType< dim, spacedim >::active_cell_iterator > find_cells_adjacent_to_vertex(const MeshType< dim, spacedim > &container, const unsigned int vertex_index)
Definition: grid_tools.cc:1483
std::vector< CellData< 2 > > boundary_quads
numbers::NumberTraits< Number >::real_type square() const
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:375
static const char N
void get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, std::vector< types::subdomain_id > &subdomain)
Definition: grid_tools.cc:3040
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2582
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void refine_global(const unsigned int times=1)
Definition: tria.cc:10198
Point< spacedim > operator()(const Point< spacedim > p) const
Definition: grid_tools.cc:918