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