Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2001 - 2021 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>
21
26
30
32#include <deal.II/fe/fe_q.h>
37
43#include <deal.II/grid/tria.h>
46
54#include <deal.II/lac/vector.h>
56
59
61#include <boost/random/mersenne_twister.hpp>
62#include <boost/random/uniform_real_distribution.hpp>
64
65#include <array>
66#include <cmath>
67#include <iostream>
68#include <list>
69#include <numeric>
70#include <set>
71#include <tuple>
72#include <unordered_map>
73
75
76
77namespace GridTools
78{
79 template <int dim, int spacedim>
80 double
82 {
83 // we can't deal with distributed meshes since we don't have all
84 // vertices locally. there is one exception, however: if the mesh has
85 // never been refined. the way to test this is not to ask
86 // tria.n_levels()==1, since this is something that can happen on one
87 // processor without being true on all. however, we can ask for the
88 // global number of active cells and use that
89#if defined(DEAL_II_WITH_P4EST) && defined(DEBUG)
91 dynamic_cast<
93 Assert(p_tria->n_global_active_cells() == tria.n_cells(0),
95#endif
96
97 // the algorithm used simply traverses all cells and picks out the
98 // boundary vertices. it may or may not be faster to simply get all
99 // vectors, don't mark boundary vertices, and compute the distances
100 // thereof, but at least as the mesh is refined, it seems better to
101 // first mark boundary nodes, as marking is O(N) in the number of
102 // cells/vertices, while computing the maximal distance is O(N*N)
103 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
104 std::vector<bool> boundary_vertices(vertices.size(), false);
105
107 tria.begin_active();
109 tria.end();
110 for (; cell != endc; ++cell)
111 for (const unsigned int face : cell->face_indices())
112 if (cell->face(face)->at_boundary())
113 for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
114 boundary_vertices[cell->face(face)->vertex_index(i)] = true;
115
116 // now traverse the list of boundary vertices and check distances.
117 // since distances are symmetric, we only have to check one half
118 double max_distance_sqr = 0;
119 std::vector<bool>::const_iterator pi = boundary_vertices.begin();
120 const unsigned int N = boundary_vertices.size();
121 for (unsigned int i = 0; i < N; ++i, ++pi)
122 {
123 std::vector<bool>::const_iterator pj = pi + 1;
124 for (unsigned int j = i + 1; j < N; ++j, ++pj)
125 if ((*pi == true) && (*pj == true) &&
126 ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
127 max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
128 }
129
130 return std::sqrt(max_distance_sqr);
131 }
132
133
134
135 template <int dim, int spacedim>
136 double
138 const Mapping<dim, spacedim> & mapping)
139 {
140 // get the degree of the mapping if possible. if not, just assume 1
141 unsigned int mapping_degree = 1;
142 if (const auto *p =
143 dynamic_cast<const MappingQGeneric<dim, spacedim> *>(&mapping))
144 mapping_degree = p->get_degree();
145 else if (const auto *p =
146 dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
147 mapping_degree = p->get_degree();
148
149 // then initialize an appropriate quadrature formula
150 const QGauss<dim> quadrature_formula(mapping_degree + 1);
151 const unsigned int n_q_points = quadrature_formula.size();
152
153 // we really want the JxW values from the FEValues object, but it
154 // wants a finite element. create a cheap element as a dummy
155 // element
157 FEValues<dim, spacedim> fe_values(mapping,
158 dummy_fe,
159 quadrature_formula,
161
163 cell = triangulation.begin_active(),
164 endc = triangulation.end();
165
166 double local_volume = 0;
167
168 // compute the integral quantities by quadrature
169 for (; cell != endc; ++cell)
170 if (cell->is_locally_owned())
171 {
172 fe_values.reinit(cell);
173 for (unsigned int q = 0; q < n_q_points; ++q)
174 local_volume += fe_values.JxW(q);
175 }
176
177 double global_volume = 0;
178
179#ifdef DEAL_II_WITH_MPI
181 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
183 global_volume =
184 Utilities::MPI::sum(local_volume, p_tria->get_communicator());
185 else
186#endif
187 global_volume = local_volume;
188
189 return global_volume;
190 }
191
192
193
194 namespace
195 {
210 template <int dim>
211 struct TransformR2UAffine
212 {
213 static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
215 };
216
217
218 /*
219 Octave code:
220 M=[0 1; 1 1];
221 K1 = transpose(M) * inverse (M*transpose(M));
222 printf ("{%f, %f},\n", K1' );
223 */
224 template <>
226 [1] = {{-1.000000}, {1.000000}};
227
228 template <>
230 {1.000000, 0.000000};
231
232
233 /*
234 Octave code:
235 M=[0 1 0 1;0 0 1 1;1 1 1 1];
236 K2 = transpose(M) * inverse (M*transpose(M));
237 printf ("{%f, %f, %f},\n", K2' );
238 */
239 template <>
241 [2] = {{-0.500000, -0.500000},
242 {0.500000, -0.500000},
243 {-0.500000, 0.500000},
244 {0.500000, 0.500000}};
245
246 /*
247 Octave code:
248 M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
249 K3 = transpose(M) * inverse (M*transpose(M))
250 printf ("{%f, %f, %f, %f},\n", K3' );
251 */
252 template <>
254 {0.750000, 0.250000, 0.250000, -0.250000};
255
256
257 template <>
259 [3] = {
260 {-0.250000, -0.250000, -0.250000},
261 {0.250000, -0.250000, -0.250000},
262 {-0.250000, 0.250000, -0.250000},
263 {0.250000, 0.250000, -0.250000},
264 {-0.250000, -0.250000, 0.250000},
265 {0.250000, -0.250000, 0.250000},
266 {-0.250000, 0.250000, 0.250000},
267 {0.250000, 0.250000, 0.250000}
268
269 };
270
271
272 template <>
274 {0.500000,
275 0.250000,
276 0.250000,
277 0.000000,
278 0.250000,
279 0.000000,
280 0.000000,
281 -0.250000};
282 } // namespace
283
284
285
286 template <int dim, int spacedim>
287 std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
289 {
291
292 // A = vertex * KA
294
295 for (unsigned int d = 0; d < spacedim; ++d)
296 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
297 for (unsigned int e = 0; e < dim; ++e)
298 A[d][e] += vertices[v][d] * TransformR2UAffine<dim>::KA[v][e];
299
300 // b = vertex * Kb
302 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
304
305 return std::make_pair(A, b);
306 }
307
308
309
310 template <int dim>
314 const Quadrature<dim> & quadrature)
315 {
317 FEValues<dim> fe_values(mapping, fe, quadrature, update_jacobians);
318
319 Vector<double> aspect_ratio_vector(triangulation.n_active_cells());
320
321 // loop over cells of processor
322 for (const auto &cell : triangulation.active_cell_iterators())
323 {
324 if (cell->is_locally_owned())
325 {
326 double aspect_ratio_cell = 0.0;
327
328 fe_values.reinit(cell);
329
330 // loop over quadrature points
331 for (unsigned int q = 0; q < quadrature.size(); ++q)
332 {
333 const Tensor<2, dim, double> jacobian =
334 Tensor<2, dim, double>(fe_values.jacobian(q));
335
336 // We intentionally do not want to throw an exception in case of
337 // inverted elements since this is not the task of this
338 // function. Instead, inf is written into the vector in case of
339 // inverted elements.
340 if (determinant(jacobian) <= 0)
341 {
342 aspect_ratio_cell = std::numeric_limits<double>::infinity();
343 }
344 else
345 {
347 for (unsigned int i = 0; i < dim; i++)
348 for (unsigned int j = 0; j < dim; j++)
349 J(i, j) = jacobian[i][j];
350
351 J.compute_svd();
352
353 double const max_sv = J.singular_value(0);
354 double const min_sv = J.singular_value(dim - 1);
355 double const ar = max_sv / min_sv;
356
357 // Take the max between the previous and the current
358 // aspect ratio value; if we had previously encountered
359 // an inverted cell, we will have placed an infinity
360 // in the aspect_ratio_cell variable, and that value
361 // will survive this max operation.
362 aspect_ratio_cell = std::max(aspect_ratio_cell, ar);
363 }
364 }
365
366 // fill vector
367 aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
368 }
369 }
370
371 return aspect_ratio_vector;
372 }
373
374
375
376 template <int dim>
377 double
380 const Quadrature<dim> & quadrature)
381 {
382 Vector<double> aspect_ratio_vector =
383 compute_aspect_ratio_of_cells(mapping, triangulation, quadrature);
384
386 aspect_ratio_vector,
388 }
389
390
391
392 template <int dim, int spacedim>
395 {
396 using iterator =
398 const auto predicate = [](const iterator &) { return true; };
399
401 tria, std::function<bool(const iterator &)>(predicate));
402 }
403
404
405
406 // Generic functions for appending face data in 2D or 3D. TODO: we can
407 // remove these once we have 'if constexpr'.
408 namespace internal
409 {
410 inline void
411 append_face_data(const CellData<1> &face_data, SubCellData &subcell_data)
412 {
413 subcell_data.boundary_lines.push_back(face_data);
414 }
415
416
417
418 inline void
419 append_face_data(const CellData<2> &face_data, SubCellData &subcell_data)
420 {
421 subcell_data.boundary_quads.push_back(face_data);
422 }
423
424
425
426 // Lexical comparison for sorting CellData objects.
427 template <int structdim>
429 {
430 bool
432 const CellData<structdim> &b) const
433 {
434 // Check vertices:
435 if (std::lexicographical_compare(std::begin(a.vertices),
437 std::begin(b.vertices),
438 std::end(b.vertices)))
439 return true;
440 // it should never be necessary to check the material or manifold
441 // ids as a 'tiebreaker' (since they must be equal if the vertex
442 // indices are equal). Assert it anyway:
443#ifdef DEBUG
444 if (std::equal(std::begin(a.vertices),
446 std::begin(b.vertices)))
447 {
448 Assert(a.material_id == b.material_id &&
449 a.manifold_id == b.manifold_id,
451 "Two CellData objects with equal vertices must "
452 "have the same material/boundary ids and manifold "
453 "ids."));
454 }
455#endif
456 return false;
457 }
458 };
459
460
470 template <int dim>
472 {
473 public:
477 template <class FaceIteratorType>
478 void
479 insert_face_data(const FaceIteratorType &face)
480 {
481 CellData<dim - 1> face_cell_data;
482 for (unsigned int vertex_n = 0; vertex_n < face->n_vertices();
483 ++vertex_n)
484 face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
485 face_cell_data.boundary_id = face->boundary_id();
486 face_cell_data.manifold_id = face->manifold_id();
487
488 face_data.insert(face_cell_data);
489 }
490
496 {
497 SubCellData subcell_data;
498
499 for (const CellData<dim - 1> &face_cell_data : face_data)
500 internal::append_face_data(face_cell_data, subcell_data);
501 return subcell_data;
502 }
503
504
505 private:
508 };
509
510
511 // Do nothing for dim=1:
512 template <>
514 {
515 public:
516 template <class FaceIteratorType>
517 void
518 insert_face_data(const FaceIteratorType &)
519 {}
520
523 {
524 return SubCellData();
525 }
526 };
527 } // namespace internal
528
529
530
531 template <int dim, int spacedim>
532 std::
533 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
535 {
536 Assert(1 <= tria.n_levels(),
537 ExcMessage("The input triangulation must be non-empty."));
538
539 std::vector<Point<spacedim>> vertices;
540 std::vector<CellData<dim>> cells;
541
542 unsigned int max_level_0_vertex_n = 0;
543 for (const auto &cell : tria.cell_iterators_on_level(0))
544 for (const unsigned int cell_vertex_n : cell->vertex_indices())
545 max_level_0_vertex_n =
546 std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
547 vertices.resize(max_level_0_vertex_n + 1);
548
550 std::set<CellData<1>, internal::CellDataComparator<1>>
551 line_data; // only used in 3D
552
553 for (const auto &cell : tria.cell_iterators_on_level(0))
554 {
555 // Save cell data
556 CellData<dim> cell_data(cell->n_vertices());
557 for (const unsigned int cell_vertex_n : cell->vertex_indices())
558 {
559 Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
561 vertices[cell->vertex_index(cell_vertex_n)] =
562 cell->vertex(cell_vertex_n);
563 cell_data.vertices[cell_vertex_n] =
564 cell->vertex_index(cell_vertex_n);
565 }
566 cell_data.material_id = cell->material_id();
567 cell_data.manifold_id = cell->manifold_id();
568 cells.push_back(cell_data);
569
570 // Save face data
571 if (dim > 1)
572 {
573 for (const unsigned int face_n : cell->face_indices())
574 face_data.insert_face_data(cell->face(face_n));
575 }
576 // Save line data
577 if (dim == 3)
578 {
579 for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
580 {
581 const auto line = cell->line(line_n);
582 CellData<1> line_cell_data;
583 for (unsigned int vertex_n = 0; vertex_n < line->n_vertices();
584 ++vertex_n)
585 line_cell_data.vertices[vertex_n] =
586 line->vertex_index(vertex_n);
587 line_cell_data.boundary_id = line->boundary_id();
588 line_cell_data.manifold_id = line->manifold_id();
589
590 line_data.insert(line_cell_data);
591 }
592 }
593 }
594
595 // Double-check that there are no unused vertices:
596#ifdef DEBUG
597 {
598 std::vector<bool> used_vertices(vertices.size());
599 for (const CellData<dim> &cell_data : cells)
600 for (const auto v : cell_data.vertices)
601 used_vertices[v] = true;
602 Assert(std::find(used_vertices.begin(), used_vertices.end(), false) ==
603 used_vertices.end(),
604 ExcMessage("The level zero vertices should form a contiguous "
605 "range."));
606 }
607#endif
608
609 SubCellData subcell_data = face_data.get();
610
611 if (dim == 3)
612 for (const CellData<1> &face_line_data : line_data)
613 subcell_data.boundary_lines.push_back(face_line_data);
614
615 return std::tuple<std::vector<Point<spacedim>>,
616 std::vector<CellData<dim>>,
617 SubCellData>(std::move(vertices),
618 std::move(cells),
619 std::move(subcell_data));
620 }
621
622
623
624 template <int dim, int spacedim>
625 void
627 std::vector<CellData<dim>> & cells,
628 SubCellData & subcelldata)
629 {
630 Assert(
631 subcelldata.check_consistency(dim),
633 "Invalid SubCellData supplied according to ::check_consistency(). "
634 "This is caused by data containing objects for the wrong dimension."));
635
636 // first check which vertices are actually used
637 std::vector<bool> vertex_used(vertices.size(), false);
638 for (unsigned int c = 0; c < cells.size(); ++c)
639 for (unsigned int v = 0; v < cells[c].vertices.size(); ++v)
640 {
641 Assert(cells[c].vertices[v] < vertices.size(),
642 ExcMessage("Invalid vertex index encountered! cells[" +
643 Utilities::int_to_string(c) + "].vertices[" +
644 Utilities::int_to_string(v) + "]=" +
646 " is invalid, because only " +
648 " vertices were supplied."));
649 vertex_used[cells[c].vertices[v]] = true;
650 }
651
652
653 // then renumber the vertices that are actually used in the same order as
654 // they were beforehand
655 const unsigned int invalid_vertex = numbers::invalid_unsigned_int;
656 std::vector<unsigned int> new_vertex_numbers(vertices.size(),
657 invalid_vertex);
658 unsigned int next_free_number = 0;
659 for (unsigned int i = 0; i < vertices.size(); ++i)
660 if (vertex_used[i] == true)
661 {
662 new_vertex_numbers[i] = next_free_number;
663 ++next_free_number;
664 }
665
666 // next replace old vertex numbers by the new ones
667 for (unsigned int c = 0; c < cells.size(); ++c)
668 for (auto &v : cells[c].vertices)
669 v = new_vertex_numbers[v];
670
671 // same for boundary data
672 for (unsigned int c = 0; c < subcelldata.boundary_lines.size(); // NOLINT
673 ++c)
674 for (unsigned int v = 0;
675 v < subcelldata.boundary_lines[c].vertices.size();
676 ++v)
677 {
678 Assert(subcelldata.boundary_lines[c].vertices[v] <
679 new_vertex_numbers.size(),
681 "Invalid vertex index in subcelldata.boundary_lines. "
682 "subcelldata.boundary_lines[" +
683 Utilities::int_to_string(c) + "].vertices[" +
684 Utilities::int_to_string(v) + "]=" +
686 subcelldata.boundary_lines[c].vertices[v]) +
687 " is invalid, because only " +
689 " vertices were supplied."));
690 subcelldata.boundary_lines[c].vertices[v] =
691 new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]];
692 }
693
694 for (unsigned int c = 0; c < subcelldata.boundary_quads.size(); // NOLINT
695 ++c)
696 for (unsigned int v = 0;
697 v < subcelldata.boundary_quads[c].vertices.size();
698 ++v)
699 {
700 Assert(subcelldata.boundary_quads[c].vertices[v] <
701 new_vertex_numbers.size(),
703 "Invalid vertex index in subcelldata.boundary_quads. "
704 "subcelldata.boundary_quads[" +
705 Utilities::int_to_string(c) + "].vertices[" +
706 Utilities::int_to_string(v) + "]=" +
708 subcelldata.boundary_quads[c].vertices[v]) +
709 " is invalid, because only " +
711 " vertices were supplied."));
712
713 subcelldata.boundary_quads[c].vertices[v] =
714 new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]];
715 }
716
717 // finally copy over the vertices which we really need to a new array and
718 // replace the old one by the new one
719 std::vector<Point<spacedim>> tmp;
720 tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(), true));
721 for (unsigned int v = 0; v < vertices.size(); ++v)
722 if (vertex_used[v] == true)
723 tmp.push_back(vertices[v]);
724 swap(vertices, tmp);
725 }
726
727
728
729 template <int dim, int spacedim>
730 void
732 std::vector<CellData<dim>> & cells,
733 SubCellData & subcelldata,
734 std::vector<unsigned int> & considered_vertices,
735 const double tol)
736 {
737 AssertIndexRange(2, vertices.size());
738 // create a vector of vertex indices. initialize it to the identity, later
739 // on change that if necessary.
740 std::vector<unsigned int> new_vertex_numbers(vertices.size());
741 std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
742
743 // if the considered_vertices vector is empty, consider all vertices
744 if (considered_vertices.size() == 0)
745 considered_vertices = new_vertex_numbers;
746 Assert(considered_vertices.size() <= vertices.size(), ExcInternalError());
747
748 // The algorithm below improves upon the naive O(n^2) algorithm by first
749 // sorting vertices by their value in one component and then only
750 // comparing vertices for equality which are nearly equal in that
751 // component. For example, if @p vertices form a cube, then we will only
752 // compare points that have the same x coordinate when we try to find
753 // duplicated vertices.
754
755 // Start by finding the longest coordinate direction. This minimizes the
756 // number of points that need to be compared against each-other in a
757 // single set for typical geometries.
759 const auto & min = bbox.get_boundary_points().first;
760 const auto & max = bbox.get_boundary_points().second;
761
762 unsigned int longest_coordinate_direction = 0;
763 double longest_coordinate_length = max[0] - min[0];
764 for (unsigned int d = 1; d < spacedim; ++d)
765 {
766 const double coordinate_length = max[d] - min[d];
767 if (longest_coordinate_length < coordinate_length)
768 {
769 longest_coordinate_length = coordinate_length;
770 longest_coordinate_direction = d;
771 }
772 }
773
774 // Sort vertices (while preserving their vertex numbers) along that
775 // coordinate direction:
776 std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
777 sorted_vertices.reserve(vertices.size());
778 for (const unsigned int vertex_n : considered_vertices)
779 {
780 AssertIndexRange(vertex_n, vertices.size());
781 sorted_vertices.emplace_back(vertex_n, vertices[vertex_n]);
782 }
783 std::sort(sorted_vertices.begin(),
784 sorted_vertices.end(),
785 [&](const std::pair<unsigned int, Point<spacedim>> &a,
786 const std::pair<unsigned int, Point<spacedim>> &b) {
787 return a.second[longest_coordinate_direction] <
788 b.second[longest_coordinate_direction];
789 });
790
791 auto within_tolerance = [=](const Point<spacedim> &a,
792 const Point<spacedim> &b) {
793 for (unsigned int d = 0; d < spacedim; ++d)
794 if (std::abs(a[d] - b[d]) > tol)
795 return false;
796 return true;
797 };
798
799 // Find a range of numbers that have the same component in the longest
800 // coordinate direction:
801 auto range_start = sorted_vertices.begin();
802 while (range_start != sorted_vertices.end())
803 {
804 auto range_end = range_start + 1;
805 while (range_end != sorted_vertices.end() &&
806 std::abs(range_end->second[longest_coordinate_direction] -
807 range_start->second[longest_coordinate_direction]) <
808 tol)
809 ++range_end;
810
811 // preserve behavior with older versions of this function by replacing
812 // higher vertex numbers by lower vertex numbers
813 std::sort(range_start,
814 range_end,
815 [](const std::pair<unsigned int, Point<spacedim>> &a,
816 const std::pair<unsigned int, Point<spacedim>> &b) {
817 return a.first < b.first;
818 });
819
820 // Now de-duplicate [range_start, range_end)
821 //
822 // We have identified all points that are within a strip of width 'tol'
823 // in one coordinate direction. Now we need to figure out which of these
824 // are also close in other coordinate directions. If two are close, we
825 // can mark the second one for deletion.
826 for (auto reference = range_start; reference != range_end; ++reference)
827 {
828 if (reference->first != numbers::invalid_unsigned_int)
829 for (auto it = reference + 1; it != range_end; ++it)
830 {
831 if (within_tolerance(reference->second, it->second))
832 {
833 new_vertex_numbers[it->first] = reference->first;
834 // skip the replaced vertex in the future
836 }
837 }
838 }
839 range_start = range_end;
840 }
841
842 // now we got a renumbering list. simply renumber all vertices
843 // (non-duplicate vertices get renumbered to themselves, so nothing bad
844 // happens). after that, the duplicate vertices will be unused, so call
845 // delete_unused_vertices() to do that part of the job.
846 for (auto &cell : cells)
847 for (auto &vertex_index : cell.vertices)
848 vertex_index = new_vertex_numbers[vertex_index];
849 for (auto &quad : subcelldata.boundary_quads)
850 for (auto &vertex_index : quad.vertices)
851 vertex_index = new_vertex_numbers[vertex_index];
852 for (auto &line : subcelldata.boundary_lines)
853 for (auto &vertex_index : line.vertices)
854 vertex_index = new_vertex_numbers[vertex_index];
855
856 delete_unused_vertices(vertices, cells, subcelldata);
857 }
858
859
860
861 template <int dim, int spacedim>
862 void
864 const std::vector<Point<spacedim>> &all_vertices,
865 std::vector<CellData<dim>> & cells)
866 {
867 if (dim == 1)
868 return;
869 if (dim == 2 && spacedim == 3)
870 Assert(false, ExcNotImplemented());
871
872 std::size_t n_negative_cells = 0;
873 for (auto &cell : cells)
874 {
875 Assert(cell.vertices.size() ==
876 ReferenceCells::get_hypercube<dim>().n_vertices(),
878 const ArrayView<const unsigned int> vertices(cell.vertices);
879 if (GridTools::cell_measure(all_vertices, vertices) < 0)
880 {
881 ++n_negative_cells;
882
883 // TODO: this only works for quads and hexes
884 if (dim == 2)
885 {
886 // flip the cell across the y = x line in 2D
887 std::swap(cell.vertices[1], cell.vertices[2]);
888 }
889 else if (dim == 3)
890 {
891 // swap the front and back faces in 3D
892 std::swap(cell.vertices[0], cell.vertices[2]);
893 std::swap(cell.vertices[1], cell.vertices[3]);
894 std::swap(cell.vertices[4], cell.vertices[6]);
895 std::swap(cell.vertices[5], cell.vertices[7]);
896 }
897
898 // Check whether the resulting cell is now ok.
899 // If not, then the grid is seriously broken and
900 // we just give up.
903 }
904 }
905
906 // We assume that all cells of a grid have
907 // either positive or negative volumes but
908 // not both mixed. Although above reordering
909 // might work also on single cells, grids
910 // with both kind of cells are very likely to
911 // be broken. Check for this here.
912 AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
914 std::string(
915 "This function assumes that either all cells have positive "
916 "volume, or that all cells have been specified in an "
917 "inverted vertex order so that their volume is negative. "
918 "(In the latter case, this class automatically inverts "
919 "every cell.) However, the mesh you have specified "
920 "appears to have both cells with positive and cells with "
921 "negative volume. You need to check your mesh which "
922 "cells these are and how they got there.\n"
923 "As a hint, of the total ") +
924 std::to_string(cells.size()) + " cells in the mesh, " +
925 std::to_string(n_negative_cells) +
926 " appear to have a negative volume."));
927 }
928
929
930
931 // Functions and classes for consistently_order_cells
932 namespace
933 {
939 struct CheapEdge
940 {
944 CheapEdge(const unsigned int v0, const unsigned int v1)
945 : v0(v0)
946 , v1(v1)
947 {}
948
953 bool
954 operator<(const CheapEdge &e) const
955 {
956 return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1)));
957 }
958
959 private:
963 const unsigned int v0, v1;
964 };
965
966
975 template <int dim>
976 bool
977 is_consistent(const std::vector<CellData<dim>> &cells)
978 {
979 std::set<CheapEdge> edges;
980
981 for (typename std::vector<CellData<dim>>::const_iterator c =
982 cells.begin();
983 c != cells.end();
984 ++c)
985 {
986 // construct the edges in reverse order. for each of them,
987 // ensure that the reverse edge is not yet in the list of
988 // edges (return false if the reverse edge already *is* in
989 // the list) and then add the actual edge to it; std::set
990 // eliminates duplicates automatically
991 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
992 {
993 const CheapEdge reverse_edge(
996 if (edges.find(reverse_edge) != edges.end())
997 return false;
998
999
1000 // ok, not. insert edge in correct order
1001 const CheapEdge correct_edge(
1004 edges.insert(correct_edge);
1005 }
1006 }
1007
1008 // no conflicts found, so return true
1009 return true;
1010 }
1011
1012
1019 template <int dim>
1020 struct ParallelEdges
1021 {
1027 static const unsigned int starter_edges[dim];
1028
1033 static const unsigned int n_other_parallel_edges = (1 << (dim - 1)) - 1;
1034 static const unsigned int
1037 };
1038
1039 template <>
1040 const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2};
1041
1042 template <>
1043 const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1},
1044 {0},
1045 {3},
1046 {2}};
1047
1048 template <>
1049 const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8};
1050
1051 template <>
1052 const unsigned int ParallelEdges<3>::parallel_edges[12][3] = {
1053 {1, 4, 5}, // line 0
1054 {0, 4, 5}, // line 1
1055 {3, 6, 7}, // line 2
1056 {2, 6, 7}, // line 3
1057 {0, 1, 5}, // line 4
1058 {0, 1, 4}, // line 5
1059 {2, 3, 7}, // line 6
1060 {2, 3, 6}, // line 7
1061 {9, 10, 11}, // line 8
1062 {8, 10, 11}, // line 9
1063 {8, 9, 11}, // line 10
1064 {8, 9, 10} // line 11
1065 };
1066
1067
1072 struct AdjacentCell
1073 {
1077 AdjacentCell()
1080 {}
1081
1085 AdjacentCell(const unsigned int cell_index,
1086 const unsigned int edge_within_cell)
1089 {}
1090
1091
1092 unsigned int cell_index;
1093 unsigned int edge_within_cell;
1094 };
1095
1096
1097
1098 template <int dim>
1099 class AdjacentCells;
1100
1106 template <>
1107 class AdjacentCells<2>
1108 {
1109 public:
1114 using const_iterator = const AdjacentCell *;
1115
1124 void
1125 push_back(const AdjacentCell &adjacent_cell)
1126 {
1128 adjacent_cells[0] = adjacent_cell;
1129 else
1130 {
1134 adjacent_cells[1] = adjacent_cell;
1135 }
1136 }
1137
1138
1143 const_iterator
1144 begin() const
1145 {
1146 return adjacent_cells;
1147 }
1148
1149
1155 const_iterator
1156 end() const
1157 {
1158 // check whether the current object stores zero, one, or two
1159 // adjacent cells, and use this to point to the element past the
1160 // last valid one
1162 return adjacent_cells;
1164 return adjacent_cells + 1;
1165 else
1166 return adjacent_cells + 2;
1167 }
1168
1169 private:
1176 AdjacentCell adjacent_cells[2];
1177 };
1178
1179
1180
1188 template <>
1189 class AdjacentCells<3> : public std::vector<AdjacentCell>
1190 {};
1191
1192
1202 template <int dim>
1203 class Edge
1204 {
1205 public:
1211 Edge(const CellData<dim> &cell, const unsigned int edge_number)
1212 : orientation_status(not_oriented)
1213 {
1216
1217 // copy vertices for this particular line
1218 vertex_indices[0] =
1219 cell
1221 vertex_indices[1] =
1222 cell
1224
1225 // bring them into standard orientation
1226 if (vertex_indices[0] > vertex_indices[1])
1228 }
1229
1234 bool
1235 operator<(const Edge<dim> &e) const
1236 {
1237 return ((vertex_indices[0] < e.vertex_indices[0]) ||
1238 ((vertex_indices[0] == e.vertex_indices[0]) &&
1239 (vertex_indices[1] < e.vertex_indices[1])));
1240 }
1241
1245 bool
1246 operator==(const Edge<dim> &e) const
1247 {
1248 return ((vertex_indices[0] == e.vertex_indices[0]) &&
1249 (vertex_indices[1] == e.vertex_indices[1]));
1250 }
1251
1256 unsigned int vertex_indices[2];
1257
1262 enum OrientationStatus
1263 {
1264 not_oriented,
1265 forward,
1266 backward
1267 };
1268
1269 OrientationStatus orientation_status;
1270
1275 AdjacentCells<dim> adjacent_cells;
1276 };
1277
1278
1279
1284 template <int dim>
1285 struct Cell
1286 {
1292 Cell(const CellData<dim> &c, const std::vector<Edge<dim>> &edge_list)
1293 {
1294 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1295 vertex_indices[i] = c.vertices[i];
1296
1297 // now for each of the edges of this cell, find the location inside the
1298 // given edge_list array and store than index
1299 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1300 {
1301 const Edge<dim> e(c, l);
1302 edge_indices[l] =
1303 (std::lower_bound(edge_list.begin(), edge_list.end(), e) -
1304 edge_list.begin());
1305 Assert(edge_indices[l] < edge_list.size(), ExcInternalError());
1306 Assert(edge_list[edge_indices[l]] == e, ExcInternalError())
1307 }
1308 }
1309
1314
1320 };
1321
1322
1323
1324 template <int dim>
1325 class EdgeDeltaSet;
1326
1336 template <>
1337 class EdgeDeltaSet<2>
1338 {
1339 public:
1343 using const_iterator = const unsigned int *;
1344
1349 EdgeDeltaSet()
1350 {
1352 }
1353
1354
1358 void
1359 clear()
1360 {
1362 }
1363
1368 void
1369 insert(const unsigned int edge_index)
1370 {
1372 edge_indices[0] = edge_index;
1373 else
1374 {
1377 edge_indices[1] = edge_index;
1378 }
1379 }
1380
1381
1385 const_iterator
1386 begin() const
1387 {
1388 return edge_indices;
1389 }
1390
1391
1395 const_iterator
1396 end() const
1397 {
1398 // check whether the current object stores zero, one, or two
1399 // indices, and use this to point to the element past the
1400 // last valid one
1402 return edge_indices;
1404 return edge_indices + 1;
1405 else
1406 return edge_indices + 2;
1407 }
1408
1409 private:
1413 unsigned int edge_indices[2];
1414 };
1415
1416
1417
1429 template <>
1430 class EdgeDeltaSet<3> : public std::set<unsigned int>
1431 {};
1432
1433
1434
1439 template <int dim>
1440 std::vector<Edge<dim>>
1441 build_edges(const std::vector<CellData<dim>> &cells)
1442 {
1443 // build the edge list for all cells. because each cell has
1444 // GeometryInfo<dim>::lines_per_cell edges, the total number
1445 // of edges is this many times the number of cells. of course
1446 // some of them will be duplicates, and we throw them out below
1447 std::vector<Edge<dim>> edge_list;
1448 edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
1449 for (unsigned int i = 0; i < cells.size(); ++i)
1450 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1451 edge_list.emplace_back(cells[i], l);
1452
1453 // next sort the edge list and then remove duplicates
1454 std::sort(edge_list.begin(), edge_list.end());
1455 edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
1456 edge_list.end());
1457
1458 return edge_list;
1459 }
1460
1461
1462
1467 template <int dim>
1468 std::vector<Cell<dim>>
1469 build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
1470 std::vector<Edge<dim>> & edges)
1471 {
1472 std::vector<Cell<dim>> cell_list;
1473 cell_list.reserve(cells.size());
1474 for (unsigned int i = 0; i < cells.size(); ++i)
1475 {
1476 // create our own data structure for the cells and let it
1477 // connect to the edges array
1478 cell_list.emplace_back(cells[i], edges);
1479
1480 // then also inform the edges that they are adjacent
1481 // to the current cell, and where within this cell
1482 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1483 edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back(
1484 AdjacentCell(i, l));
1485 }
1486 Assert(cell_list.size() == cells.size(), ExcInternalError());
1487
1488 return cell_list;
1489 }
1490
1491
1492
1497 template <int dim>
1498 unsigned int
1499 get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
1500 const std::vector<Edge<dim>> &edges,
1501 const unsigned int current_cell)
1502 {
1503 for (unsigned int c = current_cell; c < cells.size(); ++c)
1504 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1505 if (edges[cells[c].edge_indices[l]].orientation_status ==
1506 Edge<dim>::not_oriented)
1507 return c;
1508
1510 }
1511
1512
1513
1519 template <int dim>
1520 void
1521 orient_one_set_of_parallel_edges(const std::vector<Cell<dim>> &cells,
1522 std::vector<Edge<dim>> & edges,
1523 const unsigned int cell,
1524 const unsigned int local_edge)
1525 {
1526 // choose the direction of the first edge. we have free choice
1527 // here and could simply choose "forward" if that's what pleases
1528 // us. however, for backward compatibility with the previous
1529 // implementation used till 2016, let us just choose the
1530 // direction so that it matches what we have in the given cell.
1531 //
1532 // in fact, in what can only be assumed to be a bug in the
1533 // original implementation, after orienting all edges, the code
1534 // that rotates the cells so that they match edge orientations
1535 // (see the rotate_cell() function below) rotated the cell two
1536 // more times by 90 degrees. this is ok -- it simply flips all
1537 // edge orientations, which leaves them valid. rather than do
1538 // the same in the current implementation, we can achieve the
1539 // same effect by modifying the rule above to choose the
1540 // direction of the starting edge of this parallel set
1541 // *opposite* to what it looks like in the current cell
1542 //
1543 // this bug only existed in the 2d implementation since there
1544 // were different implementations for 2d and 3d. consequently,
1545 // only replicate it for the 2d case and be "intuitive" in 3d.
1546 if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
1548 local_edge, 0)])
1549 // orient initial edge *opposite* to the way it is in the cell
1550 // (see above for the reason)
1551 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1552 (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
1553 else
1554 {
1555 Assert(
1556 edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
1557 cells[cell].vertex_indices
1560 Assert(
1561 edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
1562 cells[cell].vertex_indices
1565
1566 // orient initial edge *opposite* to the way it is in the cell
1567 // (see above for the reason)
1568 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1569 (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
1570 }
1571
1572 // walk outward from the given edge as described in
1573 // the algorithm in the paper that documents all of
1574 // this
1575 //
1576 // note that in 2d, each of the Deltas can at most
1577 // contain two elements, whereas in 3d it can be arbitrarily many
1578 EdgeDeltaSet<dim> Delta_k;
1579 EdgeDeltaSet<dim> Delta_k_minus_1;
1580 Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]);
1581
1582 while (Delta_k_minus_1.begin() !=
1583 Delta_k_minus_1.end()) // while set is not empty
1584 {
1585 Delta_k.clear();
1586
1587 for (typename EdgeDeltaSet<dim>::const_iterator delta =
1588 Delta_k_minus_1.begin();
1589 delta != Delta_k_minus_1.end();
1590 ++delta)
1591 {
1592 Assert(edges[*delta].orientation_status !=
1593 Edge<dim>::not_oriented,
1595
1596 // now go through the cells adjacent to this edge
1597 for (typename AdjacentCells<dim>::const_iterator adjacent_cell =
1598 edges[*delta].adjacent_cells.begin();
1599 adjacent_cell != edges[*delta].adjacent_cells.end();
1600 ++adjacent_cell)
1601 {
1602 const unsigned int K = adjacent_cell->cell_index;
1603 const unsigned int delta_is_edge_in_K =
1604 adjacent_cell->edge_within_cell;
1605
1606 // figure out the direction of delta with respect to the cell
1607 // K (in the orientation in which the user has given it to us)
1608 const unsigned int first_edge_vertex =
1609 (edges[*delta].orientation_status == Edge<dim>::forward ?
1610 edges[*delta].vertex_indices[0] :
1611 edges[*delta].vertex_indices[1]);
1612 const unsigned int first_edge_vertex_in_K =
1613 cells[K]
1615 delta_is_edge_in_K, 0)];
1616 Assert(
1617 first_edge_vertex == first_edge_vertex_in_K ||
1618 first_edge_vertex ==
1620 dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
1622
1623 // now figure out which direction the each of the "opposite"
1624 // edges needs to be oriented into.
1625 for (unsigned int o_e = 0;
1627 ++o_e)
1628 {
1629 // get the index of the opposite edge and select which its
1630 // first vertex needs to be based on how the current edge
1631 // is oriented in the current cell
1632 const unsigned int opposite_edge =
1633 cells[K].edge_indices[ParallelEdges<
1634 dim>::parallel_edges[delta_is_edge_in_K][o_e]];
1635 const unsigned int first_opposite_edge_vertex =
1636 cells[K].vertex_indices
1638 ParallelEdges<
1639 dim>::parallel_edges[delta_is_edge_in_K][o_e],
1640 (first_edge_vertex == first_edge_vertex_in_K ? 0 :
1641 1))];
1642
1643 // then determine the orientation of the edge based on
1644 // whether the vertex we want to be the edge's first
1645 // vertex is already the first vertex of the edge, or
1646 // whether it points in the opposite direction
1647 const typename Edge<dim>::OrientationStatus
1648 opposite_edge_orientation =
1649 (edges[opposite_edge].vertex_indices[0] ==
1650 first_opposite_edge_vertex ?
1651 Edge<dim>::forward :
1652 Edge<dim>::backward);
1653
1654 // see if the opposite edge (there is only one in 2d) has
1655 // already been oriented.
1656 if (edges[opposite_edge].orientation_status ==
1657 Edge<dim>::not_oriented)
1658 {
1659 // the opposite edge is not yet oriented. do orient it
1660 // and add it to Delta_k
1661 edges[opposite_edge].orientation_status =
1662 opposite_edge_orientation;
1663 Delta_k.insert(opposite_edge);
1664 }
1665 else
1666 {
1667 // this opposite edge has already been oriented. it
1668 // should be consistent with the current one in 2d,
1669 // while in 3d it may in fact be mis-oriented, and in
1670 // that case the mesh will not be orientable. indicate
1671 // this by throwing an exception that we can catch
1672 // further up; this has the advantage that we can
1673 // propagate through a couple of functions without
1674 // having to do error checking and without modifying
1675 // the 'cells' array that the user gave us
1676 if (dim == 2)
1677 {
1678 Assert(edges[opposite_edge].orientation_status ==
1679 opposite_edge_orientation,
1681 }
1682 else if (dim == 3)
1683 {
1684 if (edges[opposite_edge].orientation_status !=
1685 opposite_edge_orientation)
1686 throw ExcMeshNotOrientable();
1687 }
1688 else
1689 Assert(false, ExcNotImplemented());
1690 }
1691 }
1692 }
1693 }
1694
1695 // finally copy the new set to the previous one
1696 // (corresponding to increasing 'k' by one in the
1697 // algorithm)
1698 Delta_k_minus_1 = Delta_k;
1699 }
1700 }
1701
1702
1710 template <int dim>
1711 void
1712 rotate_cell(const std::vector<Cell<dim>> &cell_list,
1713 const std::vector<Edge<dim>> &edge_list,
1714 const unsigned int cell_index,
1715 std::vector<CellData<dim>> & raw_cells)
1716 {
1717 // find the first vertex of the cell. this is the vertex where dim edges
1718 // originate, so for each of the edges record which the starting vertex is
1719 unsigned int starting_vertex_of_edge[GeometryInfo<dim>::lines_per_cell];
1720 for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
1721 {
1722 Assert(edge_list[cell_list[cell_index].edge_indices[e]]
1723 .orientation_status != Edge<dim>::not_oriented,
1725 if (edge_list[cell_list[cell_index].edge_indices[e]]
1726 .orientation_status == Edge<dim>::forward)
1727 starting_vertex_of_edge[e] =
1728 edge_list[cell_list[cell_index].edge_indices[e]]
1729 .vertex_indices[0];
1730 else
1731 starting_vertex_of_edge[e] =
1732 edge_list[cell_list[cell_index].edge_indices[e]]
1733 .vertex_indices[1];
1734 }
1735
1736 // find the vertex number that appears dim times. this will then be
1737 // the vertex at which we want to locate the origin of the cell's
1738 // coordinate system (i.e., vertex 0)
1739 unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int;
1740 switch (dim)
1741 {
1742 case 2:
1743 {
1744 // in 2d, we can simply enumerate the possibilities where the
1745 // origin may be located because edges zero and one don't share
1746 // any vertices, and the same for edges two and three
1747 if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
1748 (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
1749 origin_vertex_of_cell = starting_vertex_of_edge[0];
1750 else if ((starting_vertex_of_edge[1] ==
1751 starting_vertex_of_edge[2]) ||
1752 (starting_vertex_of_edge[1] ==
1753 starting_vertex_of_edge[3]))
1754 origin_vertex_of_cell = starting_vertex_of_edge[1];
1755 else
1756 Assert(false, ExcInternalError());
1757
1758 break;
1759 }
1760
1761 case 3:
1762 {
1763 // one could probably do something similar in 3d, but that seems
1764 // more complicated than one wants to write down. just go
1765 // through the list of possible starting vertices and check
1766 for (origin_vertex_of_cell = 0;
1767 origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
1768 ++origin_vertex_of_cell)
1769 if (std::count(starting_vertex_of_edge,
1770 starting_vertex_of_edge +
1772 cell_list[cell_index]
1773 .vertex_indices[origin_vertex_of_cell]) == dim)
1774 break;
1775 Assert(origin_vertex_of_cell <
1778
1779 break;
1780 }
1781
1782 default:
1783 Assert(false, ExcNotImplemented());
1784 }
1785
1786 // now rotate raw_cells[cell_index] in such a way that its orientation
1787 // matches that of cell_list[cell_index]
1788 switch (dim)
1789 {
1790 case 2:
1791 {
1792 // in 2d, we can literally rotate the cell until its origin
1793 // matches the one that we have determined above should be
1794 // the origin vertex
1795 //
1796 // when doing a rotation, take into account the ordering of
1797 // vertices (not in clockwise or counter-clockwise sense)
1798 while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell)
1799 {
1800 const unsigned int tmp = raw_cells[cell_index].vertices[0];
1801 raw_cells[cell_index].vertices[0] =
1802 raw_cells[cell_index].vertices[1];
1803 raw_cells[cell_index].vertices[1] =
1804 raw_cells[cell_index].vertices[3];
1805 raw_cells[cell_index].vertices[3] =
1806 raw_cells[cell_index].vertices[2];
1807 raw_cells[cell_index].vertices[2] = tmp;
1808 }
1809 break;
1810 }
1811
1812 case 3:
1813 {
1814 // in 3d, the situation is a bit more complicated. from above, we
1815 // now know which vertex is at the origin (because 3 edges
1816 // originate from it), but that still leaves 3 possible rotations
1817 // of the cube. the important realization is that we can choose
1818 // any of them: in all 3 rotations, all edges originate from the
1819 // one vertex, and that fixes the directions of all 12 edges in
1820 // the cube because these 3 cover all 3 equivalence classes!
1821 // consequently, we can select an arbitrary one among the
1822 // permutations -- for example the following ones:
1823 static const unsigned int cube_permutations[8][8] = {
1824 {0, 1, 2, 3, 4, 5, 6, 7},
1825 {1, 5, 3, 7, 0, 4, 2, 6},
1826 {2, 6, 0, 4, 3, 7, 1, 5},
1827 {3, 2, 1, 0, 7, 6, 5, 4},
1828 {4, 0, 6, 2, 5, 1, 7, 3},
1829 {5, 4, 7, 6, 1, 0, 3, 2},
1830 {6, 7, 4, 5, 2, 3, 0, 1},
1831 {7, 3, 5, 1, 6, 2, 4, 0}};
1832
1833 unsigned int
1834 temp_vertex_indices[GeometryInfo<dim>::vertices_per_cell];
1835 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1836 temp_vertex_indices[v] =
1837 raw_cells[cell_index]
1838 .vertices[cube_permutations[origin_vertex_of_cell][v]];
1839 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1840 raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
1841
1842 break;
1843 }
1844
1845 default:
1846 {
1847 Assert(false, ExcNotImplemented());
1848 }
1849 }
1850 }
1851
1852
1858 template <int dim>
1859 void
1860 reorient(std::vector<CellData<dim>> &cells)
1861 {
1862 // first build the arrays that connect cells to edges and the other
1863 // way around
1864 std::vector<Edge<dim>> edge_list = build_edges(cells);
1865 std::vector<Cell<dim>> cell_list =
1866 build_cells_and_connect_edges(cells, edge_list);
1867
1868 // then loop over all cells and start orienting parallel edge sets
1869 // of cells that still have non-oriented edges
1870 unsigned int next_cell_with_unoriented_edge = 0;
1871 while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
1872 cell_list, edge_list, next_cell_with_unoriented_edge)) !=
1874 {
1875 // see which edge sets are still not oriented
1876 //
1877 // we do not need to look at each edge because if we orient edge
1878 // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there
1879 // will be 3 other edges that are also oriented). there are only
1880 // dim independent sets of edges, so loop over these.
1881 //
1882 // we need to check whether each one of these starter edges may
1883 // already be oriented because the line (sheet) that connects
1884 // globally parallel edges may be self-intersecting in the
1885 // current cell
1886 for (unsigned int l = 0; l < dim; ++l)
1887 if (edge_list[cell_list[next_cell_with_unoriented_edge]
1889 .orientation_status == Edge<dim>::not_oriented)
1890 orient_one_set_of_parallel_edges(
1891 cell_list,
1892 edge_list,
1893 next_cell_with_unoriented_edge,
1895
1896 // ensure that we have really oriented all edges now, not just
1897 // the starter edges
1898 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1899 Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
1900 .edge_indices[l]]
1901 .orientation_status != Edge<dim>::not_oriented,
1903 }
1904
1905 // now that we have oriented all edges, we need to rotate cells
1906 // so that the edges point in the right direction with the now
1907 // rotated coordinate system
1908 for (unsigned int c = 0; c < cells.size(); ++c)
1909 rotate_cell(cell_list, edge_list, c, cells);
1910 }
1911
1912
1913 // overload of the function above for 1d -- there is nothing
1914 // to orient in that case
1915 void reorient(std::vector<CellData<1>> &)
1916 {}
1917 } // namespace
1918
1919 template <int dim>
1920 void
1922 {
1923 Assert(cells.size() != 0,
1924 ExcMessage(
1925 "List of elements to orient must have at least one cell"));
1926
1927 // there is nothing for us to do in 1d
1928 if (dim == 1)
1929 return;
1930
1931 // check if grids are already consistent. if so, do
1932 // nothing. if not, then do the reordering
1933 if (!is_consistent(cells))
1934 try
1935 {
1936 reorient(cells);
1937 }
1938 catch (const ExcMeshNotOrientable &)
1939 {
1940 // the mesh is not orientable. this is acceptable if we are in 3d,
1941 // as class Triangulation knows how to handle this, but it is
1942 // not in 2d; in that case, re-throw the exception
1943 if (dim < 3)
1944 throw;
1945 }
1946 }
1947
1948
1949 // define some transformations
1950 namespace internal
1951 {
1952 template <int spacedim>
1953 class Shift
1954 {
1955 public:
1957 : shift(shift)
1958 {}
1961 {
1962 return p + shift;
1963 }
1964
1965 private:
1967 };
1968
1969
1970 // Transformation to rotate around one of the cartesian axes.
1972 {
1973 public:
1974 Rotate3d(const double angle, const unsigned int axis)
1975 : angle(angle)
1976 , axis(axis)
1977 {}
1978
1979 Point<3>
1980 operator()(const Point<3> &p) const
1981 {
1982 if (axis == 0)
1983 return {p(0),
1984 std::cos(angle) * p(1) - std::sin(angle) * p(2),
1985 std::sin(angle) * p(1) + std::cos(angle) * p(2)};
1986 else if (axis == 1)
1987 return {std::cos(angle) * p(0) + std::sin(angle) * p(2),
1988 p(1),
1989 -std::sin(angle) * p(0) + std::cos(angle) * p(2)};
1990 else
1991 return {std::cos(angle) * p(0) - std::sin(angle) * p(1),
1992 std::sin(angle) * p(0) + std::cos(angle) * p(1),
1993 p(2)};
1994 }
1995
1996 private:
1997 const double angle;
1998 const unsigned int axis;
1999 };
2000
2001 template <int spacedim>
2002 class Scale
2003 {
2004 public:
2005 explicit Scale(const double factor)
2006 : factor(factor)
2007 {}
2010 {
2011 return p * factor;
2012 }
2013
2014 private:
2015 const double factor;
2016 };
2017 } // namespace internal
2018
2019
2020 template <int dim, int spacedim>
2021 void
2022 shift(const Tensor<1, spacedim> & shift_vector,
2024 {
2026 }
2027
2028
2029 template <int dim>
2030 void
2031 rotate(const double angle,
2032 const unsigned int axis,
2034 {
2035 Assert(axis < 3, ExcMessage("Invalid axis given!"));
2036
2038 }
2039
2040 template <int dim, int spacedim>
2041 void
2042 scale(const double scaling_factor,
2044 {
2045 Assert(scaling_factor > 0, ExcScalingFactorNotPositive(scaling_factor));
2047 }
2048
2049
2050 namespace internal
2051 {
2057 inline void
2059 const AffineConstraints<double> &constraints,
2060 Vector<double> & u)
2061 {
2062 const unsigned int n_dofs = S.n();
2063 const auto op = linear_operator(S);
2064 const auto SF = constrained_linear_operator(constraints, op);
2066 prec.initialize(S, 1.2);
2067
2068 SolverControl control(n_dofs, 1.e-10, false, false);
2070 SolverCG<Vector<double>> solver(control, mem);
2071
2072 Vector<double> f(n_dofs);
2073
2074 const auto constrained_rhs =
2075 constrained_right_hand_side(constraints, op, f);
2076 solver.solve(SF, u, constrained_rhs, prec);
2077
2078 constraints.distribute(u);
2079 }
2080 } // namespace internal
2081
2082
2083 // Implementation for dimensions except 1
2084 template <int dim>
2085 void
2086 laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
2088 const Function<dim> * coefficient,
2089 const bool solve_for_absolute_positions)
2090 {
2091 if (dim == 1)
2092 Assert(false, ExcNotImplemented());
2093
2094 // first provide everything that is needed for solving a Laplace
2095 // equation.
2096 FE_Q<dim> q1(1);
2097
2098 DoFHandler<dim> dof_handler(triangulation);
2099 dof_handler.distribute_dofs(q1);
2100
2101 DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
2102 DoFTools::make_sparsity_pattern(dof_handler, dsp);
2103 dsp.compress();
2104
2105 SparsityPattern sparsity_pattern;
2106 sparsity_pattern.copy_from(dsp);
2107 sparsity_pattern.compress();
2108
2109 SparseMatrix<double> S(sparsity_pattern);
2110
2111 QGauss<dim> quadrature(4);
2112
2114 StaticMappingQ1<dim>::mapping, dof_handler, quadrature, S, coefficient);
2115
2116 // set up the boundary values for the laplace problem
2117 std::array<AffineConstraints<double>, dim> constraints;
2118 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
2119 new_points.end();
2120
2121 // fill these maps using the data given by new_points
2122 for (const auto &cell : dof_handler.active_cell_iterators())
2123 {
2124 // loop over all vertices of the cell and see if it is listed in the map
2125 // given as first argument of the function
2126 for (const unsigned int vertex_no : cell->vertex_indices())
2127 {
2128 const unsigned int vertex_index = cell->vertex_index(vertex_no);
2129 const Point<dim> & vertex_point = cell->vertex(vertex_no);
2130
2131 const typename std::map<unsigned int, Point<dim>>::const_iterator
2132 map_iter = new_points.find(vertex_index);
2133
2134 if (map_iter != map_end)
2135 for (unsigned int i = 0; i < dim; ++i)
2136 {
2137 constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
2138 constraints[i].set_inhomogeneity(
2139 cell->vertex_dof_index(vertex_no, 0),
2140 (solve_for_absolute_positions ?
2141 map_iter->second(i) :
2142 map_iter->second(i) - vertex_point[i]));
2143 }
2144 }
2145 }
2146
2147 for (unsigned int i = 0; i < dim; ++i)
2148 constraints[i].close();
2149
2150 // solve the dim problems with different right hand sides.
2151 Vector<double> us[dim];
2152 for (unsigned int i = 0; i < dim; ++i)
2153 us[i].reinit(dof_handler.n_dofs());
2154
2155 // solve linear systems in parallel
2157 for (unsigned int i = 0; i < dim; ++i)
2158 tasks +=
2159 Threads::new_task(&internal::laplace_solve, S, constraints[i], us[i]);
2160 tasks.join_all();
2161
2162 // change the coordinates of the points of the triangulation
2163 // according to the computed values
2164 std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2165 for (const auto &cell : dof_handler.active_cell_iterators())
2166 for (const unsigned int vertex_no : cell->vertex_indices())
2167 if (vertex_touched[cell->vertex_index(vertex_no)] == false)
2168 {
2169 Point<dim> &v = cell->vertex(vertex_no);
2170
2171 const types::global_dof_index dof_index =
2172 cell->vertex_dof_index(vertex_no, 0);
2173 for (unsigned int i = 0; i < dim; ++i)
2174 if (solve_for_absolute_positions)
2175 v(i) = us[i](dof_index);
2176 else
2177 v(i) += us[i](dof_index);
2178
2179 vertex_touched[cell->vertex_index(vertex_no)] = true;
2180 }
2181 }
2182
2183 template <int dim, int spacedim>
2184 std::map<unsigned int, Point<spacedim>>
2186 {
2187 std::map<unsigned int, Point<spacedim>> vertex_map;
2189 cell = tria.begin_active(),
2190 endc = tria.end();
2191 for (; cell != endc; ++cell)
2192 {
2193 for (unsigned int i : cell->face_indices())
2194 {
2195 const typename Triangulation<dim, spacedim>::face_iterator &face =
2196 cell->face(i);
2197 if (face->at_boundary())
2198 {
2199 for (unsigned j = 0; j < face->n_vertices(); ++j)
2200 {
2201 const Point<spacedim> &vertex = face->vertex(j);
2202 const unsigned int vertex_index = face->vertex_index(j);
2203 vertex_map[vertex_index] = vertex;
2204 }
2205 }
2206 }
2207 }
2208 return vertex_map;
2209 }
2210
2215 template <int dim, int spacedim>
2216 void
2217 distort_random(const double factor,
2219 const bool keep_boundary,
2220 const unsigned int seed)
2221 {
2222 // if spacedim>dim we need to make sure that we perturb
2223 // points but keep them on
2224 // the manifold. however, this isn't implemented right now
2225 Assert(spacedim == dim, ExcNotImplemented());
2226
2227
2228 // find the smallest length of the
2229 // lines adjacent to the
2230 // vertex. take the initial value
2231 // to be larger than anything that
2232 // might be found: the diameter of
2233 // the triangulation, here
2234 // estimated by adding up the
2235 // diameters of the coarse grid
2236 // cells.
2237 double almost_infinite_length = 0;
2239 triangulation.begin(0);
2240 cell != triangulation.end(0);
2241 ++cell)
2242 almost_infinite_length += cell->diameter();
2243
2244 std::vector<double> minimal_length(triangulation.n_vertices(),
2245 almost_infinite_length);
2246
2247 // also note if a vertex is at the boundary
2248 std::vector<bool> at_boundary(keep_boundary ? triangulation.n_vertices() :
2249 0,
2250 false);
2251 // for parallel::shared::Triangulation we need to work on all vertices,
2252 // not just the ones related to locally owned cells;
2253 const bool is_parallel_shared =
2255 &triangulation) != nullptr);
2256 for (const auto &cell : triangulation.active_cell_iterators())
2257 if (is_parallel_shared || cell->is_locally_owned())
2258 {
2259 if (dim > 1)
2260 {
2261 for (unsigned int i = 0; i < cell->n_lines(); ++i)
2262 {
2264 line = cell->line(i);
2265
2266 if (keep_boundary && line->at_boundary())
2267 {
2268 at_boundary[line->vertex_index(0)] = true;
2269 at_boundary[line->vertex_index(1)] = true;
2270 }
2271
2272 minimal_length[line->vertex_index(0)] =
2273 std::min(line->diameter(),
2274 minimal_length[line->vertex_index(0)]);
2275 minimal_length[line->vertex_index(1)] =
2276 std::min(line->diameter(),
2277 minimal_length[line->vertex_index(1)]);
2278 }
2279 }
2280 else // dim==1
2281 {
2282 if (keep_boundary)
2283 for (unsigned int vertex = 0; vertex < 2; ++vertex)
2284 if (cell->at_boundary(vertex) == true)
2285 at_boundary[cell->vertex_index(vertex)] = true;
2286
2287 minimal_length[cell->vertex_index(0)] =
2288 std::min(cell->diameter(),
2289 minimal_length[cell->vertex_index(0)]);
2290 minimal_length[cell->vertex_index(1)] =
2291 std::min(cell->diameter(),
2292 minimal_length[cell->vertex_index(1)]);
2293 }
2294 }
2295
2296 // create a random number generator for the interval [-1,1]
2297 boost::random::mt19937 rng(seed);
2298 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
2299
2300 // If the triangulation is distributed, we need to
2301 // exchange the moved vertices across mpi processes
2302 if (auto distributed_triangulation =
2304 &triangulation))
2305 {
2306 const std::vector<bool> locally_owned_vertices =
2308 std::vector<bool> vertex_moved(triangulation.n_vertices(), false);
2309
2310 // Next move vertices on locally owned cells
2311 for (const auto &cell : triangulation.active_cell_iterators())
2312 if (cell->is_locally_owned())
2313 {
2314 for (const unsigned int vertex_no : cell->vertex_indices())
2315 {
2316 const unsigned global_vertex_no =
2317 cell->vertex_index(vertex_no);
2318
2319 // ignore this vertex if we shall keep the boundary and
2320 // this vertex *is* at the boundary, if it is already moved
2321 // or if another process moves this vertex
2322 if ((keep_boundary && at_boundary[global_vertex_no]) ||
2323 vertex_moved[global_vertex_no] ||
2324 !locally_owned_vertices[global_vertex_no])
2325 continue;
2326
2327 // first compute a random shift vector
2328 Point<spacedim> shift_vector;
2329 for (unsigned int d = 0; d < spacedim; ++d)
2330 shift_vector(d) = uniform_distribution(rng);
2331
2332 shift_vector *= factor * minimal_length[global_vertex_no] /
2333 std::sqrt(shift_vector.square());
2334
2335 // finally move the vertex
2336 cell->vertex(vertex_no) += shift_vector;
2337 vertex_moved[global_vertex_no] = true;
2338 }
2339 }
2340
2341 distributed_triangulation->communicate_locally_moved_vertices(
2342 locally_owned_vertices);
2343 }
2344 else
2345 // if this is a sequential triangulation, we could in principle
2346 // use the algorithm above, but we'll use an algorithm that we used
2347 // before the parallel::distributed::Triangulation was introduced
2348 // in order to preserve backward compatibility
2349 {
2350 // loop over all vertices and compute their new locations
2351 const unsigned int n_vertices = triangulation.n_vertices();
2352 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
2353 const std::vector<Point<spacedim>> &old_vertex_locations =
2354 triangulation.get_vertices();
2355
2356 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2357 {
2358 // ignore this vertex if we will keep the boundary and
2359 // this vertex *is* at the boundary
2360 if (keep_boundary && at_boundary[vertex])
2361 new_vertex_locations[vertex] = old_vertex_locations[vertex];
2362 else
2363 {
2364 // compute a random shift vector
2365 Point<spacedim> shift_vector;
2366 for (unsigned int d = 0; d < spacedim; ++d)
2367 shift_vector(d) = uniform_distribution(rng);
2368
2369 shift_vector *= factor * minimal_length[vertex] /
2370 std::sqrt(shift_vector.square());
2371
2372 // record new vertex location
2373 new_vertex_locations[vertex] =
2374 old_vertex_locations[vertex] + shift_vector;
2375 }
2376 }
2377
2378 // now do the actual move of the vertices
2379 for (const auto &cell : triangulation.active_cell_iterators())
2380 for (const unsigned int vertex_no : cell->vertex_indices())
2381 cell->vertex(vertex_no) =
2382 new_vertex_locations[cell->vertex_index(vertex_no)];
2383 }
2384
2385 // Correct hanging nodes if necessary
2386 if (dim >= 2)
2387 {
2388 // We do the same as in GridTools::transform
2389 //
2390 // exclude hanging nodes at the boundaries of artificial cells:
2391 // these may belong to ghost cells for which we know the exact
2392 // location of vertices, whereas the artificial cell may or may
2393 // not be further refined, and so we cannot know whether
2394 // the location of the hanging node is correct or not
2396 cell = triangulation.begin_active(),
2397 endc = triangulation.end();
2398 for (; cell != endc; ++cell)
2399 if (!cell->is_artificial())
2400 for (const unsigned int face : cell->face_indices())
2401 if (cell->face(face)->has_children() &&
2402 !cell->face(face)->at_boundary())
2403 {
2404 // this face has hanging nodes
2405 if (dim == 2)
2406 cell->face(face)->child(0)->vertex(1) =
2407 (cell->face(face)->vertex(0) +
2408 cell->face(face)->vertex(1)) /
2409 2;
2410 else if (dim == 3)
2411 {
2412 cell->face(face)->child(0)->vertex(1) =
2413 .5 * (cell->face(face)->vertex(0) +
2414 cell->face(face)->vertex(1));
2415 cell->face(face)->child(0)->vertex(2) =
2416 .5 * (cell->face(face)->vertex(0) +
2417 cell->face(face)->vertex(2));
2418 cell->face(face)->child(1)->vertex(3) =
2419 .5 * (cell->face(face)->vertex(1) +
2420 cell->face(face)->vertex(3));
2421 cell->face(face)->child(2)->vertex(3) =
2422 .5 * (cell->face(face)->vertex(2) +
2423 cell->face(face)->vertex(3));
2424
2425 // center of the face
2426 cell->face(face)->child(0)->vertex(3) =
2427 .25 * (cell->face(face)->vertex(0) +
2428 cell->face(face)->vertex(1) +
2429 cell->face(face)->vertex(2) +
2430 cell->face(face)->vertex(3));
2431 }
2432 }
2433 }
2434 }
2435
2436
2437
2438 template <int dim, template <int, int> class MeshType, int spacedim>
2439 unsigned int
2440 find_closest_vertex(const MeshType<dim, spacedim> &mesh,
2441 const Point<spacedim> & p,
2442 const std::vector<bool> & marked_vertices)
2443 {
2444 // first get the underlying triangulation from the mesh and determine
2445 // vertices and used vertices
2447
2448 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
2449
2450 Assert(tria.get_vertices().size() == marked_vertices.size() ||
2451 marked_vertices.size() == 0,
2452 ExcDimensionMismatch(tria.get_vertices().size(),
2453 marked_vertices.size()));
2454
2455 // marked_vertices is expected to be a subset of used_vertices. Thus,
2456 // comparing the range marked_vertices.begin() to marked_vertices.end() with
2457 // the range used_vertices.begin() to used_vertices.end() the element in the
2458 // second range must be valid if the element in the first range is valid.
2459 Assert(
2460 marked_vertices.size() == 0 ||
2461 std::equal(marked_vertices.begin(),
2462 marked_vertices.end(),
2463 tria.get_used_vertices().begin(),
2464 [](bool p, bool q) { return !p || q; }),
2465 ExcMessage(
2466 "marked_vertices should be a subset of used vertices in the triangulation "
2467 "but marked_vertices contains one or more vertices that are not used vertices!"));
2468
2469 // If marked_indices is empty, consider all used_vertices for finding the
2470 // closest vertex to the point. Otherwise, marked_indices is used.
2471 const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
2472 tria.get_used_vertices() :
2473 marked_vertices;
2474
2475 // At the beginning, the first used vertex is considered to be the closest
2476 // one.
2477 std::vector<bool>::const_iterator first =
2478 std::find(vertices_to_use.begin(), vertices_to_use.end(), true);
2479
2480 // Assert that at least one vertex is actually used
2481 Assert(first != vertices_to_use.end(), ExcInternalError());
2482
2483 unsigned int best_vertex = std::distance(vertices_to_use.begin(), first);
2484 double best_dist = (p - vertices[best_vertex]).norm_square();
2485
2486 // For all remaining vertices, test
2487 // whether they are any closer
2488 for (unsigned int j = best_vertex + 1; j < vertices.size(); j++)
2489 if (vertices_to_use[j])
2490 {
2491 const double dist = (p - vertices[j]).norm_square();
2492 if (dist < best_dist)
2493 {
2494 best_vertex = j;
2495 best_dist = dist;
2496 }
2497 }
2498
2499 return best_vertex;
2500 }
2501
2502
2503
2504 template <int dim, template <int, int> class MeshType, int spacedim>
2505 unsigned int
2507 const MeshType<dim, spacedim> &mesh,
2508 const Point<spacedim> & p,
2509 const std::vector<bool> & marked_vertices)
2510 {
2511 // Take a shortcut in the simple case.
2512 if (mapping.preserves_vertex_locations() == true)
2513 return find_closest_vertex(mesh, p, marked_vertices);
2514
2515 // first get the underlying triangulation from the mesh and determine
2516 // vertices and used vertices
2518
2519 auto vertices = extract_used_vertices(tria, mapping);
2520
2521 Assert(tria.get_vertices().size() == marked_vertices.size() ||
2522 marked_vertices.size() == 0,
2523 ExcDimensionMismatch(tria.get_vertices().size(),
2524 marked_vertices.size()));
2525
2526 // marked_vertices is expected to be a subset of used_vertices. Thus,
2527 // comparing the range marked_vertices.begin() to marked_vertices.end()
2528 // with the range used_vertices.begin() to used_vertices.end() the element
2529 // in the second range must be valid if the element in the first range is
2530 // valid.
2531 Assert(
2532 marked_vertices.size() == 0 ||
2533 std::equal(marked_vertices.begin(),
2534 marked_vertices.end(),
2535 tria.get_used_vertices().begin(),
2536 [](bool p, bool q) { return !p || q; }),
2537 ExcMessage(
2538 "marked_vertices should be a subset of used vertices in the triangulation "
2539 "but marked_vertices contains one or more vertices that are not used vertices!"));
2540
2541 // Remove from the map unwanted elements.
2542 if (marked_vertices.size() != 0)
2543 for (auto it = vertices.begin(); it != vertices.end();)
2544 {
2545 if (marked_vertices[it->first] == false)
2546 {
2547 it = vertices.erase(it);
2548 }
2549 else
2550 {
2551 ++it;
2552 }
2553 }
2554
2555 return find_closest_vertex(vertices, p);
2556 }
2557
2558
2559
2560 template <int dim, template <int, int> class MeshType, int spacedim>
2561#ifndef _MSC_VER
2562 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
2563#else
2564 std::vector<
2565 typename ::internal::
2566 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2567#endif
2568 find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &mesh,
2569 const unsigned int vertex)
2570 {
2571 // make sure that the given vertex is
2572 // an active vertex of the underlying
2573 // triangulation
2574 AssertIndexRange(vertex, mesh.get_triangulation().n_vertices());
2575 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
2576 ExcVertexNotUsed(vertex));
2577
2578 // use a set instead of a vector
2579 // to ensure that cells are inserted only
2580 // once
2581 std::set<typename ::internal::
2582 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2584
2585 // go through all active cells and look if the vertex is part of that cell
2586 //
2587 // in 1d, this is all we need to care about. in 2d/3d we also need to worry
2588 // that the vertex might be a hanging node on a face or edge of a cell; in
2589 // this case, we would want to add those cells as well on whose faces the
2590 // vertex is located but for which it is not a vertex itself.
2591 //
2592 // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging
2593 // node can only be in the middle of a face and we can query the neighboring
2594 // cell from the current cell. on the other hand, in 3d a hanging node
2595 // vertex can also be on an edge but there can be many other cells on
2596 // this edge and we can not access them from the cell we are currently
2597 // on.
2598 //
2599 // so, in the 3d case, if we run the algorithm as in 2d, we catch all
2600 // those cells for which the vertex we seek is on a *subface*, but we
2601 // miss the case of cells for which the vertex we seek is on a
2602 // sub-edge for which there is no corresponding sub-face (because the
2603 // immediate neighbor behind this face is not refined), see for example
2604 // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we
2605 // haven't yet found the vertex for the current cell we also need to
2606 // look at the mid-points of edges
2607 //
2608 // as a final note, deciding whether a neighbor is actually coarser is
2609 // simple in the case of isotropic refinement (we just need to look at
2610 // the level of the current and the neighboring cell). however, this
2611 // isn't so simple if we have used anisotropic refinement since then
2612 // the level of a cell is not indicative of whether it is coarser or
2613 // not than the current cell. ultimately, we want to add all cells on
2614 // which the vertex is, independent of whether they are coarser or
2615 // finer and so in the 2d case below we simply add *any* *active* neighbor.
2616 // in the worst case, we add cells multiple times to the adjacent_cells
2617 // list, but std::set throws out those cells already entered
2618 for (const auto &cell : mesh.active_cell_iterators())
2619 {
2620 for (const unsigned int v : cell->vertex_indices())
2621 if (cell->vertex_index(v) == vertex)
2622 {
2623 // OK, we found a cell that contains
2624 // the given vertex. We add it
2625 // to the list.
2626 adjacent_cells.insert(cell);
2627
2628 // as explained above, in 2+d we need to check whether
2629 // this vertex is on a face behind which there is a
2630 // (possibly) coarser neighbor. if this is the case,
2631 // then we need to also add this neighbor
2632 if (dim >= 2)
2633 for (const auto face :
2634 cell->reference_cell().faces_for_given_vertex(v))
2635 if (!cell->at_boundary(face) &&
2636 cell->neighbor(face)->is_active())
2637 {
2638 // there is a (possibly) coarser cell behind a
2639 // face to which the vertex belongs. the
2640 // vertex we are looking at is then either a
2641 // vertex of that coarser neighbor, or it is a
2642 // hanging node on one of the faces of that
2643 // cell. in either case, it is adjacent to the
2644 // vertex, so add it to the list as well (if
2645 // the cell was already in the list then the
2646 // std::set makes sure that we get it only
2647 // once)
2648 adjacent_cells.insert(cell->neighbor(face));
2649 }
2650
2651 // in any case, we have found a cell, so go to the next cell
2652 goto next_cell;
2653 }
2654
2655 // in 3d also loop over the edges
2656 if (dim >= 3)
2657 {
2658 for (unsigned int e = 0; e < cell->n_lines(); ++e)
2659 if (cell->line(e)->has_children())
2660 // the only place where this vertex could have been
2661 // hiding is on the mid-edge point of the edge we
2662 // are looking at
2663 if (cell->line(e)->child(0)->vertex_index(1) == vertex)
2664 {
2665 adjacent_cells.insert(cell);
2666
2667 // jump out of this tangle of nested loops
2668 goto next_cell;
2669 }
2670 }
2671
2672 // in more than 3d we would probably have to do the same as
2673 // above also for even lower-dimensional objects
2674 Assert(dim <= 3, ExcNotImplemented());
2675
2676 // move on to the next cell if we have found the
2677 // vertex on the current one
2678 next_cell:;
2679 }
2680
2681 // if this was an active vertex then there needs to have been
2682 // at least one cell to which it is adjacent!
2683 Assert(adjacent_cells.size() > 0, ExcInternalError());
2684
2685 // return the result as a vector, rather than the set we built above
2686 return std::vector<
2687 typename ::internal::
2688 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
2689 adjacent_cells.begin(), adjacent_cells.end());
2690 }
2691
2692
2693
2694 template <int dim, int spacedim>
2695 std::vector<std::vector<Tensor<1, spacedim>>>
2697 const Triangulation<dim, spacedim> &mesh,
2698 const std::vector<
2700 &vertex_to_cells)
2701 {
2702 const std::vector<Point<spacedim>> &vertices = mesh.get_vertices();
2703 const unsigned int n_vertices = vertex_to_cells.size();
2704
2705 AssertDimension(vertices.size(), n_vertices);
2706
2707
2708 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
2709 n_vertices);
2710 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2711 if (mesh.vertex_used(vertex))
2712 {
2713 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
2714 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
2715
2716 typename std::set<typename Triangulation<dim, spacedim>::
2717 active_cell_iterator>::iterator it =
2718 vertex_to_cells[vertex].begin();
2719 for (unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
2720 {
2721 vertex_to_cell_centers[vertex][cell] =
2722 (*it)->center() - vertices[vertex];
2723 vertex_to_cell_centers[vertex][cell] /=
2724 vertex_to_cell_centers[vertex][cell].norm();
2725 }
2726 }
2727 return vertex_to_cell_centers;
2728 }
2729
2730
2731 namespace internal
2732 {
2733 template <int spacedim>
2734 bool
2736 const unsigned int a,
2737 const unsigned int b,
2738 const Tensor<1, spacedim> & point_direction,
2739 const std::vector<Tensor<1, spacedim>> &center_directions)
2740 {
2741 const double scalar_product_a = center_directions[a] * point_direction;
2742 const double scalar_product_b = center_directions[b] * point_direction;
2743
2744 // The function is supposed to return if a is before b. We are looking
2745 // for the alignment of point direction and center direction, therefore
2746 // return if the scalar product of a is larger.
2747 return (scalar_product_a > scalar_product_b);
2748 }
2749 } // namespace internal
2750
2751 template <int dim, template <int, int> class MeshType, int spacedim>
2752#ifndef _MSC_VER
2753 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
2754#else
2755 std::pair<typename ::internal::
2756 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
2757 Point<dim>>
2758#endif
2760 const Mapping<dim, spacedim> & mapping,
2761 const MeshType<dim, spacedim> &mesh,
2762 const Point<spacedim> & p,
2763 const std::vector<
2764 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
2765 & vertex_to_cells,
2766 const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
2767 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
2768 const std::vector<bool> & marked_vertices,
2769 const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree,
2770 const double tolerance)
2771 {
2772 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2773 Point<dim>>
2774 cell_and_position;
2775 // To handle points at the border we keep track of points which are close to
2776 // the unit cell:
2777 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2778 Point<dim>>
2779 cell_and_position_approx;
2780
2781 bool found_cell = false;
2782 bool approx_cell = false;
2783
2784 unsigned int closest_vertex_index = 0;
2785 Tensor<1, spacedim> vertex_to_point;
2786 auto current_cell = cell_hint;
2787
2788 while (found_cell == false)
2789 {
2790 // First look at the vertices of the cell cell_hint. If it's an
2791 // invalid cell, then query for the closest global vertex
2792 if (current_cell.state() == IteratorState::valid)
2793 {
2794 const auto cell_vertices = mapping.get_vertices(current_cell);
2795 const unsigned int closest_vertex =
2796 find_closest_vertex_of_cell<dim, spacedim>(current_cell,
2797 p,
2798 mapping);
2799 vertex_to_point = p - cell_vertices[closest_vertex];
2800 closest_vertex_index = current_cell->vertex_index(closest_vertex);
2801 }
2802 else
2803 {
2804 if (!used_vertices_rtree.empty())
2805 {
2806 // If we have an rtree at our disposal, use it.
2807 using ValueType = std::pair<Point<spacedim>, unsigned int>;
2808 std::function<bool(const ValueType &)> marked;
2809 if (marked_vertices.size() == mesh.n_vertices())
2810 marked = [&marked_vertices](const ValueType &value) -> bool {
2811 return marked_vertices[value.second];
2812 };
2813 else
2814 marked = [](const ValueType &) -> bool { return true; };
2815
2816 std::vector<std::pair<Point<spacedim>, unsigned int>> res;
2817 used_vertices_rtree.query(
2818 boost::geometry::index::nearest(p, 1) &&
2819 boost::geometry::index::satisfies(marked),
2820 std::back_inserter(res));
2821
2822 // We should have one and only one result
2823 AssertDimension(res.size(), 1);
2824 closest_vertex_index = res[0].second;
2825 }
2826 else
2827 {
2828 closest_vertex_index = GridTools::find_closest_vertex(
2829 mapping, mesh, p, marked_vertices);
2830 }
2831 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
2832 }
2833
2834 const double vertex_point_norm = vertex_to_point.norm();
2835 if (vertex_point_norm > 0)
2836 vertex_to_point /= vertex_point_norm;
2837
2838 const unsigned int n_neighbor_cells =
2839 vertex_to_cells[closest_vertex_index].size();
2840
2841 // Create a corresponding map of vectors from vertex to cell center
2842 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
2843
2844 for (unsigned int i = 0; i < n_neighbor_cells; ++i)
2845 neighbor_permutation[i] = i;
2846
2847 auto comp = [&](const unsigned int a, const unsigned int b) -> bool {
2848 return internal::compare_point_association<spacedim>(
2849 a,
2850 b,
2851 vertex_to_point,
2852 vertex_to_cell_centers[closest_vertex_index]);
2853 };
2854
2855 std::sort(neighbor_permutation.begin(),
2856 neighbor_permutation.end(),
2857 comp);
2858 // It is possible the vertex is close
2859 // to an edge, thus we add a tolerance
2860 // to keep also the "best" cell
2861 double best_distance = tolerance;
2862
2863 // Search all of the cells adjacent to the closest vertex of the cell
2864 // hint Most likely we will find the point in them.
2865 for (unsigned int i = 0; i < n_neighbor_cells; ++i)
2866 {
2867 try
2868 {
2869 auto cell = vertex_to_cells[closest_vertex_index].begin();
2870 std::advance(cell, neighbor_permutation[i]);
2871
2872 if (!(*cell)->is_artificial())
2873 {
2874 const Point<dim> p_unit =
2875 mapping.transform_real_to_unit_cell(*cell, p);
2877 tolerance))
2878 {
2879 cell_and_position.first = *cell;
2880 cell_and_position.second = p_unit;
2881 found_cell = true;
2882 approx_cell = false;
2883 break;
2884 }
2885 // The point is not inside this cell: checking how far
2886 // outside it is and whether we want to use this cell as a
2887 // backup if we can't find a cell within which the point
2888 // lies.
2889 const double dist =
2891 if (dist < best_distance)
2892 {
2893 best_distance = dist;
2894 cell_and_position_approx.first = *cell;
2895 cell_and_position_approx.second = p_unit;
2896 approx_cell = true;
2897 }
2898 }
2899 }
2900 catch (typename Mapping<dim>::ExcTransformationFailed &)
2901 {}
2902 }
2903
2904 if (found_cell == true)
2905 return cell_and_position;
2906 else if (approx_cell == true)
2907 return cell_and_position_approx;
2908
2909 // The first time around, we check for vertices in the hint_cell. If
2910 // that does not work, we set the cell iterator to an invalid one, and
2911 // look for a global vertex close to the point. If that does not work,
2912 // we are in trouble, and just throw an exception.
2913 //
2914 // If we got here, then we did not find the point. If the
2915 // current_cell.state() here is not IteratorState::valid, it means that
2916 // the user did not provide a hint_cell, and at the beginning of the
2917 // while loop we performed an actual global search on the mesh
2918 // vertices. Not finding the point then means the point is outside the
2919 // domain, or that we've had problems with the algorithm above. Try as a
2920 // last resort the other (simpler) algorithm.
2921 if (current_cell.state() != IteratorState::valid)
2923 mapping, mesh, p, marked_vertices, tolerance);
2924
2925 current_cell = typename MeshType<dim, spacedim>::active_cell_iterator();
2926 }
2927 return cell_and_position;
2928 }
2929
2930
2931
2932 template <int dim, int spacedim>
2933 unsigned int
2936 const Point<spacedim> & position,
2937 const Mapping<dim, spacedim> & mapping)
2938 {
2939 const auto vertices = mapping.get_vertices(cell);
2940 double minimum_distance = position.distance_square(vertices[0]);
2941 unsigned int closest_vertex = 0;
2942
2943 for (unsigned int v = 1; v < cell->n_vertices(); ++v)
2944 {
2945 const double vertex_distance = position.distance_square(vertices[v]);
2946 if (vertex_distance < minimum_distance)
2947 {
2948 closest_vertex = v;
2949 minimum_distance = vertex_distance;
2950 }
2951 }
2952 return closest_vertex;
2953 }
2954
2955
2956
2957 namespace internal
2958 {
2959 namespace BoundingBoxPredicate
2960 {
2961 template <class MeshType>
2962 std::tuple<BoundingBox<MeshType::space_dimension>, bool>
2964 const typename MeshType::cell_iterator &parent_cell,
2965 const std::function<
2966 bool(const typename MeshType::active_cell_iterator &)> &predicate)
2967 {
2968 bool has_predicate =
2969 false; // Start assuming there's no cells with predicate inside
2970 std::vector<typename MeshType::active_cell_iterator> active_cells;
2971 if (parent_cell->is_active())
2972 active_cells = {parent_cell};
2973 else
2974 // Finding all active cells descendants of the current one (or the
2975 // current one if it is active)
2976 active_cells = get_active_child_cells<MeshType>(parent_cell);
2977
2978 const unsigned int spacedim = MeshType::space_dimension;
2979
2980 // Looking for the first active cell which has the property predicate
2981 unsigned int i = 0;
2982 while (i < active_cells.size() && !predicate(active_cells[i]))
2983 ++i;
2984
2985 // No active cells or no active cells with property
2986 if (active_cells.size() == 0 || i == active_cells.size())
2987 {
2989 return std::make_tuple(bbox, has_predicate);
2990 }
2991
2992 // The two boundary points defining the boundary box
2993 Point<spacedim> maxp = active_cells[i]->vertex(0);
2994 Point<spacedim> minp = active_cells[i]->vertex(0);
2995
2996 for (; i < active_cells.size(); ++i)
2997 if (predicate(active_cells[i]))
2998 for (const unsigned int v : active_cells[i]->vertex_indices())
2999 for (unsigned int d = 0; d < spacedim; ++d)
3000 {
3001 minp[d] = std::min(minp[d], active_cells[i]->vertex(v)[d]);
3002 maxp[d] = std::max(maxp[d], active_cells[i]->vertex(v)[d]);
3003 }
3004
3005 has_predicate = true;
3006 BoundingBox<spacedim> bbox(std::make_pair(minp, maxp));
3007 return std::make_tuple(bbox, has_predicate);
3008 }
3009 } // namespace BoundingBoxPredicate
3010 } // namespace internal
3011
3012
3013
3014 template <class MeshType>
3015 std::vector<BoundingBox<MeshType::space_dimension>>
3017 const MeshType &mesh,
3018 const std::function<bool(const typename MeshType::active_cell_iterator &)>
3019 & predicate,
3020 const unsigned int refinement_level,
3021 const bool allow_merge,
3022 const unsigned int max_boxes)
3023 {
3024 // Algorithm brief description: begin with creating bounding boxes of all
3025 // cells at refinement_level (and coarser levels if there are active cells)
3026 // which have the predicate property. These are then merged
3027
3028 Assert(
3029 refinement_level <= mesh.n_levels(),
3030 ExcMessage(
3031 "Error: refinement level is higher then total levels in the triangulation!"));
3032
3033 const unsigned int spacedim = MeshType::space_dimension;
3034 std::vector<BoundingBox<spacedim>> bounding_boxes;
3035
3036 // Creating a bounding box for all active cell on coarser level
3037
3038 for (unsigned int i = 0; i < refinement_level; ++i)
3039 for (const typename MeshType::cell_iterator &cell :
3040 mesh.active_cell_iterators_on_level(i))
3041 {
3042 bool has_predicate = false;
3044 std::tie(bbox, has_predicate) =
3046 MeshType>(cell, predicate);
3047 if (has_predicate)
3048 bounding_boxes.push_back(bbox);
3049 }
3050
3051 // Creating a Bounding Box for all cells on the chosen refinement_level
3052 for (const typename MeshType::cell_iterator &cell :
3053 mesh.cell_iterators_on_level(refinement_level))
3054 {
3055 bool has_predicate = false;
3057 std::tie(bbox, has_predicate) =
3059 MeshType>(cell, predicate);
3060 if (has_predicate)
3061 bounding_boxes.push_back(bbox);
3062 }
3063
3064 if (!allow_merge)
3065 // If merging is not requested return the created bounding_boxes
3066 return bounding_boxes;
3067 else
3068 {
3069 // Merging part of the algorithm
3070 // Part 1: merging neighbors
3071 // This array stores the indices of arrays we have already merged
3072 std::vector<unsigned int> merged_boxes_idx;
3073 bool found_neighbors = true;
3074
3075 // We merge only neighbors which can be expressed by a single bounding
3076 // box e.g. in 1d [0,1] and [1,2] can be described with [0,2] without
3077 // losing anything
3078 while (found_neighbors)
3079 {
3080 found_neighbors = false;
3081 for (unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
3082 {
3083 if (std::find(merged_boxes_idx.begin(),
3084 merged_boxes_idx.end(),
3085 i) == merged_boxes_idx.end())
3086 for (unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
3087 if (std::find(merged_boxes_idx.begin(),
3088 merged_boxes_idx.end(),
3089 j) == merged_boxes_idx.end() &&
3090 bounding_boxes[i].get_neighbor_type(
3091 bounding_boxes[j]) ==
3093 {
3094 bounding_boxes[i].merge_with(bounding_boxes[j]);
3095 merged_boxes_idx.push_back(j);
3096 found_neighbors = true;
3097 }
3098 }
3099 }
3100
3101 // Copying the merged boxes into merged_b_boxes
3102 std::vector<BoundingBox<spacedim>> merged_b_boxes;
3103 for (unsigned int i = 0; i < bounding_boxes.size(); ++i)
3104 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
3105 merged_boxes_idx.end())
3106 merged_b_boxes.push_back(bounding_boxes[i]);
3107
3108 // Part 2: if there are too many bounding boxes, merging smaller boxes
3109 // This has sense only in dimension 2 or greater, since in dimension 1,
3110 // neighboring intervals can always be merged without problems
3111 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
3112 {
3113 std::vector<double> volumes;
3114 for (unsigned int i = 0; i < merged_b_boxes.size(); ++i)
3115 volumes.push_back(merged_b_boxes[i].volume());
3116
3117 while (merged_b_boxes.size() > max_boxes)
3118 {
3119 unsigned int min_idx =
3120 std::min_element(volumes.begin(), volumes.end()) -
3121 volumes.begin();
3122 volumes.erase(volumes.begin() + min_idx);
3123 // Finding a neighbor
3124 bool not_removed = true;
3125 for (unsigned int i = 0;
3126 i < merged_b_boxes.size() && not_removed;
3127 ++i)
3128 // We merge boxes if we have "attached" or "mergeable"
3129 // neighbors, even though mergeable should be dealt with in
3130 // Part 1
3131 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
3132 merged_b_boxes[min_idx]) ==
3134 merged_b_boxes[i].get_neighbor_type(
3135 merged_b_boxes[min_idx]) ==
3137 {
3138 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
3139 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
3140 not_removed = false;
3141 }
3142 Assert(!not_removed,
3143 ExcMessage("Error: couldn't merge bounding boxes!"));
3144 }
3145 }
3146 Assert(merged_b_boxes.size() <= max_boxes,
3147 ExcMessage(
3148 "Error: couldn't reach target number of bounding boxes!"));
3149 return merged_b_boxes;
3150 }
3151 }
3152
3153
3154
3155 template <int spacedim>
3156#ifndef DOXYGEN
3157 std::tuple<std::vector<std::vector<unsigned int>>,
3158 std::map<unsigned int, unsigned int>,
3159 std::map<unsigned int, std::vector<unsigned int>>>
3160#else
3161 return_type
3162#endif
3164 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3165 const std::vector<Point<spacedim>> & points)
3166 {
3167 unsigned int n_procs = global_bboxes.size();
3168 std::vector<std::vector<unsigned int>> point_owners(n_procs);
3169 std::map<unsigned int, unsigned int> map_owners_found;
3170 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3171
3172 unsigned int n_points = points.size();
3173 for (unsigned int pt = 0; pt < n_points; ++pt)
3174 {
3175 // Keep track of how many processes we guess to own the point
3176 std::vector<unsigned int> owners_found;
3177 // Check in which other processes the point might be
3178 for (unsigned int rk = 0; rk < n_procs; ++rk)
3179 {
3180 for (const BoundingBox<spacedim> &bbox : global_bboxes[rk])
3181 if (bbox.point_inside(points[pt]))
3182 {
3183 point_owners[rk].emplace_back(pt);
3184 owners_found.emplace_back(rk);
3185 break; // We can check now the next process
3186 }
3187 }
3188 Assert(owners_found.size() > 0,
3189 ExcMessage("No owners found for the point " +
3190 std::to_string(pt)));
3191 if (owners_found.size() == 1)
3192 map_owners_found[pt] = owners_found[0];
3193 else
3194 // Multiple owners
3195 map_owners_guessed[pt] = owners_found;
3196 }
3197
3198 return std::make_tuple(std::move(point_owners),
3199 std::move(map_owners_found),
3200 std::move(map_owners_guessed));
3201 }
3202
3203 template <int spacedim>
3204#ifndef DOXYGEN
3205 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
3206 std::map<unsigned int, unsigned int>,
3207 std::map<unsigned int, std::vector<unsigned int>>>
3208#else
3209 return_type
3210#endif
3212 const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
3213 const std::vector<Point<spacedim>> & points)
3214 {
3215 std::map<unsigned int, std::vector<unsigned int>> point_owners;
3216 std::map<unsigned int, unsigned int> map_owners_found;
3217 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3218 std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> search_result;
3219
3220 unsigned int n_points = points.size();
3221 for (unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
3222 {
3223 search_result.clear(); // clearing last output
3224
3225 // Running tree search
3226 covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
3227 std::back_inserter(search_result));
3228
3229 // Keep track of how many processes we guess to own the point
3230 std::set<unsigned int> owners_found;
3231 // Check in which other processes the point might be
3232 for (const auto &rank_bbox : search_result)
3233 {
3234 // Try to add the owner to the owners found,
3235 // and check if it was already present
3236 const bool pt_inserted = owners_found.insert(pt_n).second;
3237 if (pt_inserted)
3238 point_owners[rank_bbox.second].emplace_back(pt_n);
3239 }
3240 Assert(owners_found.size() > 0,
3241 ExcMessage("No owners found for the point " +
3242 std::to_string(pt_n)));
3243 if (owners_found.size() == 1)
3244 map_owners_found[pt_n] = *owners_found.begin();
3245 else
3246 // Multiple owners
3247 std::copy(owners_found.begin(),
3248 owners_found.end(),
3249 std::back_inserter(map_owners_guessed[pt_n]));
3250 }
3251
3252 return std::make_tuple(std::move(point_owners),
3253 std::move(map_owners_found),
3254 std::move(map_owners_guessed));
3255 }
3256
3257
3258 template <int dim, int spacedim>
3259 std::vector<
3260 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3262 {
3263 std::vector<
3264 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3265 vertex_to_cell_map(triangulation.n_vertices());
3267 cell = triangulation.begin_active(),
3268 endc = triangulation.end();
3269 for (; cell != endc; ++cell)
3270 for (const unsigned int i : cell->vertex_indices())
3271 vertex_to_cell_map[cell->vertex_index(i)].insert(cell);
3272
3273 // Take care of hanging nodes
3274 cell = triangulation.begin_active();
3275 for (; cell != endc; ++cell)
3276 {
3277 for (unsigned int i : cell->face_indices())
3278 {
3279 if ((cell->at_boundary(i) == false) &&
3280 (cell->neighbor(i)->is_active()))
3281 {
3283 adjacent_cell = cell->neighbor(i);
3284 for (unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
3285 vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert(
3286 adjacent_cell);
3287 }
3288 }
3289
3290 // in 3d also loop over the edges
3291 if (dim == 3)
3292 {
3293 for (unsigned int i = 0; i < cell->n_lines(); ++i)
3294 if (cell->line(i)->has_children())
3295 // the only place where this vertex could have been
3296 // hiding is on the mid-edge point of the edge we
3297 // are looking at
3298 vertex_to_cell_map[cell->line(i)->child(0)->vertex_index(1)]
3299 .insert(cell);
3300 }
3301 }
3302
3303 return vertex_to_cell_map;
3304 }
3305
3306
3307
3308 template <int dim, int spacedim>
3309 std::map<unsigned int, types::global_vertex_index>
3312 {
3313 std::map<unsigned int, types::global_vertex_index>
3314 local_to_global_vertex_index;
3315
3316#ifndef DEAL_II_WITH_MPI
3317
3318 // without MPI, this function doesn't make sense because on cannot
3319 // use parallel::distributed::Triangulation in any meaningful
3320 // way
3321 (void)triangulation;
3322 Assert(false,
3323 ExcMessage("This function does not make any sense "
3324 "for parallel::distributed::Triangulation "
3325 "objects if you do not have MPI enabled."));
3326
3327#else
3328
3329 using active_cell_iterator =
3331 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
3333
3334 // Create a local index for the locally "owned" vertices
3335 types::global_vertex_index next_index = 0;
3336 unsigned int max_cellid_size = 0;
3337 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
3338 vertices_added;
3339 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
3340 std::map<types::subdomain_id,
3341 std::vector<std::tuple<types::global_vertex_index,
3343 std::string>>>
3344 vertices_to_send;
3345 active_cell_iterator cell = triangulation.begin_active(),
3346 endc = triangulation.end();
3347 std::set<active_cell_iterator> missing_vert_cells;
3348 std::set<unsigned int> used_vertex_index;
3349 for (; cell != endc; ++cell)
3350 {
3351 if (cell->is_locally_owned())
3352 {
3353 for (const unsigned int i : cell->vertex_indices())
3354 {
3355 types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
3356 typename std::set<active_cell_iterator>::iterator
3357 adjacent_cell = vertex_to_cell[cell->vertex_index(i)].begin(),
3358 end_adj_cell = vertex_to_cell[cell->vertex_index(i)].end();
3359 for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
3360 lowest_subdomain_id =
3361 std::min(lowest_subdomain_id,
3362 (*adjacent_cell)->subdomain_id());
3363
3364 // See if I "own" this vertex
3365 if (lowest_subdomain_id == cell->subdomain_id())
3366 {
3367 // Check that the vertex we are working on a vertex that has
3368 // not be dealt with yet
3369 if (used_vertex_index.find(cell->vertex_index(i)) ==
3370 used_vertex_index.end())
3371 {
3372 // Set the local index
3373 local_to_global_vertex_index[cell->vertex_index(i)] =
3374 next_index++;
3375
3376 // Store the information that will be sent to the
3377 // adjacent cells on other subdomains
3378 adjacent_cell =
3379 vertex_to_cell[cell->vertex_index(i)].begin();
3380 for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
3381 if ((*adjacent_cell)->subdomain_id() !=
3382 cell->subdomain_id())
3383 {
3384 std::pair<types::subdomain_id,
3386 tmp((*adjacent_cell)->subdomain_id(),
3387 cell->vertex_index(i));
3388 if (vertices_added.find(tmp) ==
3389 vertices_added.end())
3390 {
3391 vertices_to_send[(*adjacent_cell)
3392 ->subdomain_id()]
3393 .emplace_back(i,
3394 cell->vertex_index(i),
3395 cell->id().to_string());
3396 if (cell->id().to_string().size() >
3397 max_cellid_size)
3398 max_cellid_size =
3399 cell->id().to_string().size();
3400 vertices_added.insert(tmp);
3401 }
3402 }
3403 used_vertex_index.insert(cell->vertex_index(i));
3404 }
3405 }
3406 else
3407 {
3408 // We don't own the vertex so we will receive its global
3409 // index
3410 vertices_to_recv[lowest_subdomain_id].insert(
3411 cell->vertex_index(i));
3412 missing_vert_cells.insert(cell);
3413 }
3414 }
3415 }
3416
3417 // Some hanging nodes are vertices of ghost cells. They need to be
3418 // received.
3419 if (cell->is_ghost())
3420 {
3421 for (unsigned int i : cell->face_indices())
3422 {
3423 if (cell->at_boundary(i) == false)
3424 {
3425 if (cell->neighbor(i)->is_active())
3426 {
3427 typename Triangulation<dim,
3428 spacedim>::active_cell_iterator
3429 adjacent_cell = cell->neighbor(i);
3430 if ((adjacent_cell->is_locally_owned()))
3431 {
3432 types::subdomain_id adj_subdomain_id =
3433 adjacent_cell->subdomain_id();
3434 if (cell->subdomain_id() < adj_subdomain_id)
3435 for (unsigned int j = 0;
3436 j < cell->face(i)->n_vertices();
3437 ++j)
3438 {
3439 vertices_to_recv[cell->subdomain_id()].insert(
3440 cell->face(i)->vertex_index(j));
3441 missing_vert_cells.insert(cell);
3442 }
3443 }
3444 }
3445 }
3446 }
3447 }
3448 }
3449
3450 // Get the size of the largest CellID string
3451 max_cellid_size =
3452 Utilities::MPI::max(max_cellid_size, triangulation.get_communicator());
3453
3454 // Make indices global by getting the number of vertices owned by each
3455 // processors and shifting the indices accordingly
3457 int ierr = MPI_Exscan(&next_index,
3458 &shift,
3459 1,
3461 MPI_SUM,
3462 triangulation.get_communicator());
3463 AssertThrowMPI(ierr);
3464
3465 std::map<unsigned int, types::global_vertex_index>::iterator
3466 global_index_it = local_to_global_vertex_index.begin(),
3467 global_index_end = local_to_global_vertex_index.end();
3468 for (; global_index_it != global_index_end; ++global_index_it)
3469 global_index_it->second += shift;
3470
3471
3472 const int mpi_tag = Utilities::MPI::internal::Tags::
3474 const int mpi_tag2 = Utilities::MPI::internal::Tags::
3476
3477
3478 // In a first message, send the global ID of the vertices and the local
3479 // positions in the cells. In a second messages, send the cell ID as a
3480 // resize string. This is done in two messages so that types are not mixed
3481
3482 // Send the first message
3483 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
3484 vertices_to_send.size());
3485 std::vector<MPI_Request> first_requests(vertices_to_send.size());
3486 typename std::map<types::subdomain_id,
3487 std::vector<std::tuple<types::global_vertex_index,
3489 std::string>>>::iterator
3490 vert_to_send_it = vertices_to_send.begin(),
3491 vert_to_send_end = vertices_to_send.end();
3492 for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3493 ++vert_to_send_it, ++i)
3494 {
3495 int destination = vert_to_send_it->first;
3496 const unsigned int n_vertices = vert_to_send_it->second.size();
3497 const int buffer_size = 2 * n_vertices;
3498 vertices_send_buffers[i].resize(buffer_size);
3499
3500 // fill the buffer
3501 for (unsigned int j = 0; j < n_vertices; ++j)
3502 {
3503 vertices_send_buffers[i][2 * j] =
3504 std::get<0>(vert_to_send_it->second[j]);
3505 vertices_send_buffers[i][2 * j + 1] =
3506 local_to_global_vertex_index[std::get<1>(
3507 vert_to_send_it->second[j])];
3508 }
3509
3510 // Send the message
3511 ierr = MPI_Isend(vertices_send_buffers[i].data(),
3512 buffer_size,
3514 destination,
3515 mpi_tag,
3516 triangulation.get_communicator(),
3517 &first_requests[i]);
3518 AssertThrowMPI(ierr);
3519 }
3520
3521 // Receive the first message
3522 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
3523 vertices_to_recv.size());
3524 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
3525 vert_to_recv_it = vertices_to_recv.begin(),
3526 vert_to_recv_end = vertices_to_recv.end();
3527 for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3528 ++vert_to_recv_it, ++i)
3529 {
3530 int source = vert_to_recv_it->first;
3531 const unsigned int n_vertices = vert_to_recv_it->second.size();
3532 const int buffer_size = 2 * n_vertices;
3533 vertices_recv_buffers[i].resize(buffer_size);
3534
3535 // Receive the message
3536 ierr = MPI_Recv(vertices_recv_buffers[i].data(),
3537 buffer_size,
3539 source,
3540 mpi_tag,
3541 triangulation.get_communicator(),
3542 MPI_STATUS_IGNORE);
3543 AssertThrowMPI(ierr);
3544 }
3545
3546
3547 // Send second message
3548 std::vector<std::vector<char>> cellids_send_buffers(
3549 vertices_to_send.size());
3550 std::vector<MPI_Request> second_requests(vertices_to_send.size());
3551 vert_to_send_it = vertices_to_send.begin();
3552 for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3553 ++vert_to_send_it, ++i)
3554 {
3555 int destination = vert_to_send_it->first;
3556 const unsigned int n_vertices = vert_to_send_it->second.size();
3557 const int buffer_size = max_cellid_size * n_vertices;
3558 cellids_send_buffers[i].resize(buffer_size);
3559
3560 // fill the buffer
3561 unsigned int pos = 0;
3562 for (unsigned int j = 0; j < n_vertices; ++j)
3563 {
3564 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
3565 for (unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
3566 {
3567 if (k < cell_id.size())
3568 cellids_send_buffers[i][pos] = cell_id[k];
3569 // if necessary fill up the reserved part of the buffer with an
3570 // invalid value
3571 else
3572 cellids_send_buffers[i][pos] = '-';
3573 }
3574 }
3575
3576 // Send the message
3577 ierr = MPI_Isend(cellids_send_buffers[i].data(),
3578 buffer_size,
3579 MPI_CHAR,
3580 destination,
3581 mpi_tag2,
3582 triangulation.get_communicator(),
3583 &second_requests[i]);
3584 AssertThrowMPI(ierr);
3585 }
3586
3587 // Receive the second message
3588 std::vector<std::vector<char>> cellids_recv_buffers(
3589 vertices_to_recv.size());
3590 vert_to_recv_it = vertices_to_recv.begin();
3591 for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3592 ++vert_to_recv_it, ++i)
3593 {
3594 int source = vert_to_recv_it->first;
3595 const unsigned int n_vertices = vert_to_recv_it->second.size();
3596 const int buffer_size = max_cellid_size * n_vertices;
3597 cellids_recv_buffers[i].resize(buffer_size);
3598
3599 // Receive the message
3600 ierr = MPI_Recv(cellids_recv_buffers[i].data(),
3601 buffer_size,
3602 MPI_CHAR,
3603 source,
3604 mpi_tag2,
3605 triangulation.get_communicator(),
3606 MPI_STATUS_IGNORE);
3607 AssertThrowMPI(ierr);
3608 }
3609
3610
3611 // Match the data received with the required vertices
3612 vert_to_recv_it = vertices_to_recv.begin();
3613 for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3614 ++i, ++vert_to_recv_it)
3615 {
3616 for (unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
3617 {
3618 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
3619 const types::global_vertex_index global_id_recv =
3620 vertices_recv_buffers[i][2 * j + 1];
3621 const std::string cellid_recv(
3622 &cellids_recv_buffers[i][max_cellid_size * j],
3623 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
3624 bool found = false;
3625 typename std::set<active_cell_iterator>::iterator
3626 cell_set_it = missing_vert_cells.begin(),
3627 end_cell_set = missing_vert_cells.end();
3628 for (; (found == false) && (cell_set_it != end_cell_set);
3629 ++cell_set_it)
3630 {
3631 typename std::set<active_cell_iterator>::iterator
3632 candidate_cell =
3633 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
3634 end_cell =
3635 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
3636 for (; candidate_cell != end_cell; ++candidate_cell)
3637 {
3638 std::string current_cellid =
3639 (*candidate_cell)->id().to_string();
3640 current_cellid.resize(max_cellid_size, '-');
3641 if (current_cellid.compare(cellid_recv) == 0)
3642 {
3643 local_to_global_vertex_index
3644 [(*candidate_cell)->vertex_index(local_pos_recv)] =
3645 global_id_recv;
3646 found = true;
3647
3648 break;
3649 }
3650 }
3651 }
3652 }
3653 }
3654#endif
3655
3656 return local_to_global_vertex_index;
3657 }
3658
3659
3660
3661 template <int dim, int spacedim>
3662 void
3665 DynamicSparsityPattern & cell_connectivity)
3666 {
3667 cell_connectivity.reinit(triangulation.n_active_cells(),
3668 triangulation.n_active_cells());
3669
3670 // loop over all cells and their neighbors to build the sparsity
3671 // pattern. note that it's a bit hard to enter all the connections when a
3672 // neighbor has children since we would need to find out which of its
3673 // children is adjacent to the current cell. this problem can be omitted
3674 // if we only do something if the neighbor has no children -- in that case
3675 // it is either on the same or a coarser level than we are. in return, we
3676 // have to add entries in both directions for both cells
3677 for (const auto &cell : triangulation.active_cell_iterators())
3678 {
3679 const unsigned int index = cell->active_cell_index();
3680 cell_connectivity.add(index, index);
3681 for (auto f : cell->face_indices())
3682 if ((cell->at_boundary(f) == false) &&
3683 (cell->neighbor(f)->has_children() == false))
3684 {
3685 const unsigned int other_index =
3686 cell->neighbor(f)->active_cell_index();
3687 cell_connectivity.add(index, other_index);
3688 cell_connectivity.add(other_index, index);
3689 }
3690 }
3691 }
3692
3693
3694
3695 template <int dim, int spacedim>
3696 void
3699 DynamicSparsityPattern & cell_connectivity)
3700 {
3701 std::vector<std::vector<unsigned int>> vertex_to_cell(
3702 triangulation.n_vertices());
3703 for (const auto &cell : triangulation.active_cell_iterators())
3704 {
3705 for (const unsigned int v : cell->vertex_indices())
3706 vertex_to_cell[cell->vertex_index(v)].push_back(
3707 cell->active_cell_index());
3708 }
3709
3710 cell_connectivity.reinit(triangulation.n_active_cells(),
3711 triangulation.n_active_cells());
3712 for (const auto &cell : triangulation.active_cell_iterators())
3713 {
3714 for (const unsigned int v : cell->vertex_indices())
3715 for (unsigned int n = 0;
3716 n < vertex_to_cell[cell->vertex_index(v)].size();
3717 ++n)
3718 cell_connectivity.add(cell->active_cell_index(),
3719 vertex_to_cell[cell->vertex_index(v)][n]);
3720 }
3721 }
3722
3723
3724 template <int dim, int spacedim>
3725 void
3728 const unsigned int level,
3729 DynamicSparsityPattern & cell_connectivity)
3730 {
3731 std::vector<std::vector<unsigned int>> vertex_to_cell(
3732 triangulation.n_vertices());
3734 triangulation.begin(level);
3735 cell != triangulation.end(level);
3736 ++cell)
3737 {
3738 for (const unsigned int v : cell->vertex_indices())
3739 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
3740 }
3741
3742 cell_connectivity.reinit(triangulation.n_cells(level),
3743 triangulation.n_cells(level));
3745 triangulation.begin(level);
3746 cell != triangulation.end(level);
3747 ++cell)
3748 {
3749 for (const unsigned int v : cell->vertex_indices())
3750 for (unsigned int n = 0;
3751 n < vertex_to_cell[cell->vertex_index(v)].size();
3752 ++n)
3753 cell_connectivity.add(cell->index(),
3754 vertex_to_cell[cell->vertex_index(v)][n]);
3755 }
3756 }
3757
3758
3759
3760 template <int dim, int spacedim>
3761 void
3762 partition_triangulation(const unsigned int n_partitions,
3764 const SparsityTools::Partitioner partitioner)
3765 {
3767 &triangulation) == nullptr),
3768 ExcMessage("Objects of type parallel::distributed::Triangulation "
3769 "are already partitioned implicitly and can not be "
3770 "partitioned again explicitly."));
3771
3772 std::vector<unsigned int> cell_weights;
3773
3774 // Get cell weighting if a signal has been attached to the triangulation
3775 if (!triangulation.signals.cell_weight.empty())
3776 {
3777 cell_weights.resize(triangulation.n_active_cells(), 0U);
3778
3779 // In a first step, obtain the weights of the locally owned
3780 // cells. For all others, the weight remains at the zero the
3781 // vector was initialized with above.
3782 for (const auto &cell : triangulation.active_cell_iterators())
3783 if (cell->is_locally_owned())
3784 cell_weights[cell->active_cell_index()] =
3785 triangulation.signals.cell_weight(
3787
3788 // If this is a parallel triangulation, we then need to also
3789 // get the weights for all other cells. We have asserted above
3790 // that this function can't be used for
3791 // parallel::distribute::Triangulation objects, so the only
3792 // ones we have to worry about here are
3793 // parallel::shared::Triangulation
3794 if (const auto shared_tria =
3796 &triangulation))
3797 Utilities::MPI::sum(cell_weights,
3798 shared_tria->get_communicator(),
3799 cell_weights);
3800 }
3801
3802 // Call the other more general function
3803 partition_triangulation(n_partitions,
3804 cell_weights,
3806 partitioner);
3807 }
3808
3809
3810
3811 template <int dim, int spacedim>
3812 void
3813 partition_triangulation(const unsigned int n_partitions,
3814 const std::vector<unsigned int> &cell_weights,
3816 const SparsityTools::Partitioner partitioner)
3817 {
3819 &triangulation) == nullptr),
3820 ExcMessage("Objects of type parallel::distributed::Triangulation "
3821 "are already partitioned implicitly and can not be "
3822 "partitioned again explicitly."));
3823 Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
3824
3825 // check for an easy return
3826 if (n_partitions == 1)
3827 {
3828 for (const auto &cell : triangulation.active_cell_iterators())
3829 cell->set_subdomain_id(0);
3830 return;
3831 }
3832
3833 // we decompose the domain by first
3834 // generating the connection graph of all
3835 // cells with their neighbors, and then
3836 // passing this graph off to METIS.
3837 // finally defer to the other function for
3838 // partitioning and assigning subdomain ids
3839 DynamicSparsityPattern cell_connectivity;
3841
3842 SparsityPattern sp_cell_connectivity;
3843 sp_cell_connectivity.copy_from(cell_connectivity);
3844 partition_triangulation(n_partitions,
3845 cell_weights,
3846 sp_cell_connectivity,
3848 partitioner);
3849 }
3850
3851
3852
3853 template <int dim, int spacedim>
3854 void
3855 partition_triangulation(const unsigned int n_partitions,
3856 const SparsityPattern & cell_connection_graph,
3858 const SparsityTools::Partitioner partitioner)
3859 {
3861 &triangulation) == nullptr),
3862 ExcMessage("Objects of type parallel::distributed::Triangulation "
3863 "are already partitioned implicitly and can not be "
3864 "partitioned again explicitly."));
3865
3866 std::vector<unsigned int> cell_weights;
3867
3868 // Get cell weighting if a signal has been attached to the triangulation
3869 if (!triangulation.signals.cell_weight.empty())
3870 {
3871 cell_weights.resize(triangulation.n_active_cells(), 0U);
3872
3873 // In a first step, obtain the weights of the locally owned
3874 // cells. For all others, the weight remains at the zero the
3875 // vector was initialized with above.
3876 for (const auto &cell : triangulation.active_cell_iterators())
3877 if (cell->is_locally_owned())
3878 cell_weights[cell->active_cell_index()] =
3879 triangulation.signals.cell_weight(
3881
3882 // If this is a parallel triangulation, we then need to also
3883 // get the weights for all other cells. We have asserted above
3884 // that this function can't be used for
3885 // parallel::distribute::Triangulation objects, so the only
3886 // ones we have to worry about here are
3887 // parallel::shared::Triangulation
3888 if (const auto shared_tria =
3890 &triangulation))
3891 Utilities::MPI::sum(cell_weights,
3892 shared_tria->get_communicator(),
3893 cell_weights);
3894 }
3895
3896 // Call the other more general function
3897 partition_triangulation(n_partitions,
3898 cell_weights,
3899 cell_connection_graph,
3901 partitioner);
3902 }
3903
3904
3905
3906 template <int dim, int spacedim>
3907 void
3908 partition_triangulation(const unsigned int n_partitions,
3909 const std::vector<unsigned int> &cell_weights,
3910 const SparsityPattern & cell_connection_graph,
3912 const SparsityTools::Partitioner partitioner)
3913 {
3915 &triangulation) == nullptr),
3916 ExcMessage("Objects of type parallel::distributed::Triangulation "
3917 "are already partitioned implicitly and can not be "
3918 "partitioned again explicitly."));
3919 Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
3920 Assert(cell_connection_graph.n_rows() == triangulation.n_active_cells(),
3921 ExcMessage("Connectivity graph has wrong size"));
3922 Assert(cell_connection_graph.n_cols() == triangulation.n_active_cells(),
3923 ExcMessage("Connectivity graph has wrong size"));
3924
3925 // signal that partitioning is going to happen
3926 triangulation.signals.pre_partition();
3927
3928 // check for an easy return
3929 if (n_partitions == 1)
3930 {
3931 for (const auto &cell : triangulation.active_cell_iterators())
3932 cell->set_subdomain_id(0);
3933 return;
3934 }
3935
3936 // partition this connection graph and get
3937 // back a vector of indices, one per degree
3938 // of freedom (which is associated with a
3939 // cell)
3940 std::vector<unsigned int> partition_indices(triangulation.n_active_cells());
3941 SparsityTools::partition(cell_connection_graph,
3942 cell_weights,
3943 n_partitions,
3944 partition_indices,
3945 partitioner);
3946
3947 // finally loop over all cells and set the subdomain ids
3948 for (const auto &cell : triangulation.active_cell_iterators())
3949 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
3950 }
3951
3952
3953 namespace internal
3954 {
3958 template <class IT>
3959 void
3961 unsigned int & current_proc_idx,
3962 unsigned int & current_cell_idx,
3963 const unsigned int n_active_cells,
3964 const unsigned int n_partitions)
3965 {
3966 if (cell->is_active())
3967 {
3968 while (current_cell_idx >=
3969 std::floor(static_cast<uint_least64_t>(n_active_cells) *
3970 (current_proc_idx + 1) / n_partitions))
3971 ++current_proc_idx;
3972 cell->set_subdomain_id(current_proc_idx);
3973 ++current_cell_idx;
3974 }
3975 else
3976 {
3977 for (unsigned int n = 0; n < cell->n_children(); ++n)
3979 current_proc_idx,
3980 current_cell_idx,
3982 n_partitions);
3983 }
3984 }
3985 } // namespace internal
3986
3987 template <int dim, int spacedim>
3988 void
3989 partition_triangulation_zorder(const unsigned int n_partitions,
3991 const bool group_siblings)
3992 {
3994 &triangulation) == nullptr),
3995 ExcMessage("Objects of type parallel::distributed::Triangulation "
3996 "are already partitioned implicitly and can not be "
3997 "partitioned again explicitly."));
3998 Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
3999
4000 // signal that partitioning is going to happen
4001 triangulation.signals.pre_partition();
4002
4003 // check for an easy return
4004 if (n_partitions == 1)
4005 {
4006 for (const auto &cell : triangulation.active_cell_iterators())
4007 cell->set_subdomain_id(0);
4008 return;
4009 }
4010
4011 // Duplicate the coarse cell reordoring
4012 // as done in p4est
4013 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
4014 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
4015
4016 DynamicSparsityPattern cell_connectivity;
4018 0,
4019 cell_connectivity);
4020 coarse_cell_to_p4est_tree_permutation.resize(triangulation.n_cells(0));
4021 SparsityTools::reorder_hierarchical(cell_connectivity,
4022 coarse_cell_to_p4est_tree_permutation);
4023
4024 p4est_tree_to_coarse_cell_permutation =
4025 Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
4026
4027 unsigned int current_proc_idx = 0;
4028 unsigned int current_cell_idx = 0;
4029 const unsigned int n_active_cells = triangulation.n_active_cells();
4030
4031 // set subdomain id for active cell descendants
4032 // of each coarse cell in permuted order
4033 for (unsigned int idx = 0; idx < triangulation.n_cells(0); ++idx)
4034 {
4035 const unsigned int coarse_cell_idx =
4036 p4est_tree_to_coarse_cell_permutation[idx];
4038 &triangulation, 0, coarse_cell_idx);
4039
4041 current_proc_idx,
4042 current_cell_idx,
4044 n_partitions);
4045 }
4046
4047 // if all children of a cell are active (e.g. we
4048 // have a cell that is refined once and no part
4049 // is refined further), p4est places all of them
4050 // on the same processor. The new owner will be
4051 // the processor with the largest number of children
4052 // (ties are broken by picking the lower rank).
4053 // Duplicate this logic here.
4054 if (group_siblings)
4055 {
4057 cell = triangulation.begin(),
4058 endc = triangulation.end();
4059 for (; cell != endc; ++cell)
4060 {
4061 if (cell->is_active())
4062 continue;
4063 bool all_children_active = true;
4064 std::map<unsigned int, unsigned int> map_cpu_n_cells;
4065 for (unsigned int n = 0; n < cell->n_children(); ++n)
4066 if (!cell->child(n)->is_active())
4067 {
4068 all_children_active = false;
4069 break;
4070 }
4071 else
4072 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
4073
4074 if (!all_children_active)
4075 continue;
4076
4077 unsigned int new_owner = cell->child(0)->subdomain_id();
4078 for (std::map<unsigned int, unsigned int>::iterator it =
4079 map_cpu_n_cells.begin();
4080 it != map_cpu_n_cells.end();
4081 ++it)
4082 if (it->second > map_cpu_n_cells[new_owner])
4083 new_owner = it->first;
4084
4085 for (unsigned int n = 0; n < cell->n_children(); ++n)
4086 cell->child(n)->set_subdomain_id(new_owner);
4087 }
4088 }
4089 }
4090
4091
4092 template <int dim, int spacedim>
4093 void
4095 {
4096 unsigned int n_levels = triangulation.n_levels();
4097 for (int lvl = n_levels - 1; lvl >= 0; --lvl)
4098 {
4100 cell = triangulation.begin(lvl),
4101 endc = triangulation.end(lvl);
4102 for (; cell != endc; ++cell)
4103 {
4104 if (cell->is_active())
4105 cell->set_level_subdomain_id(cell->subdomain_id());
4106 else
4107 {
4108 Assert(cell->child(0)->level_subdomain_id() !=
4111 cell->set_level_subdomain_id(
4112 cell->child(0)->level_subdomain_id());
4113 }
4114 }
4115 }
4116 }
4117
4118
4119
4120 template <int dim, int spacedim>
4121 std::vector<types::subdomain_id>
4123 const std::vector<CellId> & cell_ids)
4124 {
4125 std::vector<types::subdomain_id> subdomain_ids;
4126 subdomain_ids.reserve(cell_ids.size());
4127
4128 if (dynamic_cast<
4130 &triangulation) != nullptr)
4131 {
4132 Assert(false, ExcNotImplemented());
4133 }
4135 *parallel_tria = dynamic_cast<
4137 &triangulation))
4138 {
4139#ifndef DEAL_II_WITH_P4EST
4140 Assert(
4141 false,
4142 ExcMessage(
4143 "You are attempting to use a functionality that is only available "
4144 "if deal.II was configured to use p4est, but cmake did not find a "
4145 "valid p4est library."));
4146#else
4147 // for parallel distributed triangulations, we will ask the p4est oracle
4148 // about the global partitioning of active cells since this information
4149 // is stored on every process
4150 for (const auto &cell_id : cell_ids)
4151 {
4152 // find descendent from coarse quadrant
4153 typename ::internal::p4est::types<dim>::quadrant p4est_cell,
4155
4156 ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
4157 for (const auto &child_index : cell_id.get_child_indices())
4158 {
4159 ::internal::p4est::init_quadrant_children<dim>(
4160 p4est_cell, p4est_children);
4161 p4est_cell =
4162 p4est_children[static_cast<unsigned int>(child_index)];
4163 }
4164
4165 // find owning process, i.e., the subdomain id
4166 const int owner =
4168 const_cast<typename ::internal::p4est::types<dim>::forest
4169 *>(parallel_tria->get_p4est()),
4170 cell_id.get_coarse_cell_id(),
4171 &p4est_cell,
4173 parallel_tria->get_communicator()));
4174
4175 Assert(owner >= 0, ExcMessage("p4est should know the owner."));
4176
4177 subdomain_ids.push_back(owner);
4178 }
4179#endif
4180 }
4181 else if (const parallel::shared::Triangulation<dim, spacedim> *shared_tria =
4183 *>(&triangulation))
4184 {
4185 // for parallel shared triangulations, we need to access true subdomain
4186 // ids which are also valid for artificial cells
4187 const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
4188 shared_tria->get_true_subdomain_ids_of_cells();
4189
4190 for (const auto &cell_id : cell_ids)
4191 {
4192 const unsigned int active_cell_index =
4193 shared_tria->create_cell_iterator(cell_id)->active_cell_index();
4194 subdomain_ids.push_back(
4195 true_subdomain_ids_of_cells[active_cell_index]);
4196 }
4197 }
4198 else
4199 {
4200 // the most general type of triangulation is the serial one. here, all
4201 // subdomain information is directly available
4202 for (const auto &cell_id : cell_ids)
4203 {
4204 subdomain_ids.push_back(
4205 triangulation.create_cell_iterator(cell_id)->subdomain_id());
4206 }
4207 }
4208
4209 return subdomain_ids;
4210 }
4211
4212
4213
4214 template <int dim, int spacedim>
4215 void
4217 std::vector<types::subdomain_id> & subdomain)
4218 {
4219 Assert(subdomain.size() == triangulation.n_active_cells(),
4220 ExcDimensionMismatch(subdomain.size(),
4221 triangulation.n_active_cells()));
4222 for (const auto &cell : triangulation.active_cell_iterators())
4223 subdomain[cell->active_cell_index()] = cell->subdomain_id();
4224 }
4225
4226
4227
4228 template <int dim, int spacedim>
4229 unsigned int
4232 const types::subdomain_id subdomain)
4233 {
4234 unsigned int count = 0;
4235 for (const auto &cell : triangulation.active_cell_iterators())
4236 if (cell->subdomain_id() == subdomain)
4237 ++count;
4238
4239 return count;
4240 }
4241
4242
4243
4244 template <int dim, int spacedim>
4245 std::vector<bool>
4247 {
4248 // start with all vertices
4249 std::vector<bool> locally_owned_vertices =
4250 triangulation.get_used_vertices();
4251
4252 // if the triangulation is distributed, eliminate those that
4253 // are owned by other processors -- either because the vertex is
4254 // on an artificial cell, or because it is on a ghost cell with
4255 // a smaller subdomain
4256 if (const auto *tr = dynamic_cast<
4258 &triangulation))
4259 for (const auto &cell : triangulation.active_cell_iterators())
4260 if (cell->is_artificial() ||
4261 (cell->is_ghost() &&
4262 (cell->subdomain_id() < tr->locally_owned_subdomain())))
4263 for (const unsigned int v : cell->vertex_indices())
4264 locally_owned_vertices[cell->vertex_index(v)] = false;
4265
4266 return locally_owned_vertices;
4267 }
4268
4269
4270
4271 template <int dim, int spacedim>
4272 double
4274 const Mapping<dim, spacedim> & mapping)
4275 {
4276 double min_diameter = std::numeric_limits<double>::max();
4277 for (const auto &cell : triangulation.active_cell_iterators())
4278 if (!cell->is_artificial())
4279 min_diameter = std::min(min_diameter, cell->diameter(mapping));
4280
4281 double global_min_diameter = 0;
4282
4283#ifdef DEAL_II_WITH_MPI
4285 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4286 &triangulation))
4287 global_min_diameter =
4288 Utilities::MPI::min(min_diameter, p_tria->get_communicator());
4289 else
4290#endif
4291 global_min_diameter = min_diameter;
4292
4293 return global_min_diameter;
4294 }
4295
4296
4297
4298 template <int dim, int spacedim>
4299 double
4301 const Mapping<dim, spacedim> & mapping)
4302 {
4303 double max_diameter = 0.;
4304 for (const auto &cell : triangulation.active_cell_iterators())
4305 if (!cell->is_artificial())
4306 max_diameter = std::max(max_diameter, cell->diameter(mapping));
4307
4308 double global_max_diameter = 0;
4309
4310#ifdef DEAL_II_WITH_MPI
4312 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4313 &triangulation))
4314 global_max_diameter =
4315 Utilities::MPI::max(max_diameter, p_tria->get_communicator());
4316 else
4317#endif
4318 global_max_diameter = max_diameter;
4319
4320 return global_max_diameter;
4321 }
4322
4323
4324
4325 namespace internal
4326 {
4327 namespace FixUpDistortedChildCells
4328 {
4329 // compute the mean square
4330 // deviation of the alternating
4331 // forms of the children of the
4332 // given object from that of
4333 // the object itself. for
4334 // objects with
4335 // structdim==spacedim, the
4336 // alternating form is the
4337 // determinant of the jacobian,
4338 // whereas for faces with
4339 // structdim==spacedim-1, the
4340 // alternating form is the
4341 // (signed and scaled) normal
4342 // vector
4343 //
4344 // this average square
4345 // deviation is computed for an
4346 // object where the center node
4347 // has been replaced by the
4348 // second argument to this
4349 // function
4350 template <typename Iterator, int spacedim>
4351 double
4352 objective_function(const Iterator & object,
4353 const Point<spacedim> &object_mid_point)
4354 {
4355 const unsigned int structdim =
4356 Iterator::AccessorType::structure_dimension;
4357 Assert(spacedim == Iterator::AccessorType::dimension,
4359
4360 // everything below is wrong
4361 // if not for the following
4362 // condition
4363 Assert(object->refinement_case() ==
4366 // first calculate the
4367 // average alternating form
4368 // for the parent cell/face
4371 Tensor<spacedim - structdim, spacedim>
4372 parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
4373
4374 for (const unsigned int i : object->vertex_indices())
4375 parent_vertices[i] = object->vertex(i);
4376
4378 parent_vertices, parent_alternating_forms);
4379
4380 const Tensor<spacedim - structdim, spacedim>
4381 average_parent_alternating_form =
4382 std::accumulate(parent_alternating_forms,
4383 parent_alternating_forms +
4386
4387 // now do the same
4388 // computation for the
4389 // children where we use the
4390 // given location for the
4391 // object mid point instead of
4392 // the one the triangulation
4393 // currently reports
4397 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4400
4401 for (unsigned int c = 0; c < object->n_children(); ++c)
4402 for (const unsigned int i : object->child(c)->vertex_indices())
4403 child_vertices[c][i] = object->child(c)->vertex(i);
4404
4405 // replace mid-object
4406 // vertex. note that for
4407 // child i, the mid-object
4408 // vertex happens to have the
4409 // number
4410 // max_children_per_cell-i
4411 for (unsigned int c = 0; c < object->n_children(); ++c)
4413 1] = object_mid_point;
4414
4415 for (unsigned int c = 0; c < object->n_children(); ++c)
4417 child_vertices[c], child_alternating_forms[c]);
4418
4419 // on a uniformly refined
4420 // hypercube object, the child
4421 // alternating forms should
4422 // all be smaller by a factor
4423 // of 2^structdim than the
4424 // ones of the parent. as a
4425 // consequence, we'll use the
4426 // squared deviation from
4427 // this ideal value as an
4428 // objective function
4429 double objective = 0;
4430 for (unsigned int c = 0; c < object->n_children(); ++c)
4431 for (const unsigned int i : object->child(c)->vertex_indices())
4432 objective +=
4433 (child_alternating_forms[c][i] -
4434 average_parent_alternating_form / std::pow(2., 1. * structdim))
4435 .norm_square();
4436
4437 return objective;
4438 }
4439
4440
4446 template <typename Iterator>
4448 get_face_midpoint(const Iterator & object,
4449 const unsigned int f,
4450 std::integral_constant<int, 1>)
4451 {
4452 return object->vertex(f);
4453 }
4454
4455
4456
4462 template <typename Iterator>
4464 get_face_midpoint(const Iterator & object,
4465 const unsigned int f,
4466 std::integral_constant<int, 2>)
4467 {
4468 return object->line(f)->center();
4469 }
4470
4471
4472
4478 template <typename Iterator>
4480 get_face_midpoint(const Iterator & object,
4481 const unsigned int f,
4482 std::integral_constant<int, 3>)
4483 {
4484 return object->face(f)->center();
4485 }
4486
4487
4488
4511 template <typename Iterator>
4512 double
4513 minimal_diameter(const Iterator &object)
4514 {
4515 const unsigned int structdim =
4516 Iterator::AccessorType::structure_dimension;
4517
4518 double diameter = object->diameter();
4519 for (const unsigned int f : object->face_indices())
4520 for (unsigned int e = f + 1; e < object->n_faces(); ++e)
4522 diameter,
4523 get_face_midpoint(object,
4524 f,
4525 std::integral_constant<int, structdim>())
4526 .distance(get_face_midpoint(
4527 object, e, std::integral_constant<int, structdim>())));
4528
4529 return diameter;
4530 }
4531
4532
4533
4538 template <typename Iterator>
4539 bool
4540 fix_up_object(const Iterator &object)
4541 {
4542 const unsigned int structdim =
4543 Iterator::AccessorType::structure_dimension;
4544 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
4545
4546 // right now we can only deal with cells that have been refined
4547 // isotropically because that is the only case where we have a cell
4548 // mid-point that can be moved around without having to consider
4549 // boundary information
4550 Assert(object->has_children(), ExcInternalError());
4551 Assert(object->refinement_case() ==
4554
4555 // get the current location of the object mid-vertex:
4556 Point<spacedim> object_mid_point = object->child(0)->vertex(
4558
4559 // now do a few steepest descent steps to reduce the objective
4560 // function. compute the diameter in the helper function above
4561 unsigned int iteration = 0;
4562 const double diameter = minimal_diameter(object);
4563
4564 // current value of objective function and initial delta
4565 double current_value = objective_function(object, object_mid_point);
4566 double initial_delta = 0;
4567
4568 do
4569 {
4570 // choose a step length that is initially 1/4 of the child
4571 // objects' diameter, and a sequence whose sum does not converge
4572 // (to avoid premature termination of the iteration)
4573 const double step_length = diameter / 4 / (iteration + 1);
4574
4575 // compute the objective function's derivative using a two-sided
4576 // difference formula with eps=step_length/10
4577 Tensor<1, spacedim> gradient;
4578 for (unsigned int d = 0; d < spacedim; ++d)
4579 {
4580 const double eps = step_length / 10;
4581
4583 h[d] = eps / 2;
4584
4585 gradient[d] =
4587 object, project_to_object(object, object_mid_point + h)) -
4589 object, project_to_object(object, object_mid_point - h))) /
4590 eps;
4591 }
4592
4593 // there is nowhere to go
4594 if (gradient.norm() == 0)
4595 break;
4596
4597 // We need to go in direction -gradient. the optimal value of the
4598 // objective function is zero, so assuming that the model is
4599 // quadratic we would have to go -2*val/||gradient|| in this
4600 // direction, make sure we go at most step_length into this
4601 // direction
4602 object_mid_point -=
4603 std::min(2 * current_value / (gradient * gradient),
4604 step_length / gradient.norm()) *
4605 gradient;
4606 object_mid_point = project_to_object(object, object_mid_point);
4607
4608 // compute current value of the objective function
4609 const double previous_value = current_value;
4610 current_value = objective_function(object, object_mid_point);
4611
4612 if (iteration == 0)
4613 initial_delta = (previous_value - current_value);
4614
4615 // stop if we aren't moving much any more
4616 if ((iteration >= 1) &&
4617 ((previous_value - current_value < 0) ||
4618 (std::fabs(previous_value - current_value) <
4619 0.001 * initial_delta)))
4620 break;
4621
4622 ++iteration;
4623 }
4624 while (iteration < 20);
4625
4626 // verify that the new
4627 // location is indeed better
4628 // than the one before. check
4629 // this by comparing whether
4630 // the minimum value of the
4631 // products of parent and
4632 // child alternating forms is
4633 // positive. for cells this
4634 // means that the
4635 // determinants have the same
4636 // sign, for faces that the
4637 // face normals of parent and
4638 // children point in the same
4639 // general direction
4640 double old_min_product, new_min_product;
4641
4644 for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
4645 parent_vertices[i] = object->vertex(i);
4646
4647 Tensor<spacedim - structdim, spacedim>
4648 parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
4650 parent_vertices, parent_alternating_forms);
4651
4655
4656 for (unsigned int c = 0; c < object->n_children(); ++c)
4657 for (const unsigned int i : object->child(c)->vertex_indices())
4658 child_vertices[c][i] = object->child(c)->vertex(i);
4659
4660 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4663
4664 for (unsigned int c = 0; c < object->n_children(); ++c)
4666 child_vertices[c], child_alternating_forms[c]);
4667
4668 old_min_product =
4669 child_alternating_forms[0][0] * parent_alternating_forms[0];
4670 for (unsigned int c = 0; c < object->n_children(); ++c)
4671 for (const unsigned int i : object->child(c)->vertex_indices())
4672 for (const unsigned int j : object->vertex_indices())
4673 old_min_product = std::min<double>(old_min_product,
4674 child_alternating_forms[c][i] *
4675 parent_alternating_forms[j]);
4676
4677 // for the new minimum value,
4678 // replace mid-object
4679 // vertex. note that for child
4680 // i, the mid-object vertex
4681 // happens to have the number
4682 // max_children_per_cell-i
4683 for (unsigned int c = 0; c < object->n_children(); ++c)
4685 1] = object_mid_point;
4686
4687 for (unsigned int c = 0; c < object->n_children(); ++c)
4689 child_vertices[c], child_alternating_forms[c]);
4690
4691 new_min_product =
4692 child_alternating_forms[0][0] * parent_alternating_forms[0];
4693 for (unsigned int c = 0; c < object->n_children(); ++c)
4694 for (const unsigned int i : object->child(c)->vertex_indices())
4695 for (const unsigned int j : object->vertex_indices())
4696 new_min_product = std::min<double>(new_min_product,
4697 child_alternating_forms[c][i] *
4698 parent_alternating_forms[j]);
4699
4700 // if new minimum value is
4701 // better than before, then set the
4702 // new mid point. otherwise
4703 // return this object as one of
4704 // those that can't apparently
4705 // be fixed
4706 if (new_min_product >= old_min_product)
4707 object->child(0)->vertex(
4709 object_mid_point;
4710
4711 // return whether after this
4712 // operation we have an object that
4713 // is well oriented
4714 return (std::max(new_min_product, old_min_product) > 0);
4715 }
4716
4717
4718
4719 // possibly fix up the faces of a cell by moving around its mid-points
4720 template <int dim, int spacedim>
4721 void
4723 const typename ::Triangulation<dim, spacedim>::cell_iterator
4724 &cell,
4725 std::integral_constant<int, dim>,
4726 std::integral_constant<int, spacedim>)
4727 {
4728 // see if we first can fix up some of the faces of this object. We can
4729 // mess with faces if and only if the neighboring cell is not even
4730 // more refined than we are (since in that case the sub-faces have
4731 // themselves children that we can't move around any more). however,
4732 // the latter case shouldn't happen anyway: if the current face is
4733 // distorted but the neighbor is even more refined, then the face had
4734 // been deformed before already, and had been ignored at the time; we
4735 // should then also be able to ignore it this time as well
4736 for (auto f : cell->face_indices())
4737 {
4738 Assert(cell->face(f)->has_children(), ExcInternalError());
4739 Assert(cell->face(f)->refinement_case() ==
4742
4743 bool subface_is_more_refined = false;
4744 for (unsigned int g = 0;
4745 g < GeometryInfo<dim>::max_children_per_face;
4746 ++g)
4747 if (cell->face(f)->child(g)->has_children())
4748 {
4749 subface_is_more_refined = true;
4750 break;
4751 }
4752
4753 if (subface_is_more_refined == true)
4754 continue;
4755
4756 // we finally know that we can do something about this face
4757 fix_up_object(cell->face(f));
4758 }
4759 }
4760 } /* namespace FixUpDistortedChildCells */
4761 } /* namespace internal */
4762
4763
4764 template <int dim, int spacedim>
4768 &distorted_cells,
4769 Triangulation<dim, spacedim> & /*triangulation*/)
4770 {
4771 static_assert(
4772 dim != 1 && spacedim != 1,
4773 "This function is only valid when dim != 1 or spacedim != 1.");
4774 typename Triangulation<dim, spacedim>::DistortedCellList unfixable_subset;
4775
4776 // loop over all cells that we have to fix up
4777 for (typename std::list<
4778 typename Triangulation<dim, spacedim>::cell_iterator>::const_iterator
4779 cell_ptr = distorted_cells.distorted_cells.begin();
4780 cell_ptr != distorted_cells.distorted_cells.end();
4781 ++cell_ptr)
4782 {
4784 *cell_ptr;
4785
4786 Assert(!cell->is_active(),
4787 ExcMessage(
4788 "This function is only valid for a list of cells that "
4789 "have children (i.e., no cell in the list may be active)."));
4790
4792 cell,
4793 std::integral_constant<int, dim>(),
4794 std::integral_constant<int, spacedim>());
4795
4796 // If possible, fix up the object.
4798 unfixable_subset.distorted_cells.push_back(cell);
4799 }
4800
4801 return unfixable_subset;
4802 }
4803
4804
4805
4806 template <int dim, int spacedim>
4807 void
4809 const bool reset_boundary_ids)
4810 {
4811 const auto src_boundary_ids = tria.get_boundary_ids();
4812 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
4813 auto m_it = dst_manifold_ids.begin();
4814 for (const auto b : src_boundary_ids)
4815 {
4816 *m_it = static_cast<types::manifold_id>(b);
4817 ++m_it;
4818 }
4819 const std::vector<types::boundary_id> reset_boundary_id =
4820 reset_boundary_ids ?
4821 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
4822 src_boundary_ids;
4823 map_boundary_to_manifold_ids(src_boundary_ids,
4824 dst_manifold_ids,
4825 tria,
4826 reset_boundary_id);
4827 }
4828
4829
4830
4831 template <int dim, int spacedim>
4832 void
4834 const std::vector<types::boundary_id> &src_boundary_ids,
4835 const std::vector<types::manifold_id> &dst_manifold_ids,
4837 const std::vector<types::boundary_id> &reset_boundary_ids_)
4838 {
4839 AssertDimension(src_boundary_ids.size(), dst_manifold_ids.size());
4840 const auto reset_boundary_ids =
4841 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
4842 AssertDimension(reset_boundary_ids.size(), src_boundary_ids.size());
4843
4844 // in 3d, we not only have to copy boundary ids of faces, but also of edges
4845 // because we see them twice (once from each adjacent boundary face),
4846 // we cannot immediately reset their boundary ids. thus, copy first
4847 // and reset later
4848 if (dim >= 3)
4849 for (const auto &cell : tria.active_cell_iterators())
4850 for (auto f : cell->face_indices())
4851 if (cell->face(f)->at_boundary())
4852 for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
4853 {
4854 const auto bid = cell->face(f)->line(e)->boundary_id();
4855 const unsigned int ind = std::find(src_boundary_ids.begin(),
4856 src_boundary_ids.end(),
4857 bid) -
4858 src_boundary_ids.begin();
4859 if (ind < src_boundary_ids.size())
4860 cell->face(f)->line(e)->set_manifold_id(
4861 dst_manifold_ids[ind]);
4862 }
4863
4864 // now do cells
4865 for (const auto &cell : tria.active_cell_iterators())
4866 for (auto f : cell->face_indices())
4867 if (cell->face(f)->at_boundary())
4868 {
4869 const auto bid = cell->face(f)->boundary_id();
4870 const unsigned int ind =
4871 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
4872 src_boundary_ids.begin();
4873
4874 if (ind < src_boundary_ids.size())
4875 {
4876 // assign the manifold id
4877 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
4878 // then reset boundary id
4879 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
4880 }
4881
4882 if (dim >= 3)
4883 for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
4884 {
4885 const auto bid = cell->face(f)->line(e)->boundary_id();
4886 const unsigned int ind = std::find(src_boundary_ids.begin(),
4887 src_boundary_ids.end(),
4888 bid) -
4889 src_boundary_ids.begin();
4890 if (ind < src_boundary_ids.size())
4891 cell->face(f)->line(e)->set_boundary_id(
4892 reset_boundary_ids[ind]);
4893 }
4894 }
4895 }
4896
4897
4898 template <int dim, int spacedim>
4899 void
4901 const bool compute_face_ids)
4902 {
4904 cell = tria.begin_active(),
4905 endc = tria.end();
4906
4907 for (; cell != endc; ++cell)
4908 {
4909 cell->set_manifold_id(cell->material_id());
4910 if (compute_face_ids == true)
4911 {
4912 for (auto f : cell->face_indices())
4913 {
4914 if (cell->at_boundary(f) == false)
4915 cell->face(f)->set_manifold_id(
4916 std::min(cell->material_id(),
4917 cell->neighbor(f)->material_id()));
4918 else
4919 cell->face(f)->set_manifold_id(cell->material_id());
4920 }
4921 }
4922 }
4923 }
4924
4925
4926 template <int dim, int spacedim>
4927 void
4930 const std::function<types::manifold_id(
4931 const std::set<types::manifold_id> &)> &disambiguation_function,
4932 bool overwrite_only_flat_manifold_ids)
4933 {
4934 // Easy case first:
4935 if (dim == 1)
4936 return;
4937 const unsigned int n_subobjects =
4938 dim == 2 ? tria.n_lines() : tria.n_lines() + tria.n_quads();
4939
4940 // If user index is zero, then it has not been set.
4941 std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
4942 std::vector<unsigned int> backup;
4943 tria.save_user_indices(backup);
4944 tria.clear_user_data();
4945
4946 unsigned next_index = 1;
4947 for (auto &cell : tria.active_cell_iterators())
4948 {
4949 if (dim > 1)
4950 for (unsigned int l = 0; l < cell->n_lines(); ++l)
4951 {
4952 if (cell->line(l)->user_index() == 0)
4953 {
4954 AssertIndexRange(next_index, n_subobjects + 1);
4955 manifold_ids[next_index].insert(cell->manifold_id());
4956 cell->line(l)->set_user_index(next_index++);
4957 }
4958 else
4959 manifold_ids[cell->line(l)->user_index()].insert(
4960 cell->manifold_id());
4961 }
4962 if (dim > 2)
4963 for (unsigned int l = 0; l < cell->n_faces(); ++l)
4964 {
4965 if (cell->quad(l)->user_index() == 0)
4966 {
4967 AssertIndexRange(next_index, n_subobjects + 1);
4968 manifold_ids[next_index].insert(cell->manifold_id());
4969 cell->quad(l)->set_user_index(next_index++);
4970 }
4971 else
4972 manifold_ids[cell->quad(l)->user_index()].insert(
4973 cell->manifold_id());
4974 }
4975 }
4976 for (auto &cell : tria.active_cell_iterators())
4977 {
4978 if (dim > 1)
4979 for (unsigned int l = 0; l < cell->n_lines(); ++l)
4980 {
4981 const auto id = cell->line(l)->user_index();
4982 // Make sure we change the manifold indicator only once
4983 if (id != 0)
4984 {
4985 if (cell->line(l)->manifold_id() ==
4987 overwrite_only_flat_manifold_ids == false)
4988 cell->line(l)->set_manifold_id(
4989 disambiguation_function(manifold_ids[id]));
4990 cell->line(l)->set_user_index(0);
4991 }
4992 }
4993 if (dim > 2)
4994 for (unsigned int l = 0; l < cell->n_faces(); ++l)
4995 {
4996 const auto id = cell->quad(l)->user_index();
4997 // Make sure we change the manifold indicator only once
4998 if (id != 0)
4999 {
5000 if (cell->quad(l)->manifold_id() ==
5002 overwrite_only_flat_manifold_ids == false)
5003 cell->quad(l)->set_manifold_id(
5004 disambiguation_function(manifold_ids[id]));
5005 cell->quad(l)->set_user_index(0);
5006 }
5007 }
5008 }
5009 tria.load_user_indices(backup);
5010 }
5011
5012
5013
5014 template <int dim, int spacedim>
5015 std::pair<unsigned int, double>
5018 {
5019 double max_ratio = 1;
5020 unsigned int index = 0;
5021
5022 for (unsigned int i = 0; i < dim; ++i)
5023 for (unsigned int j = i + 1; j < dim; ++j)
5024 {
5025 unsigned int ax = i % dim;
5026 unsigned int next_ax = j % dim;
5027
5028 double ratio =
5029 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
5030
5031 if (ratio > max_ratio)
5032 {
5033 max_ratio = ratio;
5034 index = ax;
5035 }
5036 else if (1.0 / ratio > max_ratio)
5037 {
5038 max_ratio = 1.0 / ratio;
5039 index = next_ax;
5040 }
5041 }
5042 return std::make_pair(index, max_ratio);
5043 }
5044
5045
5046 template <int dim, int spacedim>
5047 void
5049 const bool isotropic,
5050 const unsigned int max_iterations)
5051 {
5052 unsigned int iter = 0;
5053 bool continue_refinement = true;
5054
5055 while (continue_refinement && (iter < max_iterations))
5056 {
5057 if (max_iterations != numbers::invalid_unsigned_int)
5058 iter++;
5059 continue_refinement = false;
5060
5061 for (const auto &cell : tria.active_cell_iterators())
5062 for (const unsigned int j : cell->face_indices())
5063 if (cell->at_boundary(j) == false &&
5064 cell->neighbor(j)->has_children())
5065 {
5066 if (isotropic)
5067 {
5068 cell->set_refine_flag();
5069 continue_refinement = true;
5070 }
5071 else
5072 continue_refinement |= cell->flag_for_face_refinement(j);
5073 }
5074
5076 }
5077 }
5078
5079 template <int dim, int spacedim>
5080 void
5082 const double max_ratio,
5083 const unsigned int max_iterations)
5084 {
5085 unsigned int iter = 0;
5086 bool continue_refinement = true;
5087
5088 while (continue_refinement && (iter < max_iterations))
5089 {
5090 iter++;
5091 continue_refinement = false;
5092 for (const auto &cell : tria.active_cell_iterators())
5093 {
5094 std::pair<unsigned int, double> info =
5095 GridTools::get_longest_direction<dim, spacedim>(cell);
5096 if (info.second > max_ratio)
5097 {
5098 cell->set_refine_flag(
5099 RefinementCase<dim>::cut_axis(info.first));
5100 continue_refinement = true;
5101 }
5102 }
5104 }
5105 }
5106
5107
5108 template <int dim, int spacedim>
5109 void
5111 const double limit_angle_fraction)
5112 {
5113 if (dim == 1)
5114 return; // Nothing to do
5115
5116 // Check that we don't have hanging nodes
5118 ExcMessage("The input Triangulation cannot "
5119 "have hanging nodes."));
5120
5121
5122 bool has_cells_with_more_than_dim_faces_on_boundary = true;
5123 bool has_cells_with_dim_faces_on_boundary = false;
5124
5125 unsigned int refinement_cycles = 0;
5126
5127 while (has_cells_with_more_than_dim_faces_on_boundary)
5128 {
5129 has_cells_with_more_than_dim_faces_on_boundary = false;
5130
5131 for (const auto &cell : tria.active_cell_iterators())
5132 {
5133 unsigned int boundary_face_counter = 0;
5134 for (auto f : cell->face_indices())
5135 if (cell->face(f)->at_boundary())
5136 boundary_face_counter++;
5137 if (boundary_face_counter > dim)
5138 {
5139 has_cells_with_more_than_dim_faces_on_boundary = true;
5140 break;
5141 }
5142 else if (boundary_face_counter == dim)
5143 has_cells_with_dim_faces_on_boundary = true;
5144 }
5145 if (has_cells_with_more_than_dim_faces_on_boundary)
5146 {
5147 tria.refine_global(1);
5148 refinement_cycles++;
5149 }
5150 }
5151
5152 if (has_cells_with_dim_faces_on_boundary)
5153 {
5154 tria.refine_global(1);
5155 refinement_cycles++;
5156 }
5157 else
5158 {
5159 while (refinement_cycles > 0)
5160 {
5161 for (const auto &cell : tria.active_cell_iterators())
5162 cell->set_coarsen_flag();
5164 refinement_cycles--;
5165 }
5166 return;
5167 }
5168
5169 std::vector<bool> cells_to_remove(tria.n_active_cells(), false);
5170 std::vector<Point<spacedim>> vertices = tria.get_vertices();
5171
5172 std::vector<bool> faces_to_remove(tria.n_raw_faces(), false);
5173
5174 std::vector<CellData<dim>> cells_to_add;
5175 SubCellData subcelldata_to_add;
5176
5177 // Trick compiler for dimension independent things
5178 const unsigned int v0 = 0, v1 = 1, v2 = (dim > 1 ? 2 : 0),
5179 v3 = (dim > 1 ? 3 : 0);
5180
5181 for (const auto &cell : tria.active_cell_iterators())
5182 {
5183 double angle_fraction = 0;
5184 unsigned int vertex_at_corner = numbers::invalid_unsigned_int;
5185
5186 if (dim == 2)
5187 {
5189 p0[spacedim > 1 ? 1 : 0] = 1;
5191 p1[0] = 1;
5192
5193 if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
5194 {
5195 p0 = cell->vertex(v0) - cell->vertex(v2);
5196 p1 = cell->vertex(v3) - cell->vertex(v2);
5197 vertex_at_corner = v2;
5198 }
5199 else if (cell->face(v3)->at_boundary() &&
5200 cell->face(v1)->at_boundary())
5201 {
5202 p0 = cell->vertex(v2) - cell->vertex(v3);
5203 p1 = cell->vertex(v1) - cell->vertex(v3);
5204 vertex_at_corner = v3;
5205 }
5206 else if (cell->face(1)->at_boundary() &&
5207 cell->face(2)->at_boundary())
5208 {
5209 p0 = cell->vertex(v0) - cell->vertex(v1);
5210 p1 = cell->vertex(v3) - cell->vertex(v1);
5211 vertex_at_corner = v1;
5212 }
5213 else if (cell->face(2)->at_boundary() &&
5214 cell->face(0)->at_boundary())
5215 {
5216 p0 = cell->vertex(v2) - cell->vertex(v0);
5217 p1 = cell->vertex(v1) - cell->vertex(v0);
5218 vertex_at_corner = v0;
5219 }
5220 p0 /= p0.norm();
5221 p1 /= p1.norm();
5222 angle_fraction = std::acos(p0 * p1) / numbers::PI;
5223 }
5224 else
5225 {
5226 Assert(false, ExcNotImplemented());
5227 }
5228
5229 if (angle_fraction > limit_angle_fraction)
5230 {
5231 auto flags_removal = [&](unsigned int f1,
5232 unsigned int f2,
5233 unsigned int n1,
5234 unsigned int n2) -> void {
5235 cells_to_remove[cell->active_cell_index()] = true;
5236 cells_to_remove[cell->neighbor(n1)->active_cell_index()] = true;
5237 cells_to_remove[cell->neighbor(n2)->active_cell_index()] = true;
5238
5239 faces_to_remove[cell->face(f1)->index()] = true;
5240 faces_to_remove[cell->face(f2)->index()] = true;
5241
5242 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] = true;
5243 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] = true;
5244 };
5245
5246 auto cell_creation = [&](const unsigned int vv0,
5247 const unsigned int vv1,
5248 const unsigned int f0,
5249 const unsigned int f1,
5250
5251 const unsigned int n0,
5252 const unsigned int v0n0,
5253 const unsigned int v1n0,
5254
5255 const unsigned int n1,
5256 const unsigned int v0n1,
5257 const unsigned int v1n1) {
5258 CellData<dim> c1, c2;
5259 CellData<1> l1, l2;
5260
5261 c1.vertices[v0] = cell->vertex_index(vv0);
5262 c1.vertices[v1] = cell->vertex_index(vv1);
5263 c1.vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
5264 c1.vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
5265
5266 c1.manifold_id = cell->manifold_id();
5267 c1.material_id = cell->material_id();
5268
5269 c2.vertices[v0] = cell->vertex_index(vv0);
5270 c2.vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
5271 c2.vertices[v2] = cell->vertex_index(vv1);
5272 c2.vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
5273
5274 c2.manifold_id = cell->manifold_id();
5275 c2.material_id = cell->material_id();
5276
5277 l1.vertices[0] = cell->vertex_index(vv0);
5278 l1.vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
5279
5280 l1.boundary_id = cell->line(f0)->boundary_id();
5281 l1.manifold_id = cell->line(f0)->manifold_id();
5282 subcelldata_to_add.boundary_lines.push_back(l1);
5283
5284 l2.vertices[0] = cell->vertex_index(vv0);
5285 l2.vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
5286
5287 l2.boundary_id = cell->line(f1)->boundary_id();
5288 l2.manifold_id = cell->line(f1)->manifold_id();
5289 subcelldata_to_add.boundary_lines.push_back(l2);
5290
5291 cells_to_add.push_back(c1);
5292 cells_to_add.push_back(c2);
5293 };
5294
5295 if (dim == 2)
5296 {
5297 switch (vertex_at_corner)
5298 {
5299 case 0:
5300 flags_removal(0, 2, 3, 1);
5301 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
5302 break;
5303 case 1:
5304 flags_removal(1, 2, 3, 0);
5305 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
5306 break;
5307 case 2:
5308 flags_removal(3, 0, 1, 2);
5309 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
5310 break;
5311 case 3:
5312 flags_removal(3, 1, 0, 2);
5313 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
5314 break;
5315 }
5316 }
5317 else
5318 {
5319 Assert(false, ExcNotImplemented());
5320 }
5321 }
5322 }
5323
5324 // if no cells need to be added, then no regularization is necessary.
5325 // Restore things as they were before this function was called.
5326 if (cells_to_add.size() == 0)
5327 {
5328 while (refinement_cycles > 0)
5329 {
5330 for (const auto &cell : tria.active_cell_iterators())
5331 cell->set_coarsen_flag();
5333 refinement_cycles--;
5334 }
5335 return;
5336 }
5337
5338 // add the cells that were not marked as skipped
5339 for (const auto &cell : tria.active_cell_iterators())
5340 {
5341 if (cells_to_remove[cell->active_cell_index()] == false)
5342 {
5343 CellData<dim> c;
5344 for (const unsigned int v : cell->vertex_indices())
5345 c.vertices[v] = cell->vertex_index(v);
5346 c.manifold_id = cell->manifold_id();
5347 c.material_id = cell->material_id();
5348 cells_to_add.push_back(c);
5349 }
5350 }
5351
5352 // Face counter for both dim == 2 and dim == 3
5354 face = tria.begin_active_face(),
5355 endf = tria.end_face();
5356 for (; face != endf; ++face)
5357 if ((face->at_boundary() ||
5358 face->manifold_id() != numbers::flat_manifold_id) &&
5359 faces_to_remove[face->index()] == false)
5360 {
5361 for (unsigned int l = 0; l < face->n_lines(); ++l)
5362 {
5363 CellData<1> line;
5364 if (dim == 2)
5365 {
5366 for (const unsigned int v : face->vertex_indices())
5367 line.vertices[v] = face->vertex_index(v);
5368 line.boundary_id = face->boundary_id();
5369 line.manifold_id = face->manifold_id();
5370 }
5371 else
5372 {
5373 for (const unsigned int v : face->line(l)->vertex_indices())
5374 line.vertices[v] = face->line(l)->vertex_index(v);
5375 line.boundary_id = face->line(l)->boundary_id();
5376 line.manifold_id = face->line(l)->manifold_id();
5377 }
5378 subcelldata_to_add.boundary_lines.push_back(line);
5379 }
5380 if (dim == 3)
5381 {
5382 CellData<2> quad;
5383 for (const unsigned int v : face->vertex_indices())
5384 quad.vertices[v] = face->vertex_index(v);
5385 quad.boundary_id = face->boundary_id();
5386 quad.manifold_id = face->manifold_id();
5387 subcelldata_to_add.boundary_quads.push_back(quad);
5388 }
5389 }
5391 cells_to_add,
5392 subcelldata_to_add);
5394
5395 // Save manifolds
5396 auto manifold_ids = tria.get_manifold_ids();
5397 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
5398 manifolds;
5399 // Set manifolds in new Triangulation
5400 for (const auto manifold_id : manifold_ids)
5402 manifolds[manifold_id] = tria.get_manifold(manifold_id).clone();
5403
5404 tria.clear();
5405
5406 tria.create_triangulation(vertices, cells_to_add, subcelldata_to_add);
5407
5408 // Restore manifolds
5409 for (const auto manifold_id : manifold_ids)
5411 tria.set_manifold(manifold_id, *manifolds[manifold_id]);
5412 }
5413
5414
5415
5416 template <int dim, int spacedim>
5417#ifndef DOXYGEN
5418 std::tuple<
5419 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5420 std::vector<std::vector<Point<dim>>>,
5421 std::vector<std::vector<unsigned int>>>
5422#else
5423 return_type
5424#endif
5426 const Cache<dim, spacedim> & cache,
5427 const std::vector<Point<spacedim>> &points,
5429 &cell_hint)
5430 {
5431 const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
5432 // Splitting the tuple's components
5433 auto &cells = std::get<0>(cqmp);
5434 auto &qpoints = std::get<1>(cqmp);
5435 auto &maps = std::get<2>(cqmp);
5436
5437 return std::make_tuple(std::move(cells),
5438 std::move(qpoints),
5439 std::move(maps));
5440 }
5441
5442
5443
5444 template <int dim, int spacedim>
5445#ifndef DOXYGEN
5446 std::tuple<
5447 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5448 std::vector<std::vector<Point<dim>>>,
5449 std::vector<std::vector<unsigned int>>,
5450 std::vector<unsigned int>>
5451#else
5452 return_type
5453#endif
5455 const Cache<dim, spacedim> & cache,
5456 const std::vector<Point<spacedim>> &points,
5458 &cell_hint)
5459 {
5460 Assert((dim == spacedim),
5461 ExcMessage("Only implemented for dim==spacedim."));
5462
5463 // Alias
5464 namespace bgi = boost::geometry::index;
5465
5466 // Get the mapping
5467 const auto &mapping = cache.get_mapping();
5468
5469 // How many points are here?
5470 const unsigned int np = points.size();
5471
5472 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5473 cells_out;
5474 std::vector<std::vector<Point<dim>>> qpoints_out;
5475 std::vector<std::vector<unsigned int>> maps_out;
5476 std::vector<unsigned int> missing_points_out;
5477
5478 // Now the easy case.
5479 if (np == 0)
5480 return std::make_tuple(std::move(cells_out),
5481 std::move(qpoints_out),
5482 std::move(maps_out),
5483 std::move(missing_points_out));
5484
5485 // For the search we shall use the following tree
5486 const auto &b_tree = cache.get_cell_bounding_boxes_rtree();
5487
5488 // Now make a tree of indices for the points
5489 // [TODO] This would work better with pack_rtree_of_indices, but
5490 // windows does not like it. Build a tree with pairs of point and id
5491 std::vector<std::pair<Point<spacedim>, unsigned int>> points_and_ids(np);
5492 for (unsigned int i = 0; i < np; ++i)
5493 points_and_ids[i] = std::make_pair(points[i], i);
5494 const auto p_tree = pack_rtree(points_and_ids);
5495
5496 // Keep track of all found points
5497 std::vector<bool> found_points(points.size(), false);
5498
5499 // Check if a point was found
5500 const auto already_found = [&found_points](const auto &id) {
5501 AssertIndexRange(id.second, found_points.size());
5502 return found_points[id.second];
5503 };
5504
5505 // check if the given cell was already in the vector of cells before. If so,
5506 // insert in the corresponding vectors the reference point and the id.
5507 // Otherwise append a new entry to all vectors.
5508 const auto store_cell_point_and_id =
5509 [&](
5511 const Point<dim> & ref_point,
5512 const unsigned int &id) {
5513 const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
5514 if (it != cells_out.rend())
5515 {
5516 const auto cell_id =
5517 (cells_out.size() - 1 - (it - cells_out.rbegin()));
5518 qpoints_out[cell_id].emplace_back(ref_point);
5519 maps_out[cell_id].emplace_back(id);
5520 }
5521 else
5522 {
5523 cells_out.emplace_back(cell);
5524 qpoints_out.emplace_back(std::vector<Point<dim>>({ref_point}));
5525 maps_out.emplace_back(std::vector<unsigned int>({id}));
5526 }
5527 };
5528
5529 // Check all points within a given pair of box and cell
5530 const auto check_all_points_within_box = [&](const auto &leaf) {
5531 const auto &box = leaf.first;
5532 const auto &cell_hint = leaf.second;
5533
5534 for (const auto &point_and_id :
5535 p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
5536 bgi::intersects(box)))
5537 {
5538 const auto id = point_and_id.second;
5539 const auto cell_and_ref =
5541 points[id],
5542 cell_hint);
5543 const auto &cell = cell_and_ref.first;
5544 const auto &ref_point = cell_and_ref.second;
5545
5546 if (cell.state() == IteratorState::valid)
5547 store_cell_point_and_id(cell, ref_point, id);
5548 else
5549 missing_points_out.emplace_back(id);
5550
5551 // Don't look anymore for this point
5552 found_points[id] = true;
5553 }
5554 };
5555
5556 // If a hint cell was given, use it
5557 if (cell_hint.state() == IteratorState::valid)
5558 check_all_points_within_box(
5559 std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
5560
5561 // Now loop over all points that have not been found yet
5562 for (unsigned int i = 0; i < np; ++i)
5563 if (found_points[i] == false)
5564 {
5565 // Get the closest cell to this point
5566 const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
5567 // Now checks all points that fall within this box
5568 if (leaf != b_tree.qend())
5569 check_all_points_within_box(*leaf);
5570 else
5571 {
5572 // We should not get here. Throw an error.
5573 Assert(false, ExcInternalError());
5574 }
5575 }
5576 // Now make sure we send out the rest of the points that we did not find.
5577 for (unsigned int i = 0; i < np; ++i)
5578 if (found_points[i] == false)
5579 missing_points_out.emplace_back(i);
5580
5581 // Debug Checking
5582 AssertDimension(cells_out.size(), maps_out.size());
5583 AssertDimension(cells_out.size(), qpoints_out.size());
5584
5585#ifdef DEBUG
5586 unsigned int c = cells_out.size();
5587 unsigned int qps = 0;
5588 // The number of points in all
5589 // the cells must be the same as
5590 // the number of points we
5591 // started off from,
5592 // plus the points which were ignored
5593 for (unsigned int n = 0; n < c; ++n)
5594 {
5595 AssertDimension(qpoints_out[n].size(), maps_out[n].size());
5596 qps += qpoints_out[n].size();
5597 }
5598
5599 Assert(qps + missing_points_out.size() == np,
5600 ExcDimensionMismatch(qps + missing_points_out.size(), np));
5601#endif
5602
5603 return std::make_tuple(std::move(cells_out),
5604 std::move(qpoints_out),
5605 std::move(maps_out),
5606 std::move(missing_points_out));
5607 }
5608
5609
5610
5611 template <int dim, int spacedim>
5612#ifndef DOXYGEN
5613 std::tuple<
5614 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5615 std::vector<std::vector<Point<dim>>>,
5616 std::vector<std::vector<unsigned int>>,
5617 std::vector<std::vector<Point<spacedim>>>,
5618 std::vector<std::vector<unsigned int>>>
5619#else
5620 return_type
5621#endif
5623 const GridTools::Cache<dim, spacedim> & cache,
5624 const std::vector<Point<spacedim>> & points,
5625 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5626 const double tolerance)
5627 {
5628 // run internal function ...
5630 cache, points, global_bboxes, tolerance, false)
5631 .send_components;
5632
5633 // ... and reshuffle the data
5634 std::tuple<
5635 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5636 std::vector<std::vector<Point<dim>>>,
5637 std::vector<std::vector<unsigned int>>,
5638 std::vector<std::vector<Point<spacedim>>>,
5639 std::vector<std::vector<unsigned int>>>
5640 result;
5641
5642 std::pair<int, int> dummy{-1, -1};
5643
5644 for (unsigned int i = 0; i < all.size(); ++i)
5645 {
5646 if (dummy != std::get<0>(all[i]))
5647 {
5648 std::get<0>(result).push_back(
5650 &cache.get_triangulation(),
5651 std::get<0>(all[i]).first,
5652 std::get<0>(all[i]).second});
5653
5654 const unsigned int new_size = std::get<0>(result).size();
5655
5656 std::get<1>(result).resize(new_size);
5657 std::get<2>(result).resize(new_size);
5658 std::get<3>(result).resize(new_size);
5659 std::get<4>(result).resize(new_size);
5660
5661 dummy = std::get<0>(all[i]);
5662 }
5663
5664 std::get<1>(result).back().push_back(
5665 std::get<3>(all[i])); // reference point
5666 std::get<2>(result).back().push_back(std::get<2>(all[i])); // index
5667 std::get<3>(result).back().push_back(std::get<4>(all[i])); // real point
5668 std::get<4>(result).back().push_back(std::get<1>(all[i])); // rank
5669 }
5670
5671 return result;
5672 }
5673
5674
5675
5676 namespace internal
5677 {
5678 template <int spacedim>
5679 std::tuple<std::vector<unsigned int>,
5680 std::vector<unsigned int>,
5681 std::vector<unsigned int>>
5683 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5684 const std::vector<Point<spacedim>> & points)
5685 {
5686 std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
5687 ranks_and_indices.reserve(points.size());
5688
5689 for (unsigned int i = 0; i < points.size(); ++i)
5690 {
5691 const auto &point = points[i];
5692 for (unsigned rank = 0; rank < global_bboxes.size(); ++rank)
5693 for (const auto &box : global_bboxes[rank])
5694 if (box.point_inside(point))
5695 {
5696 ranks_and_indices.emplace_back(rank, i);
5697 break;
5698 }
5699 }
5700
5701 // convert to CRS
5702 std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
5703
5704 std::vector<unsigned int> ranks;
5705 std::vector<unsigned int> ptr;
5706 std::vector<unsigned int> indices;
5707
5708 unsigned int dummy_rank = numbers::invalid_unsigned_int;
5709
5710 for (const auto &i : ranks_and_indices)
5711 {
5712 if (dummy_rank != i.first)
5713 {
5714 dummy_rank = i.first;
5715 ranks.push_back(dummy_rank);
5716 ptr.push_back(indices.size());
5717 }
5718
5719 indices.push_back(i.second);
5720 }
5721 ptr.push_back(indices.size());
5722
5723 return std::make_tuple(std::move(ranks),
5724 std::move(ptr),
5725 std::move(indices));
5726 }
5727
5728
5729
5730 template <int dim, int spacedim>
5731 std::vector<
5732 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5733 Point<dim>>>
5735 const Cache<dim, spacedim> & cache,
5736 const Point<spacedim> & point,
5738 const std::vector<bool> &marked_vertices,
5739 const double tolerance)
5740 {
5741 std::vector<
5742 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5743 Point<dim>>>
5744 locally_owned_active_cells_around_point;
5745
5746 const auto first_cell = GridTools::find_active_cell_around_point(
5747 cache, point, cell_hint, marked_vertices, tolerance);
5748
5749 cell_hint = first_cell.first;
5750 if (cell_hint.state() == IteratorState::valid)
5751 {
5752 const auto active_cells_around_point =
5754 cache.get_mapping(),
5755 cache.get_triangulation(),
5756 point,
5757 tolerance,
5758 first_cell);
5759
5760 locally_owned_active_cells_around_point.reserve(
5761 active_cells_around_point.size());
5762
5763 for (const auto &cell : active_cells_around_point)
5764 if (cell.first->is_locally_owned())
5765 locally_owned_active_cells_around_point.push_back(cell);
5766 }
5767
5768 std::sort(locally_owned_active_cells_around_point.begin(),
5769 locally_owned_active_cells_around_point.end(),
5770 [](const auto &a, const auto &b) { return a.first < b.first; });
5771
5772 return locally_owned_active_cells_around_point;
5773 }
5774
5775
5776
5777 template <int dim, int spacedim>
5778 DistributedComputePointLocationsInternal<dim, spacedim>
5780 const GridTools::Cache<dim, spacedim> & cache,
5781 const std::vector<Point<spacedim>> & points,
5782 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
5783 const double tolerance,
5784 const bool perform_handshake,
5785 const bool enforce_unique_mapping)
5786 {
5787 Assert(!enforce_unique_mapping || perform_handshake, ExcInternalError());
5788
5790
5791 auto &send_components = result.send_components;
5792 auto &send_ranks = result.send_ranks;
5793 auto &send_ptrs = result.send_ptrs;
5794 auto &recv_components = result.recv_components;
5795 auto &recv_ranks = result.recv_ranks;
5796 auto &recv_ptrs = result.recv_ptrs;
5797
5798 const auto potential_owners =
5799 internal::guess_point_owner(global_bboxes, points);
5800
5801 const auto &potential_owners_ranks = std::get<0>(potential_owners);
5802 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
5803 const auto &potential_owners_indices = std::get<2>(potential_owners);
5804
5805 const std::vector<bool> marked_vertices;
5806 auto cell_hint = cache.get_triangulation().begin_active();
5807
5808 const auto translate = [&](const unsigned int other_rank) {
5809 const auto ptr = std::find(potential_owners_ranks.begin(),
5810 potential_owners_ranks.end(),
5811 other_rank);
5812
5813 Assert(ptr != potential_owners_ranks.end(), ExcInternalError());
5814
5815 const auto other_rank_index =
5816 std::distance(potential_owners_ranks.begin(), ptr);
5817
5818 return other_rank_index;
5819 };
5820
5822 [&]() { return potential_owners_ranks; },
5823 [&](const unsigned int other_rank, std::vector<char> &send_buffer) {
5824 const auto other_rank_index = translate(other_rank);
5825
5826 std::vector<std::pair<unsigned int, Point<spacedim>>> temp;
5827 temp.reserve(potential_owners_ptrs[other_rank_index + 1] -
5828 potential_owners_ptrs[other_rank_index]);
5829
5830 for (unsigned int i = potential_owners_ptrs[other_rank_index];
5831 i < potential_owners_ptrs[other_rank_index + 1];
5832 ++i)
5833 temp.emplace_back(potential_owners_indices[i],
5834 points[potential_owners_indices[i]]);
5835
5836 send_buffer = Utilities::pack(temp, false);
5837 },
5838 [&](const unsigned int & other_rank,
5839 const std::vector<char> &recv_buffer,
5840 std::vector<char> & request_buffer) {
5841 const auto recv_buffer_unpacked = Utilities::unpack<
5842 std::vector<std::pair<unsigned int, Point<spacedim>>>>(recv_buffer,
5843 false);
5844
5845 std::vector<unsigned int> request_buffer_temp(
5846 recv_buffer_unpacked.size(), 0);
5847
5848 cell_hint = cache.get_triangulation().begin_active();
5849
5850 for (unsigned int i = 0; i < recv_buffer_unpacked.size(); ++i)
5851 {
5852 const auto &index_and_point = recv_buffer_unpacked[i];
5853
5854 const auto cells_and_reference_positions =
5856 cache,
5857 index_and_point.second,
5858 cell_hint,
5859 marked_vertices,
5860 tolerance);
5861
5862 for (const auto &cell_and_reference_position :
5863 cells_and_reference_positions)
5864 {
5865 send_components.emplace_back(
5866 std::pair<int, int>(
5867 cell_and_reference_position.first->level(),
5868 cell_and_reference_position.first->index()),
5869 other_rank,
5870 index_and_point.first,
5871 cell_and_reference_position.second,
5872 index_and_point.second,
5874
5875 if (enforce_unique_mapping)
5876 break; // in the case of unique mapping, we only need a
5877 // single cell (we take the first)
5878 }
5879
5880 if (perform_handshake)
5881 request_buffer_temp[i] =
5882 enforce_unique_mapping ?
5883 std::min<unsigned int>(
5884 1, cells_and_reference_positions.size()) :
5885 cells_and_reference_positions.size();
5886 }
5887
5888 if (perform_handshake)
5889 request_buffer = Utilities::pack(request_buffer_temp, false);
5890 },
5891 [&](const unsigned int other_rank, std::vector<char> &recv_buffer) {
5892 if (perform_handshake)
5893 {
5894 const auto other_rank_index = translate(other_rank);
5895
5896 recv_buffer =
5897 Utilities::pack(std::vector<unsigned int>(
5898 potential_owners_ptrs[other_rank_index + 1] -
5899 potential_owners_ptrs[other_rank_index]),
5900 false);
5901 }
5902 },
5903 [&](const unsigned int other_rank,
5904 const std::vector<char> &recv_buffer) {
5905 if (perform_handshake)
5906 {
5907 const auto recv_buffer_unpacked =
5908 Utilities::unpack<std::vector<unsigned int>>(recv_buffer,
5909 false);
5910
5911 const auto other_rank_index = translate(other_rank);
5912
5913 for (unsigned int i = 0; i < recv_buffer_unpacked.size(); ++i)
5914 for (unsigned int j = 0; j < recv_buffer_unpacked[i]; ++j)
5915 recv_components.emplace_back(
5916 other_rank,
5917 potential_owners_indices
5918 [i + potential_owners_ptrs[other_rank_index]],
5920 }
5921 });
5922
5924 process, cache.get_triangulation().get_communicator())
5925 .run();
5926
5927 // for unique mapping, we need to modify recv_components and
5928 // send_components consistently
5929 if (enforce_unique_mapping)
5930 {
5931 std::vector<unsigned int> mask_recv(recv_components.size());
5932 std::vector<unsigned int> mask_send(send_components.size());
5933
5934 // set up new recv_components and monitor which entries (we keep
5935 // the entry of the lowest rank) have been eliminated so that we can
5936 // communicate it
5937 auto recv_components_copy = recv_components;
5938 recv_components.clear();
5939
5940 for (unsigned int i = 0; i < recv_components_copy.size(); ++i)
5941 std::get<2>(recv_components_copy[i]) = i;
5942
5943 std::sort(recv_components_copy.begin(),
5944 recv_components_copy.end(),
5945 [&](const auto &a, const auto &b) {
5946 if (std::get<0>(a) != std::get<0>(b)) // rank
5947 return std::get<0>(a) < std::get<0>(b);
5948
5949 return std::get<2>(a) < std::get<2>(b); // enumeration
5950 });
5951
5952 std::vector<bool> unique(points.size(), false);
5953
5954 std::vector<unsigned int> recv_ranks;
5955 std::vector<unsigned int> recv_ptrs;
5956
5957 for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
5958 i < recv_components_copy.size();
5959 ++i)
5960 {
5961 if (dummy != std::get<0>(recv_components_copy[i]))
5962 {
5963 dummy = std::get<0>(recv_components_copy[i]);
5964 recv_ranks.push_back(dummy);
5965 recv_ptrs.push_back(i);
5966 }
5967
5968 if (unique[std::get<1>(recv_components_copy[i])] == false)
5969 {
5970 recv_components.emplace_back(recv_components_copy[i]);
5971 mask_recv[i] = 1;
5972 unique[std::get<1>(recv_components_copy[i])] = true;
5973 }
5974 else
5975 {
5976 mask_recv[i] = 0;
5977 }
5978 }
5979 recv_ptrs.push_back(recv_components_copy.size());
5980
5981 Assert(std::all_of(unique.begin(),
5982 unique.end(),
5983 [](const auto &v) { return v; }),
5985
5986
5987 // prepare send_components so that not needed entries can be
5988 // eliminated later on
5989 auto send_components_copy = send_components;
5990 send_components.clear();
5991
5992 for (unsigned int i = 0; i < send_components_copy.size(); ++i)
5993 std::get<5>(send_components_copy[i]) = i;
5994
5995 std::sort(send_components_copy.begin(),
5996 send_components_copy.end(),
5997 [&](const auto &a, const auto &b) {
5998 if (std::get<1>(a) != std::get<1>(b)) // rank
5999 return std::get<1>(a) < std::get<1>(b);
6000
6001 return std::get<5>(a) < std::get<5>(b); // enumeration
6002 });
6003
6004 std::vector<unsigned int> send_ranks;
6005 std::vector<unsigned int> send_ptrs;
6006
6007 for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
6008 i < send_components_copy.size();
6009 ++i)
6010 {
6011 if (dummy != std::get<1>(send_components_copy[i]))
6012 {
6013 dummy = std::get<1>(send_components_copy[i]);
6014 send_ranks.push_back(dummy);
6015 send_ptrs.push_back(i);
6016 }
6017 }
6018 send_ptrs.push_back(send_components_copy.size());
6019
6020 // perform communication
6021#ifdef DEAL_II_WITH_MPI
6022 std::vector<MPI_Request> req(send_ranks.size() + recv_ranks.size());
6023
6024 for (unsigned int i = 0; i < send_ranks.size(); ++i)
6025 {
6026 const auto ierr =
6027 MPI_Irecv(mask_send.data() + send_ptrs[i],
6028 send_ptrs[i + 1] - send_ptrs[i],
6029 MPI_UNSIGNED,
6030 send_ranks[i],
6034 &req[i]);
6035 AssertThrowMPI(ierr);
6036 }
6037
6038 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
6039 {
6040 const auto ierr =
6041 MPI_Isend(mask_recv.data() + recv_ptrs[i],
6042 recv_ptrs[i + 1] - recv_ptrs[i],
6043 MPI_UNSIGNED,
6044 recv_ranks[i],
6048 &req[i] + send_ranks.size());
6049 AssertThrowMPI(ierr);
6050 }
6051
6052 auto ierr = MPI_Waitall(req.size(), req.data(), MPI_STATUSES_IGNORE);
6053 AssertThrowMPI(ierr);
6054#else
6055 mask_send = mask_recv;
6056#endif
6057
6058 // eliminate not needed entries
6059 for (unsigned int i = 0; i < send_components_copy.size(); ++i)
6060 if (mask_send[i] == 1)
6061 send_components.emplace_back(send_components_copy[i]);
6062 }
6063
6064 if (true)
6065 {
6066 // sort according to rank (and point index and cell) -> make
6067 // deterministic
6068 std::sort(send_components.begin(),
6069 send_components.end(),
6070 [&](const auto &a, const auto &