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