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