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