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