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