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