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