Reference documentation for deal.II version GIT relicensing-462-g03c82a14bd 2024-04-21 16:20:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16#include <deal.II/base/mpi.h>
17#include <deal.II/base/mpi.templates.h>
21
22#ifdef DEAL_II_WITH_ARBORX
25#else
26template <int dim, typename Number>
27class BoundingBox;
28
29namespace ArborXWrappers
30{
31 class DistributedTree
32 {
33 public:
34 template <int dim, typename Number>
36 const std::vector<BoundingBox<dim, Number>> &);
37 template <typename QueryType>
38 std::pair<std::vector<std::pair<int, int>>, std::vector<int>>
39 query(const QueryType &queries);
40 };
41 class BoundingBoxIntersectPredicate
42 {};
43} // namespace ArborXWrappers
44#endif
45
46#ifdef DEAL_II_WITH_CGAL
49#endif
50
55
59
61#include <deal.II/fe/fe_q.h>
65
70#include <deal.II/grid/tria.h>
73
82#include <deal.II/lac/vector.h>
84
87
89
90#include <boost/random/mersenne_twister.hpp>
91#include <boost/random/uniform_real_distribution.hpp>
92
93#include <array>
94#include <cmath>
95#include <iostream>
96#include <limits>
97#include <list>
98#include <numeric>
99#include <set>
100#include <tuple>
101#include <unordered_map>
102
104
105
106namespace GridTools
107{
108 // define some transformations
109 namespace internal
110 {
111 template <int spacedim>
112 class Shift
113 {
114 public:
116 : shift(shift)
117 {}
120 {
121 return p + shift;
122 }
123
124 private:
126 };
127
128
129 // Transformation to rotate around one of the cartesian z-axis in 2d.
131 {
132 public:
133 explicit Rotate2d(const double angle)
135 Physics::Transformations::Rotations::rotation_matrix_2d(angle))
136 {}
138 operator()(const Point<2> &p) const
139 {
140 return static_cast<Point<2>>(rotation_matrix * p);
141 }
142
143 private:
145 };
146
147
148 // Transformation to rotate around one of the cartesian axes.
150 {
151 public:
152 Rotate3d(const Tensor<1, 3, double> &axis, const double angle)
154 Physics::Transformations::Rotations::rotation_matrix_3d(axis,
155 angle))
156 {}
157
159 operator()(const Point<3> &p) const
160 {
161 return static_cast<Point<3>>(rotation_matrix * p);
162 }
163
164 private:
166 };
167
168
169 template <int spacedim>
170 class Scale
171 {
172 public:
173 explicit Scale(const double factor)
174 : factor(factor)
175 {}
178 {
179 return p * factor;
180 }
181
182 private:
183 const double factor;
184 };
185 } // namespace internal
186
187
188 template <int dim, int spacedim>
189 void
190 shift(const Tensor<1, spacedim> &shift_vector,
192 {
193 transform(internal::Shift<spacedim>(shift_vector), triangulation);
194 }
195
196
197
198 template <int dim, int spacedim>
199 void
201 {
202 (void)angle;
203 (void)triangulation;
204
205 AssertThrow(false,
207 "GridTools::rotate() is only available for spacedim = 2."));
208 }
209
210
211
212 template <>
213 void
215 {
216 transform(internal::Rotate2d(angle), triangulation);
217 }
218
219
220
221 template <>
222 void
224 {
225 transform(internal::Rotate2d(angle), triangulation);
226 }
227
228
229 template <int dim>
230 void
232 const double angle,
234 {
235 transform(internal::Rotate3d(axis, angle), triangulation);
236 }
237
238
239 template <int dim>
240 void
241 rotate(const double angle,
242 const unsigned int axis,
244 {
245 Assert(axis < 3, ExcMessage("Invalid axis given!"));
246
248 vector[axis] = 1.;
249
250 transform(internal::Rotate3d(vector, angle), triangulation);
251 }
252
253
254 template <int dim, int spacedim>
255 void
256 scale(const double scaling_factor,
258 {
259 Assert(scaling_factor > 0, ExcScalingFactorNotPositive(scaling_factor));
260 transform(internal::Scale<spacedim>(scaling_factor), triangulation);
261 }
262
263
264 namespace internal
265 {
271 inline void
273 const AffineConstraints<double> &constraints,
275 {
276 const unsigned int n_dofs = S.n();
277 const auto op = linear_operator(S);
278 const auto SF = constrained_linear_operator(constraints, op);
280 prec.initialize(S, 1.2);
281
282 SolverControl control(n_dofs, 1.e-10, false, false);
284 SolverCG<Vector<double>> solver(control, mem);
285
286 Vector<double> f(n_dofs);
287
288 const auto constrained_rhs =
289 constrained_right_hand_side(constraints, op, f);
290 solver.solve(SF, u, constrained_rhs, prec);
291
292 constraints.distribute(u);
293 }
294 } // namespace internal
295
296
297 // Implementation for dimensions except 1
298 template <int dim>
299 void
300 laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
302 const Function<dim> *coefficient,
303 const bool solve_for_absolute_positions)
304 {
305 if (dim == 1)
307
308 // first provide everything that is needed for solving a Laplace
309 // equation.
310 FE_Q<dim> q1(1);
311
312 DoFHandler<dim> dof_handler(triangulation);
313 dof_handler.distribute_dofs(q1);
314
315 DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
316 DoFTools::make_sparsity_pattern(dof_handler, dsp);
317 dsp.compress();
318
319 SparsityPattern sparsity_pattern;
320 sparsity_pattern.copy_from(dsp);
321 sparsity_pattern.compress();
322
323 SparseMatrix<double> S(sparsity_pattern);
324
325 QGauss<dim> quadrature(4);
326
329 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
331 reference_cell.template get_default_linear_mapping<dim, dim>(),
332 dof_handler,
333 quadrature,
334 S,
335 coefficient);
336
337 // set up the boundary values for the laplace problem
338 std::array<AffineConstraints<double>, dim> constraints;
339 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
340 new_points.end();
341
342 // Fill these maps using the data given by new_points
343 for (const auto &cell : dof_handler.active_cell_iterators())
344 {
345 // Loop over all vertices of the cell and see if it is listed in the map
346 // given as first argument of the function. We visit vertices multiple
347 // times, so also check that if we have already added a constraint, we
348 // don't do it a second time again.
349 for (const unsigned int vertex_no : cell->vertex_indices())
350 {
351 const unsigned int vertex_index = cell->vertex_index(vertex_no);
352 const Point<dim> &vertex_point = cell->vertex(vertex_no);
353
354 const typename std::map<unsigned int, Point<dim>>::const_iterator
355 map_iter = new_points.find(vertex_index);
356
357 if (map_iter != map_end)
358 for (unsigned int i = 0; i < dim; ++i)
359 if (constraints[i].is_constrained(
360 cell->vertex_dof_index(vertex_no, 0)) == false)
361 {
362 constraints[i].add_constraint(
363 cell->vertex_dof_index(vertex_no, 0),
364 {},
365 (solve_for_absolute_positions ?
366 map_iter->second[i] :
367 map_iter->second[i] - vertex_point[i]));
368 }
369 }
370 }
371
372 for (unsigned int i = 0; i < dim; ++i)
373 constraints[i].close();
374
375 // solve the dim problems with different right hand sides.
376 Vector<double> us[dim];
377 for (unsigned int i = 0; i < dim; ++i)
378 us[i].reinit(dof_handler.n_dofs());
379
380 // solve linear systems in parallel
382 for (unsigned int i = 0; i < dim; ++i)
383 tasks +=
384 Threads::new_task(&internal::laplace_solve, S, constraints[i], us[i]);
385 tasks.join_all();
386
387 // change the coordinates of the points of the triangulation
388 // according to the computed values
389 std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
390 for (const auto &cell : dof_handler.active_cell_iterators())
391 for (const unsigned int vertex_no : cell->vertex_indices())
392 if (vertex_touched[cell->vertex_index(vertex_no)] == false)
393 {
394 Point<dim> &v = cell->vertex(vertex_no);
395
396 const types::global_dof_index dof_index =
397 cell->vertex_dof_index(vertex_no, 0);
398 for (unsigned int i = 0; i < dim; ++i)
399 if (solve_for_absolute_positions)
400 v[i] = us[i](dof_index);
401 else
402 v[i] += us[i](dof_index);
403
404 vertex_touched[cell->vertex_index(vertex_no)] = true;
405 }
406 }
407
412 template <int dim, int spacedim>
413 void
414 distort_random(const double factor,
416 const bool keep_boundary,
417 const unsigned int seed)
418 {
419 // if spacedim>dim we need to make sure that we perturb
420 // points but keep them on
421 // the manifold. however, this isn't implemented right now
422 Assert(spacedim == dim, ExcNotImplemented());
423
424
425 // find the smallest length of the
426 // lines adjacent to the
427 // vertex. take the initial value
428 // to be larger than anything that
429 // might be found: the diameter of
430 // the triangulation, here
431 // estimated by adding up the
432 // diameters of the coarse grid
433 // cells.
434 double almost_infinite_length = 0;
437 cell != triangulation.end(0);
438 ++cell)
439 almost_infinite_length += cell->diameter();
440
441 std::vector<double> minimal_length(triangulation.n_vertices(),
442 almost_infinite_length);
443
444 // also note if a vertex is at the boundary
445 std::vector<bool> at_boundary(keep_boundary ? triangulation.n_vertices() :
446 0,
447 false);
448 // for parallel::shared::Triangulation we need to work on all vertices,
449 // not just the ones related to locally owned cells;
450 const bool is_parallel_shared =
452 &triangulation) != nullptr);
453 for (const auto &cell : triangulation.active_cell_iterators())
454 if (is_parallel_shared || cell->is_locally_owned())
455 {
456 if (dim > 1)
457 {
458 for (unsigned int i = 0; i < cell->n_lines(); ++i)
459 {
461 line = cell->line(i);
462
463 if (keep_boundary && line->at_boundary())
464 {
465 at_boundary[line->vertex_index(0)] = true;
466 at_boundary[line->vertex_index(1)] = true;
467 }
468
469 minimal_length[line->vertex_index(0)] =
470 std::min(line->diameter(),
471 minimal_length[line->vertex_index(0)]);
472 minimal_length[line->vertex_index(1)] =
473 std::min(line->diameter(),
474 minimal_length[line->vertex_index(1)]);
475 }
476 }
477 else // dim==1
478 {
479 if (keep_boundary)
480 for (unsigned int vertex = 0; vertex < 2; ++vertex)
481 if (cell->at_boundary(vertex) == true)
482 at_boundary[cell->vertex_index(vertex)] = true;
483
484 minimal_length[cell->vertex_index(0)] =
485 std::min(cell->diameter(),
486 minimal_length[cell->vertex_index(0)]);
487 minimal_length[cell->vertex_index(1)] =
488 std::min(cell->diameter(),
489 minimal_length[cell->vertex_index(1)]);
490 }
491 }
492
493 // create a random number generator for the interval [-1,1]
494 boost::random::mt19937 rng(seed);
495 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
496
497 // If the triangulation is distributed, we need to
498 // exchange the moved vertices across mpi processes
499 if (auto distributed_triangulation =
502 {
503 const std::vector<bool> locally_owned_vertices =
505 std::vector<bool> vertex_moved(triangulation.n_vertices(), false);
506
507 // Next move vertices on locally owned cells
508 for (const auto &cell : triangulation.active_cell_iterators())
509 if (cell->is_locally_owned())
510 {
511 for (const unsigned int vertex_no : cell->vertex_indices())
512 {
513 const unsigned global_vertex_no =
514 cell->vertex_index(vertex_no);
515
516 // ignore this vertex if we shall keep the boundary and
517 // this vertex *is* at the boundary, if it is already moved
518 // or if another process moves this vertex
519 if ((keep_boundary && at_boundary[global_vertex_no]) ||
520 vertex_moved[global_vertex_no] ||
521 !locally_owned_vertices[global_vertex_no])
522 continue;
523
524 // first compute a random shift vector
525 Point<spacedim> shift_vector;
526 for (unsigned int d = 0; d < spacedim; ++d)
527 shift_vector[d] = uniform_distribution(rng);
528
529 shift_vector *= factor * minimal_length[global_vertex_no] /
530 std::sqrt(shift_vector.square());
531
532 // finally move the vertex
533 cell->vertex(vertex_no) += shift_vector;
534 vertex_moved[global_vertex_no] = true;
535 }
536 }
537
538 distributed_triangulation->communicate_locally_moved_vertices(
539 locally_owned_vertices);
540 }
541 else
542 // if this is a sequential triangulation, we could in principle
543 // use the algorithm above, but we'll use an algorithm that we used
544 // before the parallel::distributed::Triangulation was introduced
545 // in order to preserve backward compatibility
546 {
547 // loop over all vertices and compute their new locations
548 const unsigned int n_vertices = triangulation.n_vertices();
549 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
550 const std::vector<Point<spacedim>> &old_vertex_locations =
552
553 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
554 {
555 // ignore this vertex if we will keep the boundary and
556 // this vertex *is* at the boundary
557 if (keep_boundary && at_boundary[vertex])
558 new_vertex_locations[vertex] = old_vertex_locations[vertex];
559 else
560 {
561 // compute a random shift vector
562 Point<spacedim> shift_vector;
563 for (unsigned int d = 0; d < spacedim; ++d)
564 shift_vector[d] = uniform_distribution(rng);
565
566 shift_vector *= factor * minimal_length[vertex] /
567 std::sqrt(shift_vector.square());
568
569 // record new vertex location
570 new_vertex_locations[vertex] =
571 old_vertex_locations[vertex] + shift_vector;
572 }
573 }
574
575 // now do the actual move of the vertices
576 for (const auto &cell : triangulation.active_cell_iterators())
577 for (const unsigned int vertex_no : cell->vertex_indices())
578 cell->vertex(vertex_no) =
579 new_vertex_locations[cell->vertex_index(vertex_no)];
580 }
581
582 // Correct hanging nodes if necessary
583 if (dim >= 2)
584 {
585 // We do the same as in GridTools::transform
586 //
587 // exclude hanging nodes at the boundaries of artificial cells:
588 // these may belong to ghost cells for which we know the exact
589 // location of vertices, whereas the artificial cell may or may
590 // not be further refined, and so we cannot know whether
591 // the location of the hanging node is correct or not
594 endc = triangulation.end();
595 for (; cell != endc; ++cell)
596 if (!cell->is_artificial())
597 for (const unsigned int face : cell->face_indices())
598 if (cell->face(face)->has_children() &&
599 !cell->face(face)->at_boundary())
600 {
601 // this face has hanging nodes
602 if (dim == 2)
603 cell->face(face)->child(0)->vertex(1) =
604 (cell->face(face)->vertex(0) +
605 cell->face(face)->vertex(1)) /
606 2;
607 else if (dim == 3)
608 {
609 cell->face(face)->child(0)->vertex(1) =
610 .5 * (cell->face(face)->vertex(0) +
611 cell->face(face)->vertex(1));
612 cell->face(face)->child(0)->vertex(2) =
613 .5 * (cell->face(face)->vertex(0) +
614 cell->face(face)->vertex(2));
615 cell->face(face)->child(1)->vertex(3) =
616 .5 * (cell->face(face)->vertex(1) +
617 cell->face(face)->vertex(3));
618 cell->face(face)->child(2)->vertex(3) =
619 .5 * (cell->face(face)->vertex(2) +
620 cell->face(face)->vertex(3));
621
622 // center of the face
623 cell->face(face)->child(0)->vertex(3) =
624 .25 * (cell->face(face)->vertex(0) +
625 cell->face(face)->vertex(1) +
626 cell->face(face)->vertex(2) +
627 cell->face(face)->vertex(3));
628 }
629 }
630 }
631 }
632
633
634
635 template <int dim, template <int, int> class MeshType, int spacedim>
637 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
638 unsigned int find_closest_vertex(const MeshType<dim, spacedim> &mesh,
639 const Point<spacedim> &p,
640 const std::vector<bool> &marked_vertices)
641 {
642 // first get the underlying triangulation from the mesh and determine
643 // vertices and used vertices
645
646 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
647
649 marked_vertices.empty(),
651 marked_vertices.size()));
652
653 // marked_vertices is expected to be a subset of used_vertices. Thus,
654 // comparing the range marked_vertices.begin() to marked_vertices.end() with
655 // the range used_vertices.begin() to used_vertices.end() the element in the
656 // second range must be valid if the element in the first range is valid.
658 marked_vertices.empty() ||
659 std::equal(marked_vertices.begin(),
660 marked_vertices.end(),
661 tria.get_used_vertices().begin(),
662 [](bool p, bool q) { return !p || q; }),
664 "marked_vertices should be a subset of used vertices in the triangulation "
665 "but marked_vertices contains one or more vertices that are not used vertices!"));
666
667 // If marked_indices is empty, consider all used_vertices for finding the
668 // closest vertex to the point. Otherwise, marked_indices is used.
669 const std::vector<bool> &vertices_to_use =
671
672 // At the beginning, the first used vertex is considered to be the closest
673 // one.
674 std::vector<bool>::const_iterator first =
675 std::find(vertices_to_use.begin(), vertices_to_use.end(), true);
676
677 // Assert that at least one vertex is actually used
679
680 unsigned int best_vertex = std::distance(vertices_to_use.begin(), first);
681 double best_dist = (p - vertices[best_vertex]).norm_square();
682
683 // For all remaining vertices, test
684 // whether they are any closer
685 for (unsigned int j = best_vertex + 1; j < vertices.size(); ++j)
686 if (vertices_to_use[j])
687 {
688 const double dist = (p - vertices[j]).norm_square();
689 if (dist < best_dist)
690 {
691 best_vertex = j;
692 best_dist = dist;
693 }
694 }
695
696 return best_vertex;
697 }
698
699
700
701 template <int dim, template <int, int> class MeshType, int spacedim>
703 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
704 unsigned int find_closest_vertex(const Mapping<dim, spacedim> &mapping,
705 const MeshType<dim, spacedim> &mesh,
706 const Point<spacedim> &p,
707 const std::vector<bool> &marked_vertices)
708 {
709 // Take a shortcut in the simple case.
710 if (mapping.preserves_vertex_locations() == true)
712
713 // first get the underlying triangulation from the mesh and determine
714 // vertices and used vertices
716
718
719 Assert(tria.get_vertices().size() == marked_vertices.size() ||
720 marked_vertices.empty(),
722 marked_vertices.size()));
723
724 // marked_vertices is expected to be a subset of used_vertices. Thus,
725 // comparing the range marked_vertices.begin() to marked_vertices.end()
726 // with the range used_vertices.begin() to used_vertices.end() the element
727 // in the second range must be valid if the element in the first range is
728 // valid.
729 Assert(
730 marked_vertices.empty() ||
731 std::equal(marked_vertices.begin(),
732 marked_vertices.end(),
733 tria.get_used_vertices().begin(),
734 [](bool p, bool q) { return !p || q; }),
736 "marked_vertices should be a subset of used vertices in the triangulation "
737 "but marked_vertices contains one or more vertices that are not used vertices!"));
738
739 // Remove from the map unwanted elements.
740 if (marked_vertices.size() != 0)
741 for (auto it = vertices.begin(); it != vertices.end();)
742 {
743 if (marked_vertices[it->first] == false)
744 {
745 it = vertices.erase(it);
746 }
747 else
748 {
749 ++it;
750 }
751 }
752
754 }
755
756
757
758 template <int dim, int spacedim>
759 std::vector<std::vector<Tensor<1, spacedim>>>
762 const std::vector<
764 &vertex_to_cells)
765 {
766 const std::vector<Point<spacedim>> &vertices = mesh.get_vertices();
767 const unsigned int n_vertices = vertex_to_cells.size();
768
769 AssertDimension(vertices.size(), n_vertices);
770
771
772 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
773 n_vertices);
774 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
775 if (mesh.vertex_used(vertex))
776 {
777 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
778 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
779
780 typename std::set<typename Triangulation<dim, spacedim>::
781 active_cell_iterator>::iterator it =
782 vertex_to_cells[vertex].begin();
783 for (unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
784 {
785 vertex_to_cell_centers[vertex][cell] =
786 (*it)->center() - vertices[vertex];
787 vertex_to_cell_centers[vertex][cell] /=
788 vertex_to_cell_centers[vertex][cell].norm();
789 }
790 }
791 return vertex_to_cell_centers;
792 }
793
794
795 namespace internal
796 {
797 template <int spacedim>
798 bool
800 const unsigned int a,
801 const unsigned int b,
802 const Tensor<1, spacedim> &point_direction,
803 const std::vector<Tensor<1, spacedim>> &center_directions)
804 {
805 const double scalar_product_a = center_directions[a] * point_direction;
806 const double scalar_product_b = center_directions[b] * point_direction;
807
808 // The function is supposed to return if a is before b. We are looking
809 // for the alignment of point direction and center direction, therefore
810 // return if the scalar product of a is larger.
811 return (scalar_product_a > scalar_product_b);
812 }
813 } // namespace internal
814
815 template <int dim, template <int, int> class MeshType, int spacedim>
817 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
818#ifndef _MSC_VER
819 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
820#else
821 std::pair<typename ::internal::
822 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
824#endif
827 const MeshType<dim, spacedim> &mesh,
828 const Point<spacedim> &p,
829 const std::vector<
830 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
831 &vertex_to_cells,
832 const std::vector<std::vector<Tensor<1, spacedim>>>
833 &vertex_to_cell_centers,
834 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
835 const std::vector<bool> &marked_vertices,
836 const RTree<std::pair<Point<spacedim>, unsigned int>>
837 &used_vertices_rtree,
838 const double tolerance,
839 const RTree<
840 std::pair<BoundingBox<spacedim>,
842 *relevant_cell_bounding_boxes_rtree)
843 {
844 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
846 cell_and_position;
847 cell_and_position.first = mesh.end();
848
849 // To handle points at the border we keep track of points which are close to
850 // the unit cell:
851 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
853 cell_and_position_approx;
854
855 if (relevant_cell_bounding_boxes_rtree != nullptr &&
856 !relevant_cell_bounding_boxes_rtree->empty())
857 {
858 // create a bounding box around point p with 2*tolerance as side length.
859 const auto bb = BoundingBox<spacedim>(p).create_extended(tolerance);
860
861 if (relevant_cell_bounding_boxes_rtree->qbegin(
862 boost::geometry::index::intersects(bb)) ==
863 relevant_cell_bounding_boxes_rtree->qend())
864 return cell_and_position;
865 }
866
867 bool found_cell = false;
868 bool approx_cell = false;
869
870 unsigned int closest_vertex_index = 0;
871 // ensure closest vertex index is a marked one, otherwise cell (with vertex
872 // 0) might be found even though it is not marked. This is only relevant if
873 // searching with rtree, using find_closest_vertex already can manage not
874 // finding points
875 if (marked_vertices.size() && !used_vertices_rtree.empty())
876 {
877 const auto itr =
878 std::find(marked_vertices.begin(), marked_vertices.end(), true);
879 Assert(itr != marked_vertices.end(),
880 ::ExcMessage("No vertex has been marked!"));
881 closest_vertex_index = std::distance(marked_vertices.begin(), itr);
882 }
883
884 Tensor<1, spacedim> vertex_to_point;
885 auto current_cell = cell_hint;
886
887 // check whether cell has at least one marked vertex
888 const auto cell_marked = [&mesh, &marked_vertices](const auto &cell) {
889 if (marked_vertices.empty())
890 return true;
891
892 if (cell != mesh.active_cell_iterators().end())
893 for (unsigned int i = 0; i < cell->n_vertices(); ++i)
894 if (marked_vertices[cell->vertex_index(i)])
895 return true;
896
897 return false;
898 };
899
900 // check whether any cell in collection is marked
901 const auto any_cell_marked = [&cell_marked](const auto &cells) {
902 return std::any_of(cells.begin(),
903 cells.end(),
904 [&cell_marked](const auto &cell) {
905 return cell_marked(cell);
906 });
907 };
908 (void)any_cell_marked;
909
910 while (found_cell == false)
911 {
912 // First look at the vertices of the cell cell_hint. If it's an
913 // invalid cell, then query for the closest global vertex
914 if (current_cell.state() == IteratorState::valid &&
915 cell_marked(cell_hint))
916 {
917 const auto cell_vertices = mapping.get_vertices(current_cell);
918 const unsigned int closest_vertex =
919 find_closest_vertex_of_cell<dim, spacedim>(current_cell,
920 p,
921 mapping);
922 vertex_to_point = p - cell_vertices[closest_vertex];
923 closest_vertex_index = current_cell->vertex_index(closest_vertex);
924 }
925 else
926 {
927 // For some clang-based compilers and boost versions the call to
928 // RTree::query doesn't compile. Since using an rtree here is just a
929 // performance improvement disabling this branch is OK.
930 // This is fixed in boost in
931 // https://github.com/boostorg/numeric_conversion/commit/50a1eae942effb0a9b90724323ef8f2a67e7984a
932#if defined(DEAL_II_WITH_BOOST_BUNDLED) || \
933 !(defined(__clang_major__) && __clang_major__ >= 16) || \
934 BOOST_VERSION >= 108100
935 if (!used_vertices_rtree.empty())
936 {
937 // If we have an rtree at our disposal, use it.
938 using ValueType = std::pair<Point<spacedim>, unsigned int>;
939 std::function<bool(const ValueType &)> marked;
940 if (marked_vertices.size() == mesh.n_vertices())
941 marked = [&marked_vertices](const ValueType &value) -> bool {
942 return marked_vertices[value.second];
943 };
944 else
945 marked = [](const ValueType &) -> bool { return true; };
946
947 std::vector<std::pair<Point<spacedim>, unsigned int>> res;
948 used_vertices_rtree.query(
949 boost::geometry::index::nearest(p, 1) &&
950 boost::geometry::index::satisfies(marked),
951 std::back_inserter(res));
952
953 // Searching for a point which is located outside the
954 // triangulation results in res.size() = 0
955 Assert(res.size() < 2,
956 ::ExcMessage("There can not be multiple results"));
957
958 if (res.size() > 0)
959 if (any_cell_marked(vertex_to_cells[res[0].second]))
960 closest_vertex_index = res[0].second;
961 }
962 else
963#endif
964 {
965 closest_vertex_index = GridTools::find_closest_vertex(
967 }
968 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
969 }
970
971#ifdef DEBUG
972 {
973 // Double-check if found index is at marked cell
974 Assert(any_cell_marked(vertex_to_cells[closest_vertex_index]),
975 ::ExcMessage("Found non-marked vertex"));
976 }
977#endif
978
979 const double vertex_point_norm = vertex_to_point.norm();
980 if (vertex_point_norm > 0)
981 vertex_to_point /= vertex_point_norm;
982
983 const unsigned int n_neighbor_cells =
984 vertex_to_cells[closest_vertex_index].size();
985
986 // Create a corresponding map of vectors from vertex to cell center
987 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
988
989 for (unsigned int i = 0; i < n_neighbor_cells; ++i)
990 neighbor_permutation[i] = i;
991
992 auto comp = [&](const unsigned int a, const unsigned int b) -> bool {
993 return internal::compare_point_association<spacedim>(
994 a,
995 b,
996 vertex_to_point,
997 vertex_to_cell_centers[closest_vertex_index]);
998 };
999
1000 std::sort(neighbor_permutation.begin(),
1001 neighbor_permutation.end(),
1002 comp);
1003 // It is possible the vertex is close
1004 // to an edge, thus we add a tolerance
1005 // to keep also the "best" cell
1006 double best_distance = tolerance;
1007
1008 // Search all of the cells adjacent to the closest vertex of the cell
1009 // hint. Most likely we will find the point in them.
1010 for (unsigned int i = 0; i < n_neighbor_cells; ++i)
1011 {
1012 try
1013 {
1014 auto cell = vertex_to_cells[closest_vertex_index].begin();
1015 std::advance(cell, neighbor_permutation[i]);
1016
1017 if (!(*cell)->is_artificial())
1018 {
1019 const Point<dim> p_unit =
1020 mapping.transform_real_to_unit_cell(*cell, p);
1021 if ((*cell)->reference_cell().contains_point(p_unit,
1022 tolerance))
1023 {
1024 cell_and_position.first = *cell;
1025 cell_and_position.second = p_unit;
1026 found_cell = true;
1027 approx_cell = false;
1028 break;
1029 }
1030 // The point is not inside this cell: checking how far
1031 // outside it is and whether we want to use this cell as a
1032 // backup if we can't find a cell within which the point
1033 // lies.
1034 const double dist = p_unit.distance(
1035 (*cell)->reference_cell().closest_point(p_unit));
1036 if (dist < best_distance)
1037 {
1038 best_distance = dist;
1039 cell_and_position_approx.first = *cell;
1040 cell_and_position_approx.second = p_unit;
1041 approx_cell = true;
1042 }
1043 }
1044 }
1045 catch (typename Mapping<dim>::ExcTransformationFailed &)
1046 {}
1047 }
1048
1049 if (found_cell == true)
1050 return cell_and_position;
1051 else if (approx_cell == true)
1052 return cell_and_position_approx;
1053
1054 // The first time around, we check for vertices in the hint_cell. If
1055 // that does not work, we set the cell iterator to an invalid one, and
1056 // look for a global vertex close to the point. If that does not work,
1057 // we are in trouble, and just throw an exception.
1058 //
1059 // If we got here, then we did not find the point. If the
1060 // current_cell.state() here is not IteratorState::valid, it means that
1061 // the user did not provide a hint_cell, and at the beginning of the
1062 // while loop we performed an actual global search on the mesh
1063 // vertices. Not finding the point then means the point is outside the
1064 // domain, or that we've had problems with the algorithm above. Try as a
1065 // last resort the other (simpler) algorithm.
1066 if (current_cell.state() != IteratorState::valid)
1068 mapping, mesh, p, marked_vertices, tolerance);
1069
1070 current_cell = typename MeshType<dim, spacedim>::active_cell_iterator();
1071 }
1072 return cell_and_position;
1073 }
1074
1075
1076
1077 template <int dim, int spacedim>
1078 unsigned int
1081 const Point<spacedim> &position,
1083 {
1084 const auto vertices = mapping.get_vertices(cell);
1085 double minimum_distance = position.distance_square(vertices[0]);
1086 unsigned int closest_vertex = 0;
1087 const unsigned int n_vertices = cell->n_vertices();
1088
1089 for (unsigned int v = 1; v < n_vertices; ++v)
1090 {
1091 const double vertex_distance = position.distance_square(vertices[v]);
1092 if (vertex_distance < minimum_distance)
1093 {
1094 closest_vertex = v;
1095 minimum_distance = vertex_distance;
1096 }
1097 }
1098 return closest_vertex;
1099 }
1100
1101
1102
1103 namespace internal
1104 {
1105 namespace BoundingBoxPredicate
1106 {
1107 template <typename MeshType>
1110 std::tuple<
1112 bool> compute_cell_predicate_bounding_box(const typename MeshType::
1113 cell_iterator &parent_cell,
1114 const std::function<bool(
1115 const typename MeshType::
1116 active_cell_iterator &)>
1117 &predicate)
1118 {
1119 bool has_predicate =
1120 false; // Start assuming there's no cells with predicate inside
1121 std::vector<typename MeshType::active_cell_iterator> active_cells;
1122 if (parent_cell->is_active())
1123 active_cells = {parent_cell};
1124 else
1125 // Finding all active cells descendants of the current one (or the
1126 // current one if it is active)
1127 active_cells = get_active_child_cells<MeshType>(parent_cell);
1128
1129 const unsigned int spacedim = MeshType::space_dimension;
1130
1131 // Looking for the first active cell which has the property predicate
1132 unsigned int i = 0;
1133 while (i < active_cells.size() && !predicate(active_cells[i]))
1134 ++i;
1135
1136 // No active cells or no active cells with property
1137 if (active_cells.empty() || i == active_cells.size())
1138 {
1140 return std::make_tuple(bbox, has_predicate);
1141 }
1142
1143 // The two boundary points defining the boundary box
1144 Point<spacedim> maxp = active_cells[i]->vertex(0);
1145 Point<spacedim> minp = active_cells[i]->vertex(0);
1146
1147 for (; i < active_cells.size(); ++i)
1148 if (predicate(active_cells[i]))
1149 for (const unsigned int v : active_cells[i]->vertex_indices())
1150 for (unsigned int d = 0; d < spacedim; ++d)
1151 {
1152 minp[d] = std::min(minp[d], active_cells[i]->vertex(v)[d]);
1153 maxp[d] = std::max(maxp[d], active_cells[i]->vertex(v)[d]);
1154 }
1155
1156 has_predicate = true;
1157 BoundingBox<spacedim> bbox(std::make_pair(minp, maxp));
1158 return std::make_tuple(bbox, has_predicate);
1159 }
1160 } // namespace BoundingBoxPredicate
1161 } // namespace internal
1162
1163
1164
1165 template <typename MeshType>
1167 std::
1168 vector<BoundingBox<MeshType::space_dimension>> compute_mesh_predicate_bounding_box(
1169 const MeshType &mesh,
1170 const std::function<bool(const typename MeshType::active_cell_iterator &)>
1171 &predicate,
1172 const unsigned int refinement_level,
1173 const bool allow_merge,
1174 const unsigned int max_boxes)
1175 {
1176 // Algorithm brief description: begin with creating bounding boxes of all
1177 // cells at refinement_level (and coarser levels if there are active cells)
1178 // which have the predicate property. These are then merged
1179
1180 Assert(
1181 refinement_level <= mesh.n_levels(),
1182 ExcMessage(
1183 "Error: refinement level is higher then total levels in the triangulation!"));
1184
1185 const unsigned int spacedim = MeshType::space_dimension;
1186 std::vector<BoundingBox<spacedim>> bounding_boxes;
1187
1188 // Creating a bounding box for all active cell on coarser level
1189
1190 for (unsigned int i = 0; i < refinement_level; ++i)
1191 for (const typename MeshType::cell_iterator &cell :
1192 mesh.active_cell_iterators_on_level(i))
1193 {
1194 bool has_predicate = false;
1196 std::tie(bbox, has_predicate) =
1197 internal::BoundingBoxPredicate::compute_cell_predicate_bounding_box<
1198 MeshType>(cell, predicate);
1199 if (has_predicate)
1200 bounding_boxes.push_back(bbox);
1201 }
1202
1203 // Creating a Bounding Box for all cells on the chosen refinement_level
1204 for (const typename MeshType::cell_iterator &cell :
1205 mesh.cell_iterators_on_level(refinement_level))
1206 {
1207 bool has_predicate = false;
1209 std::tie(bbox, has_predicate) =
1210 internal::BoundingBoxPredicate::compute_cell_predicate_bounding_box<
1211 MeshType>(cell, predicate);
1212 if (has_predicate)
1213 bounding_boxes.push_back(bbox);
1214 }
1215
1216 if (!allow_merge)
1217 // If merging is not requested return the created bounding_boxes
1218 return bounding_boxes;
1219 else
1220 {
1221 // Merging part of the algorithm
1222 // Part 1: merging neighbors
1223 // This array stores the indices of arrays we have already merged
1224 std::vector<unsigned int> merged_boxes_idx;
1225 bool found_neighbors = true;
1226
1227 // We merge only neighbors which can be expressed by a single bounding
1228 // box e.g. in 1d [0,1] and [1,2] can be described with [0,2] without
1229 // losing anything
1230 while (found_neighbors)
1231 {
1232 found_neighbors = false;
1233 for (unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
1234 {
1235 if (std::find(merged_boxes_idx.begin(),
1236 merged_boxes_idx.end(),
1237 i) == merged_boxes_idx.end())
1238 for (unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
1239 if (std::find(merged_boxes_idx.begin(),
1240 merged_boxes_idx.end(),
1241 j) == merged_boxes_idx.end() &&
1242 bounding_boxes[i].get_neighbor_type(
1243 bounding_boxes[j]) ==
1245 {
1246 bounding_boxes[i].merge_with(bounding_boxes[j]);
1247 merged_boxes_idx.push_back(j);
1248 found_neighbors = true;
1249 }
1250 }
1251 }
1252
1253 // Copying the merged boxes into merged_b_boxes
1254 std::vector<BoundingBox<spacedim>> merged_b_boxes;
1255 for (unsigned int i = 0; i < bounding_boxes.size(); ++i)
1256 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
1257 merged_boxes_idx.end())
1258 merged_b_boxes.push_back(bounding_boxes[i]);
1259
1260 // Part 2: if there are too many bounding boxes, merging smaller boxes
1261 // This has sense only in dimension 2 or greater, since in dimension 1,
1262 // neighboring intervals can always be merged without problems
1263 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
1264 {
1265 std::vector<double> volumes;
1266 for (unsigned int i = 0; i < merged_b_boxes.size(); ++i)
1267 volumes.push_back(merged_b_boxes[i].volume());
1268
1269 while (merged_b_boxes.size() > max_boxes)
1270 {
1271 unsigned int min_idx =
1272 std::min_element(volumes.begin(), volumes.end()) -
1273 volumes.begin();
1274 volumes.erase(volumes.begin() + min_idx);
1275 // Finding a neighbor
1276 bool not_removed = true;
1277 for (unsigned int i = 0;
1278 i < merged_b_boxes.size() && not_removed;
1279 ++i)
1280 // We merge boxes if we have "attached" or "mergeable"
1281 // neighbors, even though mergeable should be dealt with in
1282 // Part 1
1283 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
1284 merged_b_boxes[min_idx]) ==
1286 merged_b_boxes[i].get_neighbor_type(
1287 merged_b_boxes[min_idx]) ==
1289 {
1290 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
1291 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
1292 not_removed = false;
1293 }
1294 Assert(!not_removed,
1295 ExcMessage("Error: couldn't merge bounding boxes!"));
1296 }
1297 }
1298 Assert(merged_b_boxes.size() <= max_boxes,
1299 ExcMessage(
1300 "Error: couldn't reach target number of bounding boxes!"));
1301 return merged_b_boxes;
1302 }
1303 }
1304
1305
1306
1307 template <int spacedim>
1308#ifndef DOXYGEN
1309 std::tuple<std::vector<std::vector<unsigned int>>,
1310 std::map<unsigned int, unsigned int>,
1311 std::map<unsigned int, std::vector<unsigned int>>>
1312#else
1313 return_type
1314#endif
1316 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1317 const std::vector<Point<spacedim>> &points)
1318 {
1319 unsigned int n_procs = global_bboxes.size();
1320 std::vector<std::vector<unsigned int>> point_owners(n_procs);
1321 std::map<unsigned int, unsigned int> map_owners_found;
1322 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
1323
1324 unsigned int n_points = points.size();
1325 for (unsigned int pt = 0; pt < n_points; ++pt)
1326 {
1327 // Keep track of how many processes we guess to own the point
1328 std::vector<unsigned int> owners_found;
1329 // Check in which other processes the point might be
1330 for (unsigned int rk = 0; rk < n_procs; ++rk)
1331 {
1332 for (const BoundingBox<spacedim> &bbox : global_bboxes[rk])
1333 if (bbox.point_inside(points[pt]))
1334 {
1335 point_owners[rk].emplace_back(pt);
1336 owners_found.emplace_back(rk);
1337 break; // We can check now the next process
1338 }
1339 }
1340 Assert(owners_found.size() > 0,
1341 ExcMessage("No owners found for the point " +
1342 std::to_string(pt)));
1343 if (owners_found.size() == 1)
1344 map_owners_found[pt] = owners_found[0];
1345 else
1346 // Multiple owners
1347 map_owners_guessed[pt] = owners_found;
1348 }
1349
1350 return std::make_tuple(std::move(point_owners),
1351 std::move(map_owners_found),
1352 std::move(map_owners_guessed));
1353 }
1354
1355 template <int spacedim>
1356#ifndef DOXYGEN
1357 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1358 std::map<unsigned int, unsigned int>,
1359 std::map<unsigned int, std::vector<unsigned int>>>
1360#else
1361 return_type
1362#endif
1364 const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1365 const std::vector<Point<spacedim>> &points)
1366 {
1367 std::map<unsigned int, std::vector<unsigned int>> point_owners;
1368 std::map<unsigned int, unsigned int> map_owners_found;
1369 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
1370 std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> search_result;
1371
1372 unsigned int n_points = points.size();
1373 for (unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
1374 {
1375 search_result.clear(); // clearing last output
1376
1377 // Running tree search
1378 covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
1379 std::back_inserter(search_result));
1380
1381 // Keep track of how many processes we guess to own the point
1382 std::set<unsigned int> owners_found;
1383 // Check in which other processes the point might be
1384 for (const auto &rank_bbox : search_result)
1385 {
1386 // Try to add the owner to the owners found,
1387 // and check if it was already present
1388 const bool pt_inserted = owners_found.insert(pt_n).second;
1389 if (pt_inserted)
1390 point_owners[rank_bbox.second].emplace_back(pt_n);
1391 }
1392 Assert(owners_found.size() > 0,
1393 ExcMessage("No owners found for the point " +
1394 std::to_string(pt_n)));
1395 if (owners_found.size() == 1)
1396 map_owners_found[pt_n] = *owners_found.begin();
1397 else
1398 // Multiple owners
1399 std::copy(owners_found.begin(),
1400 owners_found.end(),
1401 std::back_inserter(map_owners_guessed[pt_n]));
1402 }
1403
1404 return std::make_tuple(std::move(point_owners),
1405 std::move(map_owners_found),
1406 std::move(map_owners_guessed));
1407 }
1408
1409
1410
1411 template <int dim, int spacedim>
1412 std::map<unsigned int, types::global_vertex_index>
1415 {
1416 std::map<unsigned int, types::global_vertex_index>
1417 local_to_global_vertex_index;
1418
1419#ifndef DEAL_II_WITH_MPI
1420
1421 // without MPI, this function doesn't make sense because on cannot
1422 // use parallel::distributed::Triangulation in any meaningful
1423 // way
1424 (void)triangulation;
1425 Assert(false,
1426 ExcMessage("This function does not make any sense "
1427 "for parallel::distributed::Triangulation "
1428 "objects if you do not have MPI enabled."));
1429
1430#else
1431
1432 using active_cell_iterator =
1434 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
1436
1437 // Create a local index for the locally "owned" vertices
1438 types::global_vertex_index next_index = 0;
1439 unsigned int max_cellid_size = 0;
1440 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
1441 vertices_added;
1442 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
1443 std::map<types::subdomain_id,
1444 std::vector<std::tuple<types::global_vertex_index,
1446 std::string>>>
1447 vertices_to_send;
1448 std::set<active_cell_iterator> missing_vert_cells;
1449 std::set<unsigned int> used_vertex_index;
1450 for (const auto &cell : triangulation.active_cell_iterators())
1451 {
1452 if (cell->is_locally_owned())
1453 {
1454 for (const unsigned int i : cell->vertex_indices())
1455 {
1456 types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
1457 for (const auto &adjacent_cell :
1458 vertex_to_cell[cell->vertex_index(i)])
1459 lowest_subdomain_id = std::min(lowest_subdomain_id,
1460 adjacent_cell->subdomain_id());
1461
1462 // See if this process "owns" this vertex
1463 if (lowest_subdomain_id == cell->subdomain_id())
1464 {
1465 // Check that the vertex we are working on is a vertex that
1466 // has not been dealt with yet
1467 if (used_vertex_index.find(cell->vertex_index(i)) ==
1468 used_vertex_index.end())
1469 {
1470 // Set the local index
1471 local_to_global_vertex_index[cell->vertex_index(i)] =
1472 next_index++;
1473
1474 // Store the information that will be sent to the
1475 // adjacent cells on other subdomains
1476 for (const auto &adjacent_cell :
1477 vertex_to_cell[cell->vertex_index(i)])
1478 if (adjacent_cell->subdomain_id() !=
1479 cell->subdomain_id())
1480 {
1481 std::pair<types::subdomain_id,
1483 tmp(adjacent_cell->subdomain_id(),
1484 cell->vertex_index(i));
1485 if (vertices_added.find(tmp) ==
1486 vertices_added.end())
1487 {
1488 vertices_to_send[adjacent_cell
1489 ->subdomain_id()]
1490 .emplace_back(i,
1491 cell->vertex_index(i),
1492 cell->id().to_string());
1493 if (cell->id().to_string().size() >
1494 max_cellid_size)
1495 max_cellid_size =
1496 cell->id().to_string().size();
1497 vertices_added.insert(tmp);
1498 }
1499 }
1500 used_vertex_index.insert(cell->vertex_index(i));
1501 }
1502 }
1503 else
1504 {
1505 // We don't own the vertex so we will receive its global
1506 // index
1507 vertices_to_recv[lowest_subdomain_id].insert(
1508 cell->vertex_index(i));
1509 missing_vert_cells.insert(cell);
1510 }
1511 }
1512 }
1513
1514 // Some hanging nodes are vertices of ghost cells. They need to be
1515 // received.
1516 if (cell->is_ghost())
1517 {
1518 for (const unsigned int i : cell->face_indices())
1519 {
1520 if (cell->at_boundary(i) == false)
1521 {
1522 if (cell->neighbor(i)->is_active())
1523 {
1524 typename Triangulation<dim,
1525 spacedim>::active_cell_iterator
1526 adjacent_cell = cell->neighbor(i);
1527 if ((adjacent_cell->is_locally_owned()))
1528 {
1529 types::subdomain_id adj_subdomain_id =
1530 adjacent_cell->subdomain_id();
1531 if (cell->subdomain_id() < adj_subdomain_id)
1532 for (unsigned int j = 0;
1533 j < cell->face(i)->n_vertices();
1534 ++j)
1535 {
1536 vertices_to_recv[cell->subdomain_id()].insert(
1537 cell->face(i)->vertex_index(j));
1538 missing_vert_cells.insert(cell);
1539 }
1540 }
1541 }
1542 }
1543 }
1544 }
1545 }
1546
1547 // Get the size of the largest CellID string
1548 max_cellid_size =
1550
1551 // Make indices global by getting the number of vertices owned by each
1552 // processors and shifting the indices accordingly
1554 int ierr = MPI_Exscan(&next_index,
1555 &shift,
1556 1,
1558 MPI_SUM,
1560 AssertThrowMPI(ierr);
1561
1562 for (auto &global_index_it : local_to_global_vertex_index)
1563 global_index_it.second += shift;
1564
1565
1566 const int mpi_tag = Utilities::MPI::internal::Tags::
1568 const int mpi_tag2 = Utilities::MPI::internal::Tags::
1570
1571
1572 // In a first message, send the global ID of the vertices and the local
1573 // positions in the cells. In a second messages, send the cell ID as a
1574 // resize string. This is done in two messages so that types are not mixed
1575
1576 // Send the first message
1577 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
1578 vertices_to_send.size());
1579 std::vector<MPI_Request> first_requests(vertices_to_send.size());
1580 typename std::map<types::subdomain_id,
1581 std::vector<std::tuple<types::global_vertex_index,
1583 std::string>>>::iterator
1584 vert_to_send_it = vertices_to_send.begin(),
1585 vert_to_send_end = vertices_to_send.end();
1586 for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
1587 ++vert_to_send_it, ++i)
1588 {
1589 int destination = vert_to_send_it->first;
1590 const unsigned int n_vertices = vert_to_send_it->second.size();
1591 const int buffer_size = 2 * n_vertices;
1592 vertices_send_buffers[i].resize(buffer_size);
1593
1594 // fill the buffer
1595 for (unsigned int j = 0; j < n_vertices; ++j)
1596 {
1597 vertices_send_buffers[i][2 * j] =
1598 std::get<0>(vert_to_send_it->second[j]);
1599 vertices_send_buffers[i][2 * j + 1] =
1600 local_to_global_vertex_index[std::get<1>(
1601 vert_to_send_it->second[j])];
1602 }
1603
1604 // Send the message
1605 ierr = MPI_Isend(vertices_send_buffers[i].data(),
1606 buffer_size,
1608 destination,
1609 mpi_tag,
1611 &first_requests[i]);
1612 AssertThrowMPI(ierr);
1613 }
1614
1615 // Receive the first message
1616 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
1617 vertices_to_recv.size());
1618 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
1619 vert_to_recv_it = vertices_to_recv.begin(),
1620 vert_to_recv_end = vertices_to_recv.end();
1621 for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1622 ++vert_to_recv_it, ++i)
1623 {
1624 int source = vert_to_recv_it->first;
1625 const unsigned int n_vertices = vert_to_recv_it->second.size();
1626 const int buffer_size = 2 * n_vertices;
1627 vertices_recv_buffers[i].resize(buffer_size);
1628
1629 // Receive the message
1630 ierr = MPI_Recv(vertices_recv_buffers[i].data(),
1631 buffer_size,
1633 source,
1634 mpi_tag,
1636 MPI_STATUS_IGNORE);
1637 AssertThrowMPI(ierr);
1638 }
1639
1640 // At this point, wait for all of the isend operations to finish:
1641 MPI_Waitall(first_requests.size(),
1642 first_requests.data(),
1643 MPI_STATUSES_IGNORE);
1644
1645
1646 // Send second message
1647 std::vector<std::vector<char>> cellids_send_buffers(
1648 vertices_to_send.size());
1649 std::vector<MPI_Request> second_requests(vertices_to_send.size());
1650 vert_to_send_it = vertices_to_send.begin();
1651 for (unsigned int i = 0; vert_to_send_it != vert_to_send_end;
1652 ++vert_to_send_it, ++i)
1653 {
1654 int destination = vert_to_send_it->first;
1655 const unsigned int n_vertices = vert_to_send_it->second.size();
1656 const int buffer_size = max_cellid_size * n_vertices;
1657 cellids_send_buffers[i].resize(buffer_size);
1658
1659 // fill the buffer
1660 unsigned int pos = 0;
1661 for (unsigned int j = 0; j < n_vertices; ++j)
1662 {
1663 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
1664 for (unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
1665 {
1666 if (k < cell_id.size())
1667 cellids_send_buffers[i][pos] = cell_id[k];
1668 // if necessary fill up the reserved part of the buffer with an
1669 // invalid value
1670 else
1671 cellids_send_buffers[i][pos] = '-';
1672 }
1673 }
1674
1675 // Send the message
1676 ierr = MPI_Isend(cellids_send_buffers[i].data(),
1677 buffer_size,
1678 MPI_CHAR,
1679 destination,
1680 mpi_tag2,
1682 &second_requests[i]);
1683 AssertThrowMPI(ierr);
1684 }
1685
1686 // Receive the second message
1687 std::vector<std::vector<char>> cellids_recv_buffers(
1688 vertices_to_recv.size());
1689 vert_to_recv_it = vertices_to_recv.begin();
1690 for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1691 ++vert_to_recv_it, ++i)
1692 {
1693 int source = vert_to_recv_it->first;
1694 const unsigned int n_vertices = vert_to_recv_it->second.size();
1695 const int buffer_size = max_cellid_size * n_vertices;
1696 cellids_recv_buffers[i].resize(buffer_size);
1697
1698 // Receive the message
1699 ierr = MPI_Recv(cellids_recv_buffers[i].data(),
1700 buffer_size,
1701 MPI_CHAR,
1702 source,
1703 mpi_tag2,
1705 MPI_STATUS_IGNORE);
1706 AssertThrowMPI(ierr);
1707 }
1708
1709
1710 // Match the data received with the required vertices
1711 vert_to_recv_it = vertices_to_recv.begin();
1712 for (unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1713 ++i, ++vert_to_recv_it)
1714 {
1715 for (unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
1716 {
1717 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
1718 const types::global_vertex_index global_id_recv =
1719 vertices_recv_buffers[i][2 * j + 1];
1720 const std::string cellid_recv(
1721 &cellids_recv_buffers[i][max_cellid_size * j],
1722 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
1723 bool found = false;
1724 typename std::set<active_cell_iterator>::iterator
1725 cell_set_it = missing_vert_cells.begin(),
1726 end_cell_set = missing_vert_cells.end();
1727 for (; (found == false) && (cell_set_it != end_cell_set);
1728 ++cell_set_it)
1729 {
1730 typename std::set<active_cell_iterator>::iterator
1731 candidate_cell =
1732 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
1733 end_cell =
1734 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
1735 for (; candidate_cell != end_cell; ++candidate_cell)
1736 {
1737 std::string current_cellid =
1738 (*candidate_cell)->id().to_string();
1739 current_cellid.resize(max_cellid_size, '-');
1740 if (current_cellid.compare(cellid_recv) == 0)
1741 {
1742 local_to_global_vertex_index
1743 [(*candidate_cell)->vertex_index(local_pos_recv)] =
1744 global_id_recv;
1745 found = true;
1746
1747 break;
1748 }
1749 }
1750 }
1751 }
1752 }
1753
1754 // At this point, wait for all of the isend operations of the second round
1755 // to finish:
1756 MPI_Waitall(second_requests.size(),
1757 second_requests.data(),
1758 MPI_STATUSES_IGNORE);
1759#endif
1760
1761 return local_to_global_vertex_index;
1762 }
1763
1764
1765
1766 template <int dim, int spacedim>
1767 void
1768 partition_triangulation(const unsigned int n_partitions,
1770 const SparsityTools::Partitioner partitioner)
1771 {
1773 &triangulation) == nullptr),
1774 ExcMessage("Objects of type parallel::distributed::Triangulation "
1775 "are already partitioned implicitly and can not be "
1776 "partitioned again explicitly."));
1777
1778 std::vector<unsigned int> cell_weights;
1779
1780 // Get cell weighting if a signal has been attached to the triangulation
1781 if (!triangulation.signals.weight.empty())
1782 {
1783 cell_weights.resize(triangulation.n_active_cells(), 0U);
1784
1785 // In a first step, obtain the weights of the locally owned
1786 // cells. For all others, the weight remains at the zero the
1787 // vector was initialized with above.
1788 for (const auto &cell : triangulation.active_cell_iterators())
1789 if (cell->is_locally_owned())
1790 cell_weights[cell->active_cell_index()] =
1792
1793 // If this is a parallel triangulation, we then need to also
1794 // get the weights for all other cells. We have asserted above
1795 // that this function can't be used for
1796 // parallel::distributed::Triangulation objects, so the only
1797 // ones we have to worry about here are
1798 // parallel::shared::Triangulation
1799 if (const auto shared_tria =
1801 &triangulation))
1802 Utilities::MPI::sum(cell_weights,
1803 shared_tria->get_communicator(),
1804 cell_weights);
1805
1806 // verify that the global sum of weights is larger than 0
1807 Assert(std::accumulate(cell_weights.begin(),
1808 cell_weights.end(),
1809 std::uint64_t(0)) > 0,
1810 ExcMessage("The global sum of weights over all active cells "
1811 "is zero. Please verify how you generate weights."));
1812 }
1813
1814 // Call the other more general function
1815 partition_triangulation(n_partitions,
1816 cell_weights,
1818 partitioner);
1819 }
1820
1821
1822
1823 template <int dim, int spacedim>
1824 void
1825 partition_triangulation(const unsigned int n_partitions,
1826 const std::vector<unsigned int> &cell_weights,
1828 const SparsityTools::Partitioner partitioner)
1829 {
1831 &triangulation) == nullptr),
1832 ExcMessage("Objects of type parallel::distributed::Triangulation "
1833 "are already partitioned implicitly and can not be "
1834 "partitioned again explicitly."));
1835 Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
1836
1837 // check for an easy return
1838 if (n_partitions == 1)
1839 {
1840 for (const auto &cell : triangulation.active_cell_iterators())
1841 cell->set_subdomain_id(0);
1842 return;
1843 }
1844
1845 // we decompose the domain by first
1846 // generating the connection graph of all
1847 // cells with their neighbors, and then
1848 // passing this graph off to METIS.
1849 // finally defer to the other function for
1850 // partitioning and assigning subdomain ids
1851 DynamicSparsityPattern cell_connectivity;
1853
1854 SparsityPattern sp_cell_connectivity;
1855 sp_cell_connectivity.copy_from(cell_connectivity);
1856 partition_triangulation(n_partitions,
1857 cell_weights,
1858 sp_cell_connectivity,
1860 partitioner);
1861 }
1862
1863
1864
1865 template <int dim, int spacedim>
1866 void
1867 partition_triangulation(const unsigned int n_partitions,
1868 const SparsityPattern &cell_connection_graph,
1870 const SparsityTools::Partitioner partitioner)
1871 {
1873 &triangulation) == nullptr),
1874 ExcMessage("Objects of type parallel::distributed::Triangulation "
1875 "are already partitioned implicitly and can not be "
1876 "partitioned again explicitly."));
1877
1878 std::vector<unsigned int> cell_weights;
1879
1880 // Get cell weighting if a signal has been attached to the triangulation
1881 if (!triangulation.signals.weight.empty())
1882 {
1883 cell_weights.resize(triangulation.n_active_cells(), 0U);
1884
1885 // In a first step, obtain the weights of the locally owned
1886 // cells. For all others, the weight remains at the zero the
1887 // vector was initialized with above.
1888 for (const auto &cell : triangulation.active_cell_iterators() |
1890 cell_weights[cell->active_cell_index()] =
1892
1893 // If this is a parallel triangulation, we then need to also
1894 // get the weights for all other cells. We have asserted above
1895 // that this function can't be used for
1896 // parallel::distribute::Triangulation objects, so the only
1897 // ones we have to worry about here are
1898 // parallel::shared::Triangulation
1899 if (const auto shared_tria =
1901 &triangulation))
1902 Utilities::MPI::sum(cell_weights,
1903 shared_tria->get_communicator(),
1904 cell_weights);
1905
1906 // verify that the global sum of weights is larger than 0
1907 Assert(std::accumulate(cell_weights.begin(),
1908 cell_weights.end(),
1909 std::uint64_t(0)) > 0,
1910 ExcMessage("The global sum of weights over all active cells "
1911 "is zero. Please verify how you generate weights."));
1912 }
1913
1914 // Call the other more general function
1915 partition_triangulation(n_partitions,
1916 cell_weights,
1917 cell_connection_graph,
1919 partitioner);
1920 }
1921
1922
1923
1924 template <int dim, int spacedim>
1925 void
1926 partition_triangulation(const unsigned int n_partitions,
1927 const std::vector<unsigned int> &cell_weights,
1928 const SparsityPattern &cell_connection_graph,
1930 const SparsityTools::Partitioner partitioner)
1931 {
1933 &triangulation) == nullptr),
1934 ExcMessage("Objects of type parallel::distributed::Triangulation "
1935 "are already partitioned implicitly and can not be "
1936 "partitioned again explicitly."));
1937 Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
1938 Assert(cell_connection_graph.n_rows() == triangulation.n_active_cells(),
1939 ExcMessage("Connectivity graph has wrong size"));
1940 Assert(cell_connection_graph.n_cols() == triangulation.n_active_cells(),
1941 ExcMessage("Connectivity graph has wrong size"));
1942
1943 // signal that partitioning is going to happen
1944 triangulation.signals.pre_partition();
1945
1946 // check for an easy return
1947 if (n_partitions == 1)
1948 {
1949 for (const auto &cell : triangulation.active_cell_iterators())
1950 cell->set_subdomain_id(0);
1951 return;
1952 }
1953
1954 // partition this connection graph and get
1955 // back a vector of indices, one per degree
1956 // of freedom (which is associated with a
1957 // cell)
1958 std::vector<unsigned int> partition_indices(triangulation.n_active_cells());
1959 SparsityTools::partition(cell_connection_graph,
1960 cell_weights,
1961 n_partitions,
1962 partition_indices,
1963 partitioner);
1964
1965 // finally loop over all cells and set the subdomain ids
1966 for (const auto &cell : triangulation.active_cell_iterators())
1967 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
1968 }
1969
1970
1971 namespace internal
1972 {
1976 template <class IT>
1977 void
1979 unsigned int &current_proc_idx,
1980 unsigned int &current_cell_idx,
1981 const unsigned int n_active_cells,
1982 const unsigned int n_partitions)
1983 {
1984 if (cell->is_active())
1985 {
1986 while (current_cell_idx >=
1987 std::floor(static_cast<uint_least64_t>(n_active_cells) *
1988 (current_proc_idx + 1) / n_partitions))
1989 ++current_proc_idx;
1990 cell->set_subdomain_id(current_proc_idx);
1991 ++current_cell_idx;
1992 }
1993 else
1994 {
1995 for (unsigned int n = 0; n < cell->n_children(); ++n)
1997 current_proc_idx,
1998 current_cell_idx,
1999 n_active_cells,
2000 n_partitions);
2001 }
2002 }
2003 } // namespace internal
2004
2005 template <int dim, int spacedim>
2006 void
2007 partition_triangulation_zorder(const unsigned int n_partitions,
2009 const bool group_siblings)
2010 {
2012 &triangulation) == nullptr),
2013 ExcMessage("Objects of type parallel::distributed::Triangulation "
2014 "are already partitioned implicitly and can not be "
2015 "partitioned again explicitly."));
2016 Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
2017 Assert(triangulation.signals.weight.empty(), ExcNotImplemented());
2018
2019 // signal that partitioning is going to happen
2020 triangulation.signals.pre_partition();
2021
2022 // check for an easy return
2023 if (n_partitions == 1)
2024 {
2025 for (const auto &cell : triangulation.active_cell_iterators())
2026 cell->set_subdomain_id(0);
2027 return;
2028 }
2029
2030 // Duplicate the coarse cell reordoring
2031 // as done in p4est
2032 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
2033 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
2034
2035 DynamicSparsityPattern cell_connectivity;
2037 0,
2038 cell_connectivity);
2039 coarse_cell_to_p4est_tree_permutation.resize(triangulation.n_cells(0));
2040 SparsityTools::reorder_hierarchical(cell_connectivity,
2041 coarse_cell_to_p4est_tree_permutation);
2042
2043 p4est_tree_to_coarse_cell_permutation =
2044 Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
2045
2046 unsigned int current_proc_idx = 0;
2047 unsigned int current_cell_idx = 0;
2048 const unsigned int n_active_cells = triangulation.n_active_cells();
2049
2050 // set subdomain id for active cell descendants
2051 // of each coarse cell in permuted order
2052 for (unsigned int idx = 0; idx < triangulation.n_cells(0); ++idx)
2053 {
2054 const unsigned int coarse_cell_idx =
2055 p4est_tree_to_coarse_cell_permutation[idx];
2057 &triangulation, 0, coarse_cell_idx);
2058
2060 current_proc_idx,
2061 current_cell_idx,
2062 n_active_cells,
2063 n_partitions);
2064 }
2065
2066 // if all children of a cell are active (e.g. we
2067 // have a cell that is refined once and no part
2068 // is refined further), p4est places all of them
2069 // on the same processor. The new owner will be
2070 // the processor with the largest number of children
2071 // (ties are broken by picking the lower rank).
2072 // Duplicate this logic here.
2073 if (group_siblings)
2074 {
2076 cell = triangulation.begin(),
2077 endc = triangulation.end();
2078 for (; cell != endc; ++cell)
2079 {
2080 if (cell->is_active())
2081 continue;
2082 bool all_children_active = true;
2083 std::map<unsigned int, unsigned int> map_cpu_n_cells;
2084 for (unsigned int n = 0; n < cell->n_children(); ++n)
2085 if (!cell->child(n)->is_active())
2086 {
2087 all_children_active = false;
2088 break;
2089 }
2090 else
2091 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
2092
2093 if (!all_children_active)
2094 continue;
2095
2096 unsigned int new_owner = cell->child(0)->subdomain_id();
2097 for (std::map<unsigned int, unsigned int>::iterator it =
2098 map_cpu_n_cells.begin();
2099 it != map_cpu_n_cells.end();
2100 ++it)
2101 if (it->second > map_cpu_n_cells[new_owner])
2102 new_owner = it->first;
2103
2104 for (unsigned int n = 0; n < cell->n_children(); ++n)
2105 cell->child(n)->set_subdomain_id(new_owner);
2106 }
2107 }
2108 }
2109
2110
2111 template <int dim, int spacedim>
2112 void
2114 {
2115 unsigned int n_levels = triangulation.n_levels();
2116 for (int lvl = n_levels - 1; lvl >= 0; --lvl)
2117 {
2118 for (const auto &cell : triangulation.cell_iterators_on_level(lvl))
2119 {
2120 if (cell->is_active())
2121 cell->set_level_subdomain_id(cell->subdomain_id());
2122 else
2123 {
2124 Assert(cell->child(0)->level_subdomain_id() !=
2127 cell->set_level_subdomain_id(
2128 cell->child(0)->level_subdomain_id());
2129 }
2130 }
2131 }
2132 }
2133
2134 namespace internal
2135 {
2136 namespace
2137 {
2138 // Split get_subdomain_association() for p::d::T since we want to compile
2139 // it in 1d but none of the p4est stuff is available in 1d.
2140 template <int dim, int spacedim>
2141 void
2145 const std::vector<CellId> &cell_ids,
2146 std::vector<types::subdomain_id> &subdomain_ids)
2147 {
2148#ifndef DEAL_II_WITH_P4EST
2149 (void)triangulation;
2150 (void)cell_ids;
2151 (void)subdomain_ids;
2152 Assert(
2153 false,
2154 ExcMessage(
2155 "You are attempting to use a functionality that is only available "
2156 "if deal.II was configured to use p4est, but cmake did not find a "
2157 "valid p4est library."));
2158#else
2159 // for parallel distributed triangulations, we will ask the p4est oracle
2160 // about the global partitioning of active cells since this information
2161 // is stored on every process
2162 for (const auto &cell_id : cell_ids)
2163 {
2164 // find descendent from coarse quadrant
2165 typename ::internal::p4est::types<dim>::quadrant p4est_cell,
2167
2168 ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
2169 for (const auto &child_index : cell_id.get_child_indices())
2170 {
2171 ::internal::p4est::init_quadrant_children<dim>(
2172 p4est_cell, p4est_children);
2173 p4est_cell =
2174 p4est_children[static_cast<unsigned int>(child_index)];
2175 }
2176
2177 // find owning process, i.e., the subdomain id
2178 const int owner =
2180 const_cast<typename ::internal::p4est::types<dim>::forest
2181 *>(triangulation.get_p4est()),
2182 cell_id.get_coarse_cell_id(),
2183 &p4est_cell,
2186
2187 Assert(owner >= 0, ExcMessage("p4est should know the owner."));
2188
2189 subdomain_ids.push_back(owner);
2190 }
2191#endif
2192 }
2193
2194
2195
2196 template <int spacedim>
2197 void
2200 const std::vector<CellId> &,
2201 std::vector<types::subdomain_id> &)
2202 {
2204 }
2205 } // anonymous namespace
2206 } // namespace internal
2207
2208
2209
2210 template <int dim, int spacedim>
2211 std::vector<types::subdomain_id>
2213 const std::vector<CellId> &cell_ids)
2214 {
2215 std::vector<types::subdomain_id> subdomain_ids;
2216 subdomain_ids.reserve(cell_ids.size());
2217
2218 if (dynamic_cast<
2220 &triangulation) != nullptr)
2221 {
2223 }
2225 *parallel_tria = dynamic_cast<
2227 &triangulation))
2228 {
2229 internal::get_subdomain_association(*parallel_tria,
2230 cell_ids,
2231 subdomain_ids);
2232 }
2233 else if (const parallel::shared::Triangulation<dim, spacedim> *shared_tria =
2235 *>(&triangulation))
2236 {
2237 // for parallel shared triangulations, we need to access true subdomain
2238 // ids which are also valid for artificial cells
2239 const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
2240 shared_tria->get_true_subdomain_ids_of_cells();
2241
2242 for (const auto &cell_id : cell_ids)
2243 {
2244 const unsigned int active_cell_index =
2245 shared_tria->create_cell_iterator(cell_id)->active_cell_index();
2246 subdomain_ids.push_back(
2247 true_subdomain_ids_of_cells[active_cell_index]);
2248 }
2249 }
2250 else
2251 {
2252 // the most general type of triangulation is the serial one. here, all
2253 // subdomain information is directly available
2254 for (const auto &cell_id : cell_ids)
2255 {
2256 subdomain_ids.push_back(
2257 triangulation.create_cell_iterator(cell_id)->subdomain_id());
2258 }
2259 }
2260
2261 return subdomain_ids;
2262 }
2263
2264
2265
2266 template <int dim, int spacedim>
2267 void
2269 std::vector<types::subdomain_id> &subdomain)
2270 {
2271 Assert(subdomain.size() == triangulation.n_active_cells(),
2272 ExcDimensionMismatch(subdomain.size(),
2274 for (const auto &cell : triangulation.active_cell_iterators())
2275 subdomain[cell->active_cell_index()] = cell->subdomain_id();
2276 }
2277
2278
2279
2280 template <int dim, int spacedim>
2281 unsigned int
2284 const types::subdomain_id subdomain)
2285 {
2286 unsigned int count = 0;
2287 for (const auto &cell : triangulation.active_cell_iterators())
2288 if (cell->subdomain_id() == subdomain)
2289 ++count;
2290
2291 return count;
2292 }
2293
2294
2295
2296 template <int dim, int spacedim>
2297 std::vector<bool>
2299 {
2300 // start with all vertices
2301 std::vector<bool> locally_owned_vertices =
2303
2304 // if the triangulation is distributed, eliminate those that
2305 // are owned by other processors -- either because the vertex is
2306 // on an artificial cell, or because it is on a ghost cell with
2307 // a smaller subdomain
2308 if (const auto *tr = dynamic_cast<
2310 &triangulation))
2311 for (const auto &cell : triangulation.active_cell_iterators())
2312 if (cell->is_artificial() ||
2313 (cell->is_ghost() &&
2314 (cell->subdomain_id() < tr->locally_owned_subdomain())))
2315 for (const unsigned int v : cell->vertex_indices())
2316 locally_owned_vertices[cell->vertex_index(v)] = false;
2317
2318 return locally_owned_vertices;
2319 }
2320
2321
2322
2323 namespace internal
2324 {
2325 namespace FixUpDistortedChildCells
2326 {
2327 // compute the mean square
2328 // deviation of the alternating
2329 // forms of the children of the
2330 // given object from that of
2331 // the object itself. for
2332 // objects with
2333 // structdim==spacedim, the
2334 // alternating form is the
2335 // determinant of the jacobian,
2336 // whereas for faces with
2337 // structdim==spacedim-1, the
2338 // alternating form is the
2339 // (signed and scaled) normal
2340 // vector
2341 //
2342 // this average square
2343 // deviation is computed for an
2344 // object where the center node
2345 // has been replaced by the
2346 // second argument to this
2347 // function
2348 template <typename Iterator, int spacedim>
2349 double
2350 objective_function(const Iterator &object,
2351 const Point<spacedim> &object_mid_point)
2352 {
2353 const unsigned int structdim =
2354 Iterator::AccessorType::structure_dimension;
2355 Assert(spacedim == Iterator::AccessorType::dimension,
2357
2358 // everything below is wrong
2359 // if not for the following
2360 // condition
2361 Assert(object->refinement_case() ==
2364 // first calculate the
2365 // average alternating form
2366 // for the parent cell/face
2369 Tensor<spacedim - structdim, spacedim>
2370 parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
2371
2372 for (const unsigned int i : object->vertex_indices())
2373 parent_vertices[i] = object->vertex(i);
2374
2376 parent_vertices, parent_alternating_forms);
2377
2378 const Tensor<spacedim - structdim, spacedim>
2379 average_parent_alternating_form =
2380 std::accumulate(parent_alternating_forms,
2381 parent_alternating_forms +
2384
2385 // now do the same
2386 // computation for the
2387 // children where we use the
2388 // given location for the
2389 // object mid point instead of
2390 // the one the triangulation
2391 // currently reports
2395 Tensor<spacedim - structdim, spacedim> child_alternating_forms
2398
2399 for (unsigned int c = 0; c < object->n_children(); ++c)
2400 for (const unsigned int i : object->child(c)->vertex_indices())
2401 child_vertices[c][i] = object->child(c)->vertex(i);
2402
2403 // replace mid-object
2404 // vertex. note that for
2405 // child i, the mid-object
2406 // vertex happens to have the
2407 // number
2408 // max_children_per_cell-i
2409 for (unsigned int c = 0; c < object->n_children(); ++c)
2411 1] = object_mid_point;
2412
2413 for (unsigned int c = 0; c < object->n_children(); ++c)
2415 child_vertices[c], child_alternating_forms[c]);
2416
2417 // on a uniformly refined
2418 // hypercube object, the child
2419 // alternating forms should
2420 // all be smaller by a factor
2421 // of 2^structdim than the
2422 // ones of the parent. as a
2423 // consequence, we'll use the
2424 // squared deviation from
2425 // this ideal value as an
2426 // objective function
2427 double objective = 0;
2428 for (unsigned int c = 0; c < object->n_children(); ++c)
2429 for (const unsigned int i : object->child(c)->vertex_indices())
2430 objective += (child_alternating_forms[c][i] -
2431 average_parent_alternating_form /
2432 Utilities::fixed_power<structdim>(2))
2433 .norm_square();
2434
2435 return objective;
2436 }
2437
2438
2444 template <typename Iterator>
2446 get_face_midpoint(const Iterator &object,
2447 const unsigned int f,
2448 std::integral_constant<int, 1>)
2449 {
2450 return object->vertex(f);
2451 }
2452
2453
2454
2460 template <typename Iterator>
2462 get_face_midpoint(const Iterator &object,
2463 const unsigned int f,
2464 std::integral_constant<int, 2>)
2465 {
2466 return object->line(f)->center();
2467 }
2468
2469
2470
2476 template <typename Iterator>
2478 get_face_midpoint(const Iterator &object,
2479 const unsigned int f,
2480 std::integral_constant<int, 3>)
2481 {
2482 return object->face(f)->center();
2483 }
2484
2485
2486
2509 template <typename Iterator>
2510 double
2511 minimal_diameter(const Iterator &object)
2512 {
2513 const unsigned int structdim =
2514 Iterator::AccessorType::structure_dimension;
2515
2516 double diameter = object->diameter();
2517 for (const unsigned int f : object->face_indices())
2518 for (unsigned int e = f + 1; e < object->n_faces(); ++e)
2520 diameter,
2521 get_face_midpoint(object,
2522 f,
2523 std::integral_constant<int, structdim>())
2524 .distance(get_face_midpoint(
2525 object, e, std::integral_constant<int, structdim>())));
2526
2527 return diameter;
2528 }
2529
2530
2531
2536 template <typename Iterator>
2537 bool
2538 fix_up_object(const Iterator &object)
2539 {
2540 const unsigned int structdim =
2541 Iterator::AccessorType::structure_dimension;
2542 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
2543
2544 // right now we can only deal with cells that have been refined
2545 // isotropically because that is the only case where we have a cell
2546 // mid-point that can be moved around without having to consider
2547 // boundary information
2548 Assert(object->has_children(), ExcInternalError());
2549 Assert(object->refinement_case() ==
2552
2553 // get the current location of the object mid-vertex:
2554 Point<spacedim> object_mid_point = object->child(0)->vertex(
2556
2557 // now do a few steepest descent steps to reduce the objective
2558 // function. compute the diameter in the helper function above
2559 unsigned int iteration = 0;
2560 const double diameter = minimal_diameter(object);
2561
2562 // current value of objective function and initial delta
2563 double current_value = objective_function(object, object_mid_point);
2564 double initial_delta = 0;
2565
2566 do
2567 {
2568 // choose a step length that is initially 1/4 of the child
2569 // objects' diameter, and a sequence whose sum does not converge
2570 // (to avoid premature termination of the iteration)
2571 const double step_length = diameter / 4 / (iteration + 1);
2572
2573 // compute the objective function's derivative using a two-sided
2574 // difference formula with eps=step_length/10
2575 Tensor<1, spacedim> gradient;
2576 for (unsigned int d = 0; d < spacedim; ++d)
2577 {
2578 const double eps = step_length / 10;
2579
2581 h[d] = eps / 2;
2582
2583 gradient[d] =
2585 object, project_to_object(object, object_mid_point + h)) -
2587 object, project_to_object(object, object_mid_point - h))) /
2588 eps;
2589 }
2590
2591 // there is nowhere to go
2592 if (gradient.norm() == 0)
2593 break;
2594
2595 // We need to go in direction -gradient. the optimal value of the
2596 // objective function is zero, so assuming that the model is
2597 // quadratic we would have to go -2*val/||gradient|| in this
2598 // direction, make sure we go at most step_length into this
2599 // direction
2600 object_mid_point -=
2601 std::min(2 * current_value / (gradient * gradient),
2602 step_length / gradient.norm()) *
2603 gradient;
2604 object_mid_point = project_to_object(object, object_mid_point);
2605
2606 // compute current value of the objective function
2607 const double previous_value = current_value;
2608 current_value = objective_function(object, object_mid_point);
2609
2610 if (iteration == 0)
2611 initial_delta = (previous_value - current_value);
2612
2613 // stop if we aren't moving much any more
2614 if ((iteration >= 1) &&
2615 ((previous_value - current_value < 0) ||
2616 (std::fabs(previous_value - current_value) <
2617 0.001 * initial_delta)))
2618 break;
2619
2620 ++iteration;
2621 }
2622 while (iteration < 20);
2623
2624 // verify that the new
2625 // location is indeed better
2626 // than the one before. check
2627 // this by comparing whether
2628 // the minimum value of the
2629 // products of parent and
2630 // child alternating forms is
2631 // positive. for cells this
2632 // means that the
2633 // determinants have the same
2634 // sign, for faces that the
2635 // face normals of parent and
2636 // children point in the same
2637 // general direction
2638 double old_min_product, new_min_product;
2639
2642 for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
2643 parent_vertices[i] = object->vertex(i);
2644
2645 Tensor<spacedim - structdim, spacedim>
2646 parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
2648 parent_vertices, parent_alternating_forms);
2649
2653
2654 for (unsigned int c = 0; c < object->n_children(); ++c)
2655 for (const unsigned int i : object->child(c)->vertex_indices())
2656 child_vertices[c][i] = object->child(c)->vertex(i);
2657
2658 Tensor<spacedim - structdim, spacedim> child_alternating_forms
2661
2662 for (unsigned int c = 0; c < object->n_children(); ++c)
2664 child_vertices[c], child_alternating_forms[c]);
2665
2666 old_min_product =
2667 child_alternating_forms[0][0] * parent_alternating_forms[0];
2668 for (unsigned int c = 0; c < object->n_children(); ++c)
2669 for (const unsigned int i : object->child(c)->vertex_indices())
2670 for (const unsigned int j : object->vertex_indices())
2671 old_min_product = std::min<double>(old_min_product,
2672 child_alternating_forms[c][i] *
2673 parent_alternating_forms[j]);
2674
2675 // for the new minimum value,
2676 // replace mid-object
2677 // vertex. note that for child
2678 // i, the mid-object vertex
2679 // happens to have the number
2680 // max_children_per_cell-i
2681 for (unsigned int c = 0; c < object->n_children(); ++c)
2683 1] = object_mid_point;
2684
2685 for (unsigned int c = 0; c < object->n_children(); ++c)
2687 child_vertices[c], child_alternating_forms[c]);
2688
2689 new_min_product =
2690 child_alternating_forms[0][0] * parent_alternating_forms[0];
2691 for (unsigned int c = 0; c < object->n_children(); ++c)
2692 for (const unsigned int i : object->child(c)->vertex_indices())
2693 for (const unsigned int j : object->vertex_indices())
2694 new_min_product = std::min<double>(new_min_product,
2695 child_alternating_forms[c][i] *
2696 parent_alternating_forms[j]);
2697
2698 // if new minimum value is
2699 // better than before, then set the
2700 // new mid point. otherwise
2701 // return this object as one of
2702 // those that can't apparently
2703 // be fixed
2704 if (new_min_product >= old_min_product)
2705 object->child(0)->vertex(
2707 object_mid_point;
2708
2709 // return whether after this
2710 // operation we have an object that
2711 // is well oriented
2712 return (std::max(new_min_product, old_min_product) > 0);
2713 }
2714
2715
2716
2717 // possibly fix up the faces of a cell by moving around its mid-points
2718 template <int dim, int spacedim>
2719 void
2721 const typename ::Triangulation<dim, spacedim>::cell_iterator
2722 &cell,
2723 std::integral_constant<int, dim>,
2724 std::integral_constant<int, spacedim>)
2725 {
2726 // see if we first can fix up some of the faces of this object. We can
2727 // mess with faces if and only if the neighboring cell is not even
2728 // more refined than we are (since in that case the sub-faces have
2729 // themselves children that we can't move around any more). however,
2730 // the latter case shouldn't happen anyway: if the current face is
2731 // distorted but the neighbor is even more refined, then the face had
2732 // been deformed before already, and had been ignored at the time; we
2733 // should then also be able to ignore it this time as well
2734 for (auto f : cell->face_indices())
2735 {
2736 Assert(cell->face(f)->has_children(), ExcInternalError());
2737 Assert(cell->face(f)->refinement_case() ==
2740
2741 bool subface_is_more_refined = false;
2742 for (unsigned int g = 0;
2743 g < GeometryInfo<dim>::max_children_per_face;
2744 ++g)
2745 if (cell->face(f)->child(g)->has_children())
2746 {
2747 subface_is_more_refined = true;
2748 break;
2749 }
2750
2751 if (subface_is_more_refined == true)
2752 continue;
2753
2754 // we finally know that we can do something about this face
2755 fix_up_object(cell->face(f));
2756 }
2757 }
2758 } /* namespace FixUpDistortedChildCells */
2759 } /* namespace internal */
2760
2761
2762 template <int dim, int spacedim>
2766 &distorted_cells,
2767 Triangulation<dim, spacedim> & /*triangulation*/)
2768 {
2769 static_assert(
2770 dim != 1 && spacedim != 1,
2771 "This function is only valid when dim != 1 or spacedim != 1.");
2772 typename Triangulation<dim, spacedim>::DistortedCellList unfixable_subset;
2773
2774 // loop over all cells that we have to fix up
2775 for (typename std::list<
2776 typename Triangulation<dim, spacedim>::cell_iterator>::const_iterator
2777 cell_ptr = distorted_cells.distorted_cells.begin();
2778 cell_ptr != distorted_cells.distorted_cells.end();
2779 ++cell_ptr)
2780 {
2781 const typename Triangulation<dim, spacedim>::cell_iterator &cell =
2782 *cell_ptr;
2783
2784 Assert(!cell->is_active(),
2785 ExcMessage(
2786 "This function is only valid for a list of cells that "
2787 "have children (i.e., no cell in the list may be active)."));
2788
2790 cell,
2791 std::integral_constant<int, dim>(),
2792 std::integral_constant<int, spacedim>());
2793
2794 // If possible, fix up the object.
2796 unfixable_subset.distorted_cells.push_back(cell);
2797 }
2798
2799 return unfixable_subset;
2800 }
2801
2802
2803
2804 template <int dim, int spacedim>
2805 void
2807 const bool reset_boundary_ids)
2808 {
2809 const auto src_boundary_ids = tria.get_boundary_ids();
2810 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
2811 auto m_it = dst_manifold_ids.begin();
2812 for (const auto b : src_boundary_ids)
2813 {
2814 *m_it = static_cast<types::manifold_id>(b);
2815 ++m_it;
2816 }
2817 const std::vector<types::boundary_id> reset_boundary_id =
2818 reset_boundary_ids ?
2819 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
2820 src_boundary_ids;
2821 map_boundary_to_manifold_ids(src_boundary_ids,
2822 dst_manifold_ids,
2823 tria,
2824 reset_boundary_id);
2825 }
2826
2827
2828
2829 template <int dim, int spacedim>
2830 void
2832 const std::vector<types::boundary_id> &src_boundary_ids,
2833 const std::vector<types::manifold_id> &dst_manifold_ids,
2835 const std::vector<types::boundary_id> &reset_boundary_ids_)
2836 {
2837 AssertDimension(src_boundary_ids.size(), dst_manifold_ids.size());
2838 const auto reset_boundary_ids =
2839 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
2840 AssertDimension(reset_boundary_ids.size(), src_boundary_ids.size());
2841
2842 // in 3d, we not only have to copy boundary ids of faces, but also of edges
2843 // because we see them twice (once from each adjacent boundary face),
2844 // we cannot immediately reset their boundary ids. thus, copy first
2845 // and reset later
2846 if (dim >= 3)
2847 for (const auto &cell : tria.active_cell_iterators())
2848 for (auto f : cell->face_indices())
2849 if (cell->face(f)->at_boundary())
2850 for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
2851 {
2852 const auto bid = cell->face(f)->line(e)->boundary_id();
2853 const unsigned int ind = std::find(src_boundary_ids.begin(),
2854 src_boundary_ids.end(),
2855 bid) -
2856 src_boundary_ids.begin();
2857 if (ind < src_boundary_ids.size())
2858 cell->face(f)->line(e)->set_manifold_id(
2859 dst_manifold_ids[ind]);
2860 }
2861
2862 // now do cells
2863 for (const auto &cell : tria.active_cell_iterators())
2864 for (auto f : cell->face_indices())
2865 if (cell->face(f)->at_boundary())
2866 {
2867 const auto bid = cell->face(f)->boundary_id();
2868 const unsigned int ind =
2869 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
2870 src_boundary_ids.begin();
2871
2872 if (ind < src_boundary_ids.size())
2873 {
2874 // assign the manifold id
2875 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
2876 // then reset boundary id
2877 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
2878 }
2879
2880 if (dim >= 3)
2881 for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
2882 {
2883 const auto bid = cell->face(f)->line(e)->boundary_id();
2884 const unsigned int ind = std::find(src_boundary_ids.begin(),
2885 src_boundary_ids.end(),
2886 bid) -
2887 src_boundary_ids.begin();
2888 if (ind < src_boundary_ids.size())
2889 cell->face(f)->line(e)->set_boundary_id(
2890 reset_boundary_ids[ind]);
2891 }
2892 }
2893 }
2894
2895
2896 template <int dim, int spacedim>
2897 void
2899 const bool compute_face_ids)
2900 {
2902 cell = tria.begin_active(),
2903 endc = tria.end();
2904
2905 for (; cell != endc; ++cell)
2906 {
2907 cell->set_manifold_id(cell->material_id());
2908 if (compute_face_ids == true)
2909 {
2910 for (auto f : cell->face_indices())
2911 {
2912 if (cell->at_boundary(f) == false)
2913 cell->face(f)->set_manifold_id(
2914 std::min(cell->material_id(),
2915 cell->neighbor(f)->material_id()));
2916 else
2917 cell->face(f)->set_manifold_id(cell->material_id());
2918 }
2919 }
2920 }
2921 }
2922
2923
2924 template <int dim, int spacedim>
2925 void
2928 const std::function<types::manifold_id(
2929 const std::set<types::manifold_id> &)> &disambiguation_function,
2930 bool overwrite_only_flat_manifold_ids)
2931 {
2932 // Easy case first:
2933 if (dim == 1)
2934 return;
2935 const unsigned int n_subobjects =
2936 dim == 2 ? tria.n_lines() : tria.n_lines() + tria.n_quads();
2937
2938 // If user index is zero, then it has not been set.
2939 std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
2940 std::vector<unsigned int> backup;
2941 tria.save_user_indices(backup);
2943
2944 unsigned next_index = 1;
2945 for (auto &cell : tria.active_cell_iterators())
2946 {
2947 if (dim > 1)
2948 for (unsigned int l = 0; l < cell->n_lines(); ++l)
2949 {
2950 if (cell->line(l)->user_index() == 0)
2951 {
2952 AssertIndexRange(next_index, n_subobjects + 1);
2953 manifold_ids[next_index].insert(cell->manifold_id());
2954 cell->line(l)->set_user_index(next_index++);
2955 }
2956 else
2957 manifold_ids[cell->line(l)->user_index()].insert(
2958 cell->manifold_id());
2959 }
2960 if (dim > 2)
2961 for (unsigned int l = 0; l < cell->n_faces(); ++l)
2962 {
2963 if (cell->quad(l)->user_index() == 0)
2964 {
2965 AssertIndexRange(next_index, n_subobjects + 1);
2966 manifold_ids[next_index].insert(cell->manifold_id());
2967 cell->quad(l)->set_user_index(next_index++);
2968 }
2969 else
2970 manifold_ids[cell->quad(l)->user_index()].insert(
2971 cell->manifold_id());
2972 }
2973 }
2974 for (auto &cell : tria.active_cell_iterators())
2975 {
2976 if (dim > 1)
2977 for (unsigned int l = 0; l < cell->n_lines(); ++l)
2978 {
2979 const auto id = cell->line(l)->user_index();
2980 // Make sure we change the manifold indicator only once
2981 if (id != 0)
2982 {
2983 if (cell->line(l)->manifold_id() ==
2985 overwrite_only_flat_manifold_ids == false)
2986 cell->line(l)->set_manifold_id(
2987 disambiguation_function(manifold_ids[id]));
2988 cell->line(l)->set_user_index(0);
2989 }
2990 }
2991 if (dim > 2)
2992 for (unsigned int l = 0; l < cell->n_faces(); ++l)
2993 {
2994 const auto id = cell->quad(l)->user_index();
2995 // Make sure we change the manifold indicator only once
2996 if (id != 0)
2997 {
2998 if (cell->quad(l)->manifold_id() ==
3000 overwrite_only_flat_manifold_ids == false)
3001 cell->quad(l)->set_manifold_id(
3002 disambiguation_function(manifold_ids[id]));
3003 cell->quad(l)->set_user_index(0);
3004 }
3005 }
3006 }
3007 tria.load_user_indices(backup);
3008 }
3009
3010
3011
3012 template <int dim, int spacedim>
3013 void
3015 const double limit_angle_fraction)
3016 {
3017 if (dim == 1)
3018 return; // Nothing to do
3019
3020 // Check that we don't have hanging nodes
3022 ExcMessage("The input Triangulation cannot "
3023 "have hanging nodes."));
3024
3026
3027 bool has_cells_with_more_than_dim_faces_on_boundary = true;
3028 bool has_cells_with_dim_faces_on_boundary = false;
3029
3030 unsigned int refinement_cycles = 0;
3031
3032 while (has_cells_with_more_than_dim_faces_on_boundary)
3033 {
3034 has_cells_with_more_than_dim_faces_on_boundary = false;
3035
3036 for (const auto &cell : tria.active_cell_iterators())
3037 {
3038 unsigned int boundary_face_counter = 0;
3039 for (auto f : cell->face_indices())
3040 if (cell->face(f)->at_boundary())
3041 ++boundary_face_counter;
3042 if (boundary_face_counter > dim)
3043 {
3044 has_cells_with_more_than_dim_faces_on_boundary = true;
3045 break;
3046 }
3047 else if (boundary_face_counter == dim)
3048 has_cells_with_dim_faces_on_boundary = true;
3049 }
3050 if (has_cells_with_more_than_dim_faces_on_boundary)
3051 {
3053 ++refinement_cycles;
3054 }
3055 }
3056
3057 if (has_cells_with_dim_faces_on_boundary)
3058 {
3060 ++refinement_cycles;
3061 }
3062 else
3063 {
3064 while (refinement_cycles > 0)
3065 {
3066 for (const auto &cell : tria.active_cell_iterators())
3067 cell->set_coarsen_flag();
3069 refinement_cycles--;
3070 }
3071 return;
3072 }
3073
3074 std::vector<bool> cells_to_remove(tria.n_active_cells(), false);
3075 std::vector<Point<spacedim>> vertices = tria.get_vertices();
3076
3077 std::vector<bool> faces_to_remove(tria.n_raw_faces(), false);
3078
3079 std::vector<CellData<dim>> cells_to_add;
3080 SubCellData subcelldata_to_add;
3081
3082 // Trick compiler for dimension independent things
3083 const unsigned int v0 = 0, v1 = 1, v2 = (dim > 1 ? 2 : 0),
3084 v3 = (dim > 1 ? 3 : 0);
3085
3086 for (const auto &cell : tria.active_cell_iterators())
3087 {
3088 double angle_fraction = 0;
3089 unsigned int vertex_at_corner = numbers::invalid_unsigned_int;
3090
3091 if (dim == 2)
3092 {
3094 p0[spacedim > 1 ? 1 : 0] = 1;
3096 p1[0] = 1;
3097
3098 if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
3099 {
3100 p0 = cell->vertex(v0) - cell->vertex(v2);
3101 p1 = cell->vertex(v3) - cell->vertex(v2);
3102 vertex_at_corner = v2;
3103 }
3104 else if (cell->face(v3)->at_boundary() &&
3105 cell->face(v1)->at_boundary())
3106 {
3107 p0 = cell->vertex(v2) - cell->vertex(v3);
3108 p1 = cell->vertex(v1) - cell->vertex(v3);
3109 vertex_at_corner = v3;
3110 }
3111 else if (cell->face(1)->at_boundary() &&
3112 cell->face(2)->at_boundary())
3113 {
3114 p0 = cell->vertex(v0) - cell->vertex(v1);
3115 p1 = cell->vertex(v3) - cell->vertex(v1);
3116 vertex_at_corner = v1;
3117 }
3118 else if (cell->face(2)->at_boundary() &&
3119 cell->face(0)->at_boundary())
3120 {
3121 p0 = cell->vertex(v2) - cell->vertex(v0);
3122 p1 = cell->vertex(v1) - cell->vertex(v0);
3123 vertex_at_corner = v0;
3124 }
3125 p0 /= p0.norm();
3126 p1 /= p1.norm();
3127 angle_fraction = std::acos(p0 * p1) / numbers::PI;
3128 }
3129 else
3130 {
3132 }
3133
3134 if (angle_fraction > limit_angle_fraction)
3135 {
3136 auto flags_removal = [&](unsigned int f1,
3137 unsigned int f2,
3138 unsigned int n1,
3139 unsigned int n2) -> void {
3140 cells_to_remove[cell->active_cell_index()] = true;
3141 cells_to_remove[cell->neighbor(n1)->active_cell_index()] = true;
3142 cells_to_remove[cell->neighbor(n2)->active_cell_index()] = true;
3143
3144 faces_to_remove[cell->face(f1)->index()] = true;
3145 faces_to_remove[cell->face(f2)->index()] = true;
3146
3147 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] = true;
3148 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] = true;
3149 };
3150
3151 auto cell_creation = [&](const unsigned int vv0,
3152 const unsigned int vv1,
3153 const unsigned int f0,
3154 const unsigned int f1,
3155
3156 const unsigned int n0,
3157 const unsigned int v0n0,
3158 const unsigned int v1n0,
3159
3160 const unsigned int n1,
3161 const unsigned int v0n1,
3162 const unsigned int v1n1) {
3163 CellData<dim> c1, c2;
3164 CellData<1> l1, l2;
3165
3166 c1.vertices[v0] = cell->vertex_index(vv0);
3167 c1.vertices[v1] = cell->vertex_index(vv1);
3168 c1.vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
3169 c1.vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
3170
3171 c1.manifold_id = cell->manifold_id();
3172 c1.material_id = cell->material_id();
3173
3174 c2.vertices[v0] = cell->vertex_index(vv0);
3175 c2.vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
3176 c2.vertices[v2] = cell->vertex_index(vv1);
3177 c2.vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
3178
3179 c2.manifold_id = cell->manifold_id();
3180 c2.material_id = cell->material_id();
3181
3182 l1.vertices[0] = cell->vertex_index(vv0);
3183 l1.vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
3184
3185 l1.boundary_id = cell->line(f0)->boundary_id();
3186 l1.manifold_id = cell->line(f0)->manifold_id();
3187 subcelldata_to_add.boundary_lines.push_back(l1);
3188
3189 l2.vertices[0] = cell->vertex_index(vv0);
3190 l2.vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
3191
3192 l2.boundary_id = cell->line(f1)->boundary_id();
3193 l2.manifold_id = cell->line(f1)->manifold_id();
3194 subcelldata_to_add.boundary_lines.push_back(l2);
3195
3196 cells_to_add.push_back(c1);
3197 cells_to_add.push_back(c2);
3198 };
3199
3200 if (dim == 2)
3201 {
3202 switch (vertex_at_corner)
3203 {
3204 case 0:
3205 flags_removal(0, 2, 3, 1);
3206 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
3207 break;
3208 case 1:
3209 flags_removal(1, 2, 3, 0);
3210 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
3211 break;
3212 case 2:
3213 flags_removal(3, 0, 1, 2);
3214 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
3215 break;
3216 case 3:
3217 flags_removal(3, 1, 0, 2);
3218 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
3219 break;
3220 }
3221 }
3222 else
3223 {
3225 }
3226 }
3227 }
3228
3229 // if no cells need to be added, then no regularization is necessary.
3230 // Restore things as they were before this function was called.
3231 if (cells_to_add.empty())
3232 {
3233 while (refinement_cycles > 0)
3234 {
3235 for (const auto &cell : tria.active_cell_iterators())
3236 cell->set_coarsen_flag();
3238 refinement_cycles--;
3239 }
3240 return;
3241 }
3242
3243 // add the cells that were not marked as skipped
3244 for (const auto &cell : tria.active_cell_iterators())
3245 {
3246 if (cells_to_remove[cell->active_cell_index()] == false)
3247 {
3248 CellData<dim> c(cell->n_vertices());
3249 for (const unsigned int v : cell->vertex_indices())
3250 c.vertices[v] = cell->vertex_index(v);
3251 c.manifold_id = cell->manifold_id();
3252 c.material_id = cell->material_id();
3253 cells_to_add.push_back(c);
3254 }
3255 }
3256
3257 // Face counter for both dim == 2 and dim == 3
3259 face = tria.begin_active_face(),
3260 endf = tria.end_face();
3261 for (; face != endf; ++face)
3262 if ((face->at_boundary() ||
3263 face->manifold_id() != numbers::flat_manifold_id) &&
3264 faces_to_remove[face->index()] == false)
3265 {
3266 for (unsigned int l = 0; l < face->n_lines(); ++l)
3267 {
3268 CellData<1> line;
3269 if (dim == 2)
3270 {
3271 for (const unsigned int v : face->vertex_indices())
3272 line.vertices[v] = face->vertex_index(v);
3273 line.boundary_id = face->boundary_id();
3274 line.manifold_id = face->manifold_id();
3275 }
3276 else
3277 {
3278 for (const unsigned int v : face->line(l)->vertex_indices())
3279 line.vertices[v] = face->line(l)->vertex_index(v);
3280 line.boundary_id = face->line(l)->boundary_id();
3281 line.manifold_id = face->line(l)->manifold_id();
3282 }
3283 subcelldata_to_add.boundary_lines.push_back(line);
3284 }
3285 if (dim == 3)
3286 {
3287 CellData<2> quad(face->n_vertices());
3288 for (const unsigned int v : face->vertex_indices())
3289 quad.vertices[v] = face->vertex_index(v);
3290 quad.boundary_id = face->boundary_id();
3291 quad.manifold_id = face->manifold_id();
3292 subcelldata_to_add.boundary_quads.push_back(quad);
3293 }
3294 }
3296 cells_to_add,
3297 subcelldata_to_add);
3299
3300 // Save manifolds
3301 auto manifold_ids = tria.get_manifold_ids();
3302 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
3303 manifolds;
3304 // Set manifolds in new Triangulation
3305 for (const auto manifold_id : manifold_ids)
3306 if (manifold_id != numbers::flat_manifold_id)
3307 manifolds[manifold_id] = tria.get_manifold(manifold_id).clone();
3308
3309 tria.clear();
3310
3311 tria.create_triangulation(vertices, cells_to_add, subcelldata_to_add);
3312
3313 // Restore manifolds
3314 for (const auto manifold_id : manifold_ids)
3315 if (manifold_id != numbers::flat_manifold_id)
3316 tria.set_manifold(manifold_id, *manifolds[manifold_id]);
3317 }
3318
3319
3320
3321 template <int dim, int spacedim>
3322#ifndef DOXYGEN
3323 std::tuple<
3324 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3325 std::vector<std::vector<Point<dim>>>,
3326 std::vector<std::vector<unsigned int>>>
3327#else
3328 return_type
3329#endif
3331 const Cache<dim, spacedim> &cache,
3332 const std::vector<Point<spacedim>> &points,
3334 &cell_hint)
3335 {
3336 const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
3337 // Splitting the tuple's components
3338 auto &cells = std::get<0>(cqmp);
3339 auto &qpoints = std::get<1>(cqmp);
3340 auto &maps = std::get<2>(cqmp);
3341
3342 return std::make_tuple(std::move(cells),
3343 std::move(qpoints),
3344 std::move(maps));
3345 }
3346
3347
3348
3349 template <int dim, int spacedim>
3350#ifndef DOXYGEN
3351 std::tuple<
3352 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3353 std::vector<std::vector<Point<dim>>>,
3354 std::vector<std::vector<unsigned int>>,
3355 std::vector<unsigned int>>
3356#else
3357 return_type
3358#endif
3360 const Cache<dim, spacedim> &cache,
3361 const std::vector<Point<spacedim>> &points,
3363 &cell_hint)
3364 {
3365 Assert((dim == spacedim),
3366 ExcMessage("Only implemented for dim==spacedim."));
3367
3368 // Alias
3369 namespace bgi = boost::geometry::index;
3370
3371 // Get the mapping
3372 const auto &mapping = cache.get_mapping();
3373
3374 // How many points are here?
3375 const unsigned int np = points.size();
3376
3377 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
3378 cells_out;
3379 std::vector<std::vector<Point<dim>>> qpoints_out;
3380 std::vector<std::vector<unsigned int>> maps_out;
3381 std::vector<unsigned int> missing_points_out;
3382
3383 // Now the easy case.
3384 if (np == 0)
3385 return std::make_tuple(std::move(cells_out),
3386 std::move(qpoints_out),
3387 std::move(maps_out),
3388 std::move(missing_points_out));
3389
3390 // For the search we shall use the following tree
3391 const auto &b_tree = cache.get_cell_bounding_boxes_rtree();
3392
3393 // Now make a tree of indices for the points
3394 // [TODO] This would work better with pack_rtree_of_indices, but
3395 // windows does not like it. Build a tree with pairs of point and id
3396 std::vector<std::pair<Point<spacedim>, unsigned int>> points_and_ids(np);
3397 for (unsigned int i = 0; i < np; ++i)
3398 points_and_ids[i] = std::make_pair(points[i], i);
3399 const auto p_tree = pack_rtree(points_and_ids);
3400
3401 // Keep track of all found points
3402 std::vector<bool> found_points(points.size(), false);
3403
3404 // Check if a point was found
3405 const auto already_found = [&found_points](const auto &id) {
3406 AssertIndexRange(id.second, found_points.size());
3407 return found_points[id.second];
3408 };
3409
3410 // check if the given cell was already in the vector of cells before. If so,
3411 // insert in the corresponding vectors the reference point and the id.
3412 // Otherwise append a new entry to all vectors.
3413 const auto store_cell_point_and_id =
3414 [&](
3416 const Point<dim> &ref_point,
3417 const unsigned int &id) {
3418 const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
3419 if (it != cells_out.rend())
3420 {
3421 const auto cell_id =
3422 (cells_out.size() - 1 - (it - cells_out.rbegin()));
3423 qpoints_out[cell_id].emplace_back(ref_point);
3424 maps_out[cell_id].emplace_back(id);
3425 }
3426 else
3427 {
3428 cells_out.emplace_back(cell);
3429 qpoints_out.emplace_back(std::vector<Point<dim>>({ref_point}));
3430 maps_out.emplace_back(std::vector<unsigned int>({id}));
3431 }
3432 };
3433
3434 // Check all points within a given pair of box and cell
3435 const auto check_all_points_within_box = [&](const auto &leaf) {
3436 const double relative_tolerance = 1e-12;
3437 const BoundingBox<spacedim> box =
3438 leaf.first.create_extended_relative(relative_tolerance);
3439 const auto &cell_hint = leaf.second;
3440
3441 for (const auto &point_and_id :
3442 p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
3443 bgi::intersects(box)))
3444 {
3445 const auto id = point_and_id.second;
3446 const auto cell_and_ref =
3448 points[id],
3449 cell_hint);
3450 const auto &cell = cell_and_ref.first;
3451 const auto &ref_point = cell_and_ref.second;
3452
3453 if (cell.state() == IteratorState::valid)
3454 store_cell_point_and_id(cell, ref_point, id);
3455 else
3456 missing_points_out.emplace_back(id);
3457
3458 // Don't look anymore for this point
3459 found_points[id] = true;
3460 }
3461 };
3462
3463 // If a hint cell was given, use it
3464 if (cell_hint.state() == IteratorState::valid)
3465 check_all_points_within_box(
3466 std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
3467
3468 // Now loop over all points that have not been found yet
3469 for (unsigned int i = 0; i < np; ++i)
3470 if (found_points[i] == false)
3471 {
3472 // Get the closest cell to this point
3473 const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
3474 // Now checks all points that fall within this box
3475 if (leaf != b_tree.qend())
3476 check_all_points_within_box(*leaf);
3477 else
3478 {
3479 // We should not get here. Throw an error.
3481 }
3482 }
3483 // Now make sure we send out the rest of the points that we did not find.
3484 for (unsigned int i = 0; i < np; ++i)
3485 if (found_points[i] == false)
3486 missing_points_out.emplace_back(i);
3487
3488 // Debug Checking
3489 AssertDimension(cells_out.size(), maps_out.size());
3490 AssertDimension(cells_out.size(), qpoints_out.size());
3491
3492#ifdef DEBUG
3493 unsigned int c = cells_out.size();
3494 unsigned int qps = 0;
3495 // The number of points in all
3496 // the cells must be the same as
3497 // the number of points we
3498 // started off from,
3499 // plus the points which were ignored
3500 for (unsigned int n = 0; n < c; ++n)
3501 {
3502 AssertDimension(qpoints_out[n].size(), maps_out[n].size());
3503 qps += qpoints_out[n].size();
3504 }
3505
3506 Assert(qps + missing_points_out.size() == np,
3507 ExcDimensionMismatch(qps + missing_points_out.size(), np));
3508#endif
3509
3510 return std::make_tuple(std::move(cells_out),
3511 std::move(qpoints_out),
3512 std::move(maps_out),
3513 std::move(missing_points_out));
3514 }
3515
3516
3517
3518 template <int dim, int spacedim>
3519#ifndef DOXYGEN
3520 std::tuple<
3521 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3522 std::vector<std::vector<Point<dim>>>,
3523 std::vector<std::vector<unsigned int>>,
3524 std::vector<std::vector<Point<spacedim>>>,
3525 std::vector<std::vector<unsigned int>>>
3526#else
3527 return_type
3528#endif
3531 const std::vector<Point<spacedim>> &points,
3532 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3533 const double tolerance,
3534 const std::vector<bool> &marked_vertices,
3535 const bool enforce_unique_mapping)
3536 {
3537 // run internal function ...
3538 const auto all =
3540 points,
3541 global_bboxes,
3543 tolerance,
3544 false,
3545 enforce_unique_mapping)
3546 .send_components;
3547
3548 // ... and reshuffle the data
3549 std::tuple<
3550 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3551 std::vector<std::vector<Point<dim>>>,
3552 std::vector<std::vector<unsigned int>>,
3553 std::vector<std::vector<Point<spacedim>>>,
3554 std::vector<std::vector<unsigned int>>>
3555 result;
3556
3557 std::pair<int, int> dummy{-1, -1};
3558
3559 for (unsigned int i = 0; i < all.size(); ++i)
3560 {
3561 if (dummy != std::get<0>(all[i]))
3562 {
3563 std::get<0>(result).push_back(
3565 &cache.get_triangulation(),
3566 std::get<0>(all[i]).first,
3567 std::get<0>(all[i]).second});
3568
3569 const unsigned int new_size = std::get<0>(result).size();
3570
3571 std::get<1>(result).resize(new_size);
3572 std::get<2>(result).resize(new_size);
3573 std::get<3>(result).resize(new_size);
3574 std::get<4>(result).resize(new_size);
3575
3576 dummy = std::get<0>(all[i]);
3577 }
3578
3579 std::get<1>(result).back().push_back(
3580 std::get<3>(all[i])); // reference point
3581 std::get<2>(result).back().push_back(std::get<2>(all[i])); // index
3582 std::get<3>(result).back().push_back(std::get<4>(all[i])); // real point
3583 std::get<4>(result).back().push_back(std::get<1>(all[i])); // rank
3584 }
3585
3586 return result;
3587 }
3588
3589
3590
3591 namespace internal
3592 {
3599 template <int spacedim, typename T>
3600 std::tuple<std::vector<unsigned int>,
3601 std::vector<unsigned int>,
3602 std::vector<unsigned int>>
3604 const MPI_Comm comm,
3605 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3606 const std::vector<T> &entities,
3607 const double tolerance)
3608 {
3609 std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
3610 ranks_and_indices.reserve(entities.size());
3611
3612#if defined(DEAL_II_WITH_ARBORX)
3613 static constexpr bool use_arborx = true;
3614#else
3615 static constexpr bool use_arborx = false;
3616#endif
3617 // Lambda to process bboxes if global_bboxes.size()>1 or ArborX not
3618 // available
3619 const auto process_bboxes = [&]() -> void {
3620 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes_temp;
3621 auto *global_bboxes_to_be_used = &global_bboxes;
3622
3623 if (global_bboxes.size() == 1 && use_arborx == false)
3624 {
3625 global_bboxes_temp =
3626 Utilities::MPI::all_gather(comm, global_bboxes[0]);
3627 global_bboxes_to_be_used = &global_bboxes_temp;
3628 }
3629
3630 // helper function to determine if a bounding box is valid
3631 const auto is_valid = [](const auto &bb) {
3632 for (unsigned int i = 0; i < spacedim; ++i)
3633 if (bb.get_boundary_points().first[i] >
3634 bb.get_boundary_points().second[i])
3635 return false;
3636
3637 return true;
3638 };
3639
3640 // linearize vector of vectors
3641 std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
3642 boxes_and_ranks;
3643
3644 for (unsigned rank = 0; rank < global_bboxes_to_be_used->size(); ++rank)
3645 for (const auto &box : (*global_bboxes_to_be_used)[rank])
3646 if (is_valid(box))
3647 boxes_and_ranks.emplace_back(box, rank);
3648
3649 // pack boxes into r-tree
3650 const auto tree = pack_rtree(boxes_and_ranks);
3651
3652 // loop over all entities
3653 for (unsigned int i = 0; i < entities.size(); ++i)
3654 {
3655 // create a bounding box with tolerance
3656 const auto bb =
3657 BoundingBox<spacedim>(entities[i]).create_extended(tolerance);
3658
3659 // determine ranks potentially owning point/bounding box
3660 std::set<unsigned int> my_ranks;
3661
3662 for (const auto &box_and_rank :
3663 tree | boost::geometry::index::adaptors::queried(
3664 boost::geometry::index::intersects(bb)))
3665 my_ranks.insert(box_and_rank.second);
3666
3667 for (const auto rank : my_ranks)
3668 ranks_and_indices.emplace_back(rank, i);
3669 }
3670 };
3671
3672 if constexpr (use_arborx)
3673 {
3674 if (global_bboxes.size() == 1)
3675 {
3676 ArborXWrappers::DistributedTree distributed_tree(
3677 comm, global_bboxes[0]);
3678 std::vector<BoundingBox<spacedim>> query_bounding_boxes;
3679 for (const auto &entity : entities)
3680 query_bounding_boxes.emplace_back(
3681 BoundingBox<spacedim>(entity).create_extended(tolerance));
3682
3684 query_bounding_boxes);
3685 const auto &[indices_ranks, offsets] =
3686 distributed_tree.query(bb_intersect);
3687
3688 for (unsigned long int i = 0; i < offsets.size() - 1; ++i)
3689 {
3690 std::set<unsigned int> my_ranks;
3691 for (int j = offsets[i]; j < offsets[i + 1]; ++j)
3692 my_ranks.insert(indices_ranks[j].second);
3693
3694 for (const auto rank : my_ranks)
3695 ranks_and_indices.emplace_back(rank, i);
3696 }
3697 }
3698 else
3699 {
3700 // global_bboxes.size()>1
3701 process_bboxes();
3702 }
3703 }
3704 else
3705 {
3706 // No ArborX
3707 process_bboxes();
3708 }
3709
3710
3711 // convert to CRS
3712 std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
3713
3714 std::vector<unsigned int> ranks;
3715 std::vector<unsigned int> ptr;
3716 std::vector<unsigned int> indices;
3717
3718 unsigned int current_rank = numbers::invalid_unsigned_int;
3719
3720 for (const std::pair<unsigned int, unsigned int> &i : ranks_and_indices)
3721 {
3722 if (current_rank != i.first)
3723 {
3724 current_rank = i.first;
3725 ranks.push_back(current_rank);
3726 ptr.push_back(indices.size());
3727 }
3728
3729 indices.push_back(i.second);
3730 }
3731 ptr.push_back(indices.size());
3732
3733 return {std::move(ranks), std::move(ptr), std::move(indices)};
3734 }
3735
3736
3737
3738 template <int dim, int spacedim>
3739 std::vector<
3740 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
3741 Point<dim>>>
3743 const Cache<dim, spacedim> &cache,
3744 const Point<spacedim> &point,
3746 const std::vector<bool> &marked_vertices,
3747 const double tolerance,
3748 const bool enforce_unique_mapping)
3749 {
3750 std::vector<
3751 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
3752 Point<dim>>>
3753 locally_owned_active_cells_around_point;
3754
3755 const auto first_cell = GridTools::find_active_cell_around_point(
3756 cache.get_mapping(),
3757 cache.get_triangulation(),
3758 point,
3759 cache.get_vertex_to_cell_map(),
3761 cell_hint,
3764 tolerance,
3766
3767 const unsigned int my_rank = Utilities::MPI::this_mpi_process(
3769
3770 cell_hint = first_cell.first;
3771 if (cell_hint.state() == IteratorState::valid)
3772 {
3773 const auto active_cells_around_point =
3774 GridTools::find_all_active_cells_around_point(
3775 cache.get_mapping(),
3776 cache.get_triangulation(),
3777 point,
3778 tolerance,
3779 first_cell,
3780 &cache.get_vertex_to_cell_map());
3781
3782 if (enforce_unique_mapping)
3783 {
3784 // check if the rank of this process is the lowest of all cells
3785 // if not, the other process will handle this cell and we don't
3786 // have to do here anything in the case of unique mapping
3787 unsigned int lowes_rank = numbers::invalid_unsigned_int;
3788
3789 for (const auto &cell : active_cells_around_point)
3790 lowes_rank = std::min(lowes_rank, cell.first->subdomain_id());
3791
3792 if (lowes_rank != my_rank)
3793 return {};
3794 }
3795
3796 locally_owned_active_cells_around_point.reserve(
3797 active_cells_around_point.size());
3798
3799 for (const auto &cell : active_cells_around_point)
3800 if (cell.first->is_locally_owned())
3801 locally_owned_active_cells_around_point.push_back(cell);
3802 }
3803
3804 std::sort(locally_owned_active_cells_around_point.begin(),
3805 locally_owned_active_cells_around_point.end(),
3806 [](const auto &a, const auto &b) { return a.first < b.first; });
3807
3808 if (enforce_unique_mapping &&
3809 locally_owned_active_cells_around_point.size() > 1)
3810 // in the case of unique mapping, we only need a single cell
3811 return {locally_owned_active_cells_around_point.front()};
3812 else
3813 return locally_owned_active_cells_around_point;
3814 }
3815
3816 template <int dim, int spacedim>
3821
3822 template <int dim, int spacedim>
3823 void
3825 {
3826 // before reshuffeling the data check if data.recv_components and
3827 // n_searched_points are in a valid state.
3828 Assert(n_searched_points != numbers::invalid_unsigned_int,
3830 Assert(recv_components.empty() ||
3831 std::get<1>(*std::max_element(recv_components.begin(),
3832 recv_components.end(),
3833 [](const auto &a, const auto &b) {
3834 return std::get<1>(a) <
3835 std::get<1>(b);
3836 })) < n_searched_points,
3838
3839 send_ranks.clear();
3840 recv_ranks.clear();
3841 send_ptrs.clear();
3842 recv_ptrs.clear();
3843
3844 if (true)
3845 {
3846 // sort according to rank (and point index and cell) -> make
3847 // deterministic
3848 std::sort(send_components.begin(),
3849 send_components.end(),
3850 [&](const auto &a, const auto &b) {
3851 if (std::get<1>(a) != std::get<1>(b)) // rank
3852 return std::get<1>(a) < std::get<1>(b);
3853
3854 if (std::get<2>(a) != std::get<2>(b)) // point index
3855 return std::get<2>(a) < std::get<2>(b);
3856
3857 return std::get<0>(a) < std::get<0>(b); // cell
3858 });
3859
3860 // perform enumeration and extract rank information
3861 for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
3862 i < send_components.size();
3863 ++i)
3864 {
3865 std::get<5>(send_components[i]) = i;
3866
3867 if (dummy != std::get<1>(send_components[i]))
3868 {
3869 dummy = std::get<1>(send_components[i]);
3870 send_ranks.push_back(dummy);
3871 send_ptrs.push_back(i);
3872 }
3873 }
3874 send_ptrs.push_back(send_components.size());
3875
3876 // sort according to cell, rank, point index (while keeping
3877 // partial ordering)
3878 std::sort(send_components.begin(),
3879 send_components.end(),
3880 [&](const auto &a, const auto &b) {
3881 if (std::get<0>(a) != std::get<0>(b))
3882 return std::get<0>(a) < std::get<0>(b); // cell
3883
3884 if (std::get<1>(a) != std::get<1>(b))
3885 return std::get<1>(a) < std::get<1>(b); // rank
3886
3887 if (std::get<2>(a) != std::get<2>(b))
3888 return std::get<2>(a) < std::get<2>(b); // point index
3889
3890 return std::get<5>(a) < std::get<5>(b); // enumeration
3891 });
3892 }
3893
3894 if (recv_components.size() > 0)
3895 {
3896 // sort according to rank (and point index) -> make deterministic
3897 std::sort(recv_components.begin(),
3898 recv_components.end(),
3899 [&](const auto &a, const auto &b) {
3900 if (std::get<0>(a) != std::get<0>(b))
3901 return std::get<0>(a) < std::get<0>(b); // rank
3902
3903 return std::get<1>(a) < std::get<1>(b); // point index
3904 });
3905
3906 // perform enumeration and extract rank information
3907 for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
3908 i < recv_components.size();
3909 ++i)
3910 {
3911 std::get<2>(recv_components[i]) = i;
3912
3913 if (dummy != std::get<0>(recv_components[i]))
3914 {
3915 dummy = std::get<0>(recv_components[i]);
3916 recv_ranks.push_back(dummy);
3917 recv_ptrs.push_back(i);
3918 }
3919 }
3920 recv_ptrs.push_back(recv_components.size());
3921
3922 // sort according to point index and rank (while keeping partial
3923 // ordering)
3924 std::sort(recv_components.begin(),
3925 recv_components.end(),
3926 [&](const auto &a, const auto &b) {
3927 if (std::get<1>(a) != std::get<1>(b))
3928 return std::get<1>(a) < std::get<1>(b); // point index
3929
3930 if (std::get<0>(a) != std::get<0>(b))
3931 return std::get<0>(a) < std::get<0>(b); // rank
3932
3933 return std::get<2>(a) < std::get<2>(b); // enumeration
3934 });
3935 }
3936 }
3937
3938
3939
3940 template <int dim, int spacedim>
3944 const std::vector<Point<spacedim>> &points,
3945 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
3946 const std::vector<bool> &marked_vertices,
3947 const double tolerance,
3948 const bool perform_handshake,
3949 const bool enforce_unique_mapping)
3950 {
3952 result.n_searched_points = points.size();
3953
3954 auto &send_components = result.send_components;
3955 auto &recv_components = result.recv_components;
3956
3957 const auto comm = cache.get_triangulation().get_communicator();
3958
3959 const auto potential_owners = internal::guess_owners_of_entities(
3960 comm, global_bboxes, points, tolerance);
3961
3962 const auto &potential_owners_ranks = std::get<0>(potential_owners);
3963 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
3964 const auto &potential_owners_indices = std::get<2>(potential_owners);
3965
3966 auto cell_hint = cache.get_triangulation().begin_active();
3967
3968 const auto translate = [&](const unsigned int other_rank) {
3969 const auto ptr = std::find(potential_owners_ranks.begin(),
3970 potential_owners_ranks.end(),
3971 other_rank);
3972
3973 Assert(ptr != potential_owners_ranks.end(), ExcInternalError());
3974
3975 const auto other_rank_index =
3976 std::distance(potential_owners_ranks.begin(), ptr);
3977
3978 return other_rank_index;
3979 };
3980
3981 Assert(
3982 (marked_vertices.empty()) ||
3983 (marked_vertices.size() == cache.get_triangulation().n_vertices()),
3984 ExcMessage(
3985 "The marked_vertices vector has to be either empty or its size has "
3986 "to equal the number of vertices of the triangulation."));
3987
3988 using RequestType = std::vector<std::pair<unsigned int, Point<spacedim>>>;
3989 using AnswerType = std::vector<unsigned int>;
3990
3991 // In the case that a marked_vertices vector has been given and none
3992 // of its entries is true, we know that this process does not own
3993 // any of the incoming points (and it will not send any data) so
3994 // that we can take a short cut.
3995 const bool has_relevant_vertices =
3996 (marked_vertices.empty()) ||
3997 (std::find(marked_vertices.begin(), marked_vertices.end(), true) !=
3998 marked_vertices.end());
3999
4000 const auto create_request = [&](const unsigned int other_rank) {
4001 const auto other_rank_index = translate(other_rank);
4002
4003 RequestType request;
4004 request.reserve(potential_owners_ptrs[other_rank_index + 1] -
4005 potential_owners_ptrs[other_rank_index]);
4006
4007 for (unsigned int i = potential_owners_ptrs[other_rank_index];
4008 i < potential_owners_ptrs[other_rank_index + 1];
4009 ++i)
4010 request.emplace_back(potential_owners_indices[i],
4011 points[potential_owners_indices[i]]);
4012
4013 return request;
4014 };
4015
4016 const auto answer_request =
4017 [&](const unsigned int &other_rank,
4018 const RequestType &request) -> AnswerType {
4019 AnswerType answer(request.size(), 0);
4020
4021 if (has_relevant_vertices)
4022 {
4023 cell_hint = cache.get_triangulation().begin_active();
4024
4025 for (unsigned int i = 0; i < request.size(); ++i)
4026 {
4027 const auto &index_and_point = request[i];
4028
4029 const auto cells_and_reference_positions =
4031 cache,
4032 index_and_point.second,
4033 cell_hint,
4035 tolerance,
4036 enforce_unique_mapping);
4037
4038 if (cell_hint.state() != IteratorState::valid)
4039 cell_hint = cache.get_triangulation().begin_active();
4040
4041 for (const auto &cell_and_reference_position :
4042 cells_and_reference_positions)
4043 {
4044 const auto cell = cell_and_reference_position.first;
4045 auto reference_position =
4046 cell_and_reference_position.second;
4047
4048 reference_position =
4049 cell->reference_cell().closest_point(reference_position);
4050
4051 send_components.emplace_back(
4052 std::pair<int, int>(cell->level(), cell->index()),
4053 other_rank,
4054 index_and_point.first,
4055 reference_position,
4056 index_and_point.second,
4058 }
4059
4060 answer[i] = cells_and_reference_positions.size();
4061 }
4062 }
4063
4064 if (perform_handshake)
4065 return answer;
4066 else
4067 return {};
4068 };
4069
4070 const auto process_answer = [&](const unsigned int other_rank,
4071 const AnswerType &answer) {
4072 if (perform_handshake)
4073 {
4074 const auto other_rank_index = translate(other_rank);
4075
4076 for (unsigned int i = 0; i < answer.size(); ++i)
4077 for (unsigned int j = 0; j < answer[i]; ++j)
4078 recv_components.emplace_back(
4079 other_rank,
4080 potential_owners_indices
4081 [i + potential_owners_ptrs[other_rank_index]],
4083 }
4084 };
4085
4086 Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
4087 potential_owners_ranks,
4088 create_request,
4089 answer_request,
4090 process_answer,
4091 comm);
4092
4093 result.finalize_setup();
4094
4095 return result;
4096 }
4097
4098
4099
4100 template <int structdim, int spacedim>
4101 template <int dim>
4102 DistributedComputePointLocationsInternal<dim, spacedim>
4105 const unsigned int n_points_1D,
4108 std::vector<Quadrature<spacedim>> *mapped_quadratures_recv_comp,
4109 const bool consistent_numbering_of_sender_and_receiver) const
4110 {
4111 using CellIterator =
4113
4114 if (mapped_quadratures_recv_comp != nullptr)
4115 {
4116 AssertDimension(mapped_quadratures_recv_comp->size(), 0);
4117 mapped_quadratures_recv_comp->reserve(recv_components.size());
4118 }
4119
4121 spacedim>
4122 result;
4123
4124 // We need quadrature rules for the intersections. We are using a
4125 // QGaussSimplex quadrature rule since CGAL always returns simplices
4126 // as intersections.
4127 const QGaussSimplex<structdim> quadrature(n_points_1D);
4128
4129 // Resulting quadrature points get different indices. In the case the
4130 // requested intersections are unique also the resulting quadrature
4131 // points are unique and we can simply number the points in an
4132 // ascending way.
4133 for (const auto &recv_component : recv_components)
4134 {
4135 // dependent on the size of the intersection an empty quadrature
4136 // is returned. Therefore, we have to compute the quadrature also
4137 // here.
4138 const Quadrature<spacedim> &quad =
4140 std::get<2>(recv_component));
4141
4142 for (unsigned int i = 0; i < quad.size(); ++i)
4143 {
4144 // the third component of result.recv_components is not needed
4145 // before finalize_setup.
4146 result.recv_components.emplace_back(
4147 std::get<0>(recv_component),
4148 result.recv_components.size(), // number of point
4150 }
4151
4152 // append quadrature
4153 if (mapped_quadratures_recv_comp != nullptr)
4154 mapped_quadratures_recv_comp->push_back(quad);
4155 }
4156
4157 // since empty quadratures might be present we have to set the number
4158 // of searched points after inserting the point indices into
4159 // recv_components
4160 result.n_searched_points = result.recv_components.size();
4161
4162 // send_ranks, counter, and indices_of_rank is only needed if
4163 // consistent_numbering_of_sender_and_receiver==true
4164 // indices_of_rank is always empty if deal.II is compiled without MPI
4165 std::map<unsigned int, std::vector<unsigned int>> indices_of_rank;
4166 std::map<unsigned int, unsigned int> counter;
4167 std::set<unsigned int> send_ranks;
4168 if (consistent_numbering_of_sender_and_receiver)
4169 {
4170 for (const auto &sc : send_components)
4171 send_ranks.insert(std::get<1>(sc));
4172
4173 for (const auto rank : send_ranks)
4174 counter[rank] = 0;
4175
4176 // indices assigned at recv side needed to fill send_components
4177 indices_of_rank = communicate_indices(result.recv_components,
4179 }
4180
4181 for (const auto &send_component : send_components)
4182 {
4183 const CellIterator cell(&tria,
4184 std::get<0>(send_component).first,
4185 std::get<0>(send_component).second);
4186
4187 const Quadrature<spacedim> &quad =
4189 std::get<3>(send_component));
4190
4191 const auto rank = std::get<1>(send_component);
4192
4193 for (unsigned int q = 0; q < quad.size(); ++q)
4194 {
4195 // the fifth component of result.send_components is filled
4196 // during sorting the data and initializing the CRS structures
4197 result.send_components.emplace_back(std::make_tuple(
4198 std::get<0>(send_component),
4199 rank,
4200 indices_of_rank.empty() ?
4201 result.send_components.size() :
4202 indices_of_rank.at(rank)[counter.at(rank)],
4203 mapping.transform_real_to_unit_cell(cell, quad.point(q)),
4204 quad.point(q),
4206
4207 if (!indices_of_rank.empty())
4208 ++counter[rank];
4209 }
4210 }
4211
4212 result.finalize_setup();
4213
4214 return result;
4215 }
4216
4217
4218
4219 template <int structdim, int spacedim>
4220 std::map<unsigned int, std::vector<unsigned int>>
4223 const std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
4224 &point_recv_components,
4225 const MPI_Comm comm) const
4226 {
4227#ifndef DEAL_II_WITH_MPI
4228 Assert(false, ExcNeedsMPI());
4229 (void)point_recv_components;
4230 (void)comm;
4231 return {};
4232#else
4233 // since we are converting to DistributedComputePointLocationsInternal
4234 // we use the RPE tag
4235 const auto mpi_tag =
4237
4238 const unsigned int my_rank = Utilities::MPI::this_mpi_process(comm);
4239
4240 std::set<unsigned int> send_ranks;
4241 for (const auto &sc : send_components)
4242 send_ranks.insert(std::get<1>(sc));
4243 std::set<unsigned int> recv_ranks;
4244 for (const auto &rc : recv_components)
4245 recv_ranks.insert(std::get<0>(rc));
4246
4247 std::vector<MPI_Request> requests;
4248 requests.reserve(send_ranks.size());
4249
4250 // rank to used indices on the rank needed on sending side
4251 std::map<unsigned int, std::vector<unsigned int>> indices_of_rank;
4252 indices_of_rank[my_rank] = std::vector<unsigned int>();
4253
4254 // rank to used indices on the rank known on recv side
4255 std::map<unsigned int, std::vector<unsigned int>> send_indices_of_rank;
4256 for (const auto rank : recv_ranks)
4257 if (rank != my_rank)
4258 send_indices_of_rank[rank] = std::vector<unsigned int>();
4259
4260 // fill the maps
4261 for (const auto &point_recv_component : point_recv_components)
4262 {
4263 const auto rank = std::get<0>(point_recv_component);
4264 const auto idx = std::get<1>(point_recv_component);
4265
4266 if (rank == my_rank)
4267 indices_of_rank[rank].emplace_back(idx);
4268 else
4269 send_indices_of_rank[rank].emplace_back(idx);
4270 }
4271
4272 // send indices to the ranks we normally receive from
4273 for (const auto rank : recv_ranks)
4274 {
4275 if (rank == my_rank)
4276 continue;
4277
4278 auto buffer = Utilities::pack(send_indices_of_rank[rank], false);
4279
4280 requests.push_back(MPI_Request());
4281
4282 const int ierr = MPI_Isend(buffer.data(),
4283 buffer.size(),
4284 MPI_CHAR,
4285 rank,
4286 mpi_tag,
4287 comm,
4288 &requests.back());
4289 AssertThrowMPI(ierr);
4290 }
4291
4292 // receive indices at the ranks we normally send from
4293 for (const auto rank : send_ranks)
4294 {
4295 if (rank == my_rank)
4296 continue;
4297
4298 MPI_Status status;
4299 int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, comm, &status);
4300 AssertThrowMPI(ierr);
4301
4302 int message_length;
4303 ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
4304 AssertThrowMPI(ierr);
4305
4306 std::vector<char> buffer(message_length);
4307
4308 ierr = MPI_Recv(buffer.data(),
4309 buffer.size(),
4310 MPI_CHAR,
4311 status.MPI_SOURCE,
4312 mpi_tag,
4313 comm,
4314 MPI_STATUS_IGNORE);
4315 AssertThrowMPI(ierr);
4316
4317 indices_of_rank[status.MPI_SOURCE] =
4318 Utilities::unpack<std::vector<unsigned int>>(buffer, false);
4319 }
4320
4321 // make sure all messages have been sent
4322 const int ierr =
4323 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4324 AssertThrowMPI(ierr);
4325
4326 return indices_of_rank;
4327#endif
4328 }
4329
4330
4331
4332 template <int structdim, int dim, int spacedim>
4335 const Cache<dim, spacedim> &cache,
4336 const std::vector<std::vector<Point<spacedim>>> &intersection_requests,
4337 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
4338 const std::vector<bool> &marked_vertices,
4339 const double tolerance)
4340 {
4341 using IntersectionRequest = std::vector<Point<spacedim>>;
4342
4343 using IntersectionAnswer =
4345 structdim,
4346 spacedim>::IntersectionType;
4347
4348 const auto comm = cache.get_triangulation().get_communicator();
4349
4351 result;
4352
4353 auto &send_components = result.send_components;
4354 auto &recv_components = result.recv_components;
4355 auto &recv_ptrs = result.recv_ptrs;
4356
4357 // search for potential owners
4358 const auto potential_owners = internal::guess_owners_of_entities(
4359 comm, global_bboxes, intersection_requests, tolerance);
4360
4361 const auto &potential_owners_ranks = std::get<0>(potential_owners);
4362 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
4363 const auto &potential_owners_indices = std::get<2>(potential_owners);
4364
4365 const auto translate = [&](const unsigned int other_rank) {
4366 const auto ptr = std::find(potential_owners_ranks.begin(),
4367 potential_owners_ranks.end(),
4368 other_rank);
4369
4370 Assert(ptr != potential_owners_ranks.end(), ExcInternalError());
4371
4372 const auto other_rank_index =
4373 std::distance(potential_owners_ranks.begin(), ptr);
4374
4375 return other_rank_index;
4376 };
4377
4378 Assert(
4379 (marked_vertices.empty()) ||
4380 (marked_vertices.size() == cache.get_triangulation().n_vertices()),
4381 ExcMessage(
4382 "The marked_vertices vector has to be either empty or its size has "
4383 "to equal the number of vertices of the triangulation."));
4384
4385 // In the case that a marked_vertices vector has been given and none
4386 // of its entries is true, we know that this process does not own
4387 // any of the incoming points (and it will not send any data) so
4388 // that we can take a short cut.
4389 const bool has_relevant_vertices =
4390 (marked_vertices.empty()) ||
4391 (std::find(marked_vertices.begin(), marked_vertices.end(), true) !=
4392 marked_vertices.end());
4393
4394 // intersection between two cells:
4395 // One rank requests all intersections of owning cell:
4396 // owning cell index, cgal vertices of cell
4397 using RequestType =
4398 std::vector<std::pair<unsigned int, IntersectionRequest>>;
4399 // Other ranks send back all found intersections for requesting cell:
4400 // requesting cell index, cgal vertices of found intersections
4401 using AnswerType =
4402 std::vector<std::pair<unsigned int, IntersectionAnswer>>;
4403
4404 const auto create_request = [&](const unsigned int other_rank) {
4405 const auto other_rank_index = translate(other_rank);
4406
4407 RequestType request;
4408 request.reserve(potential_owners_ptrs[other_rank_index + 1] -
4409 potential_owners_ptrs[other_rank_index]);
4410
4411 for (unsigned int i = potential_owners_ptrs[other_rank_index];
4412 i < potential_owners_ptrs[other_rank_index + 1];
4413 ++i)
4414 request.emplace_back(
4415 potential_owners_indices[i],
4416 intersection_requests[potential_owners_indices[i]]);
4417
4418 return request;
4419 };
4420
4421
4422 // TODO: this is potentially useful in many cases and it would be nice to
4423 // have cache.get_locally_owned_cell_bounding_boxes_rtree(marked_vertices)
4424 const auto construct_locally_owned_cell_bounding_boxes_rtree =
4425 [&cache](const std::vector<bool> &marked_verts) {
4426 const auto cell_marked = [&marked_verts](const auto &cell) {
4427 for (const unsigned int v : cell->vertex_indices())
4428 if (marked_verts[cell->vertex_index(v)])
4429 return true;
4430 return false;
4431 };
4432
4433 const auto &boxes_and_cells =
4435
4436 if (marked_verts.empty())
4437 return boxes_and_cells;
4438
4439 std::vector<std::pair<
4442 potential_boxes_and_cells;
4443
4444 for (const auto &box_and_cell : boxes_and_cells)
4445 if (cell_marked(box_and_cell.second))
4446 potential_boxes_and_cells.emplace_back(box_and_cell);
4447
4448 return pack_rtree(potential_boxes_and_cells);
4449 };
4450
4451
4452 RTree<
4453 std::pair<BoundingBox<spacedim>,
4455 marked_cell_tree;
4456
4457 const auto answer_request =
4458 [&](const unsigned int &other_rank,
4459 const RequestType &request) -> AnswerType {
4460 AnswerType answer;
4461
4462 if (has_relevant_vertices)
4463 {
4464 if (marked_cell_tree.empty())
4465 {
4466 marked_cell_tree =
4467 construct_locally_owned_cell_bounding_boxes_rtree(
4469 }
4470
4471 // process requests
4472 for (unsigned int i = 0; i < request.size(); ++i)
4473 {
4474 // create a bounding box with tolerance
4475 const auto bb = BoundingBox<spacedim>(request[i].second)
4476 .create_extended(tolerance);
4477
4478 for (const auto &box_cell :
4479 marked_cell_tree |
4480 boost::geometry::index::adaptors::queried(
4481 boost::geometry::index::intersects(bb)))
4482 {
4483#ifdef DEAL_II_WITH_CGAL
4484 const auto &cell = box_cell.second;
4485 const auto &request_index = request[i].first;
4486 auto requested_intersection = request[i].second;
4487 CGALWrappers::resort_dealii_vertices_to_cgal_order(
4488 structdim, requested_intersection);
4489
4490 const auto &try_intersection =
4491 CGALWrappers::get_vertices_in_cgal_order(
4492 cell, cache.get_mapping());
4493
4494 const auto &found_intersections = CGALWrappers::
4495 compute_intersection_of_cells<dim, structdim, spacedim>(
4496 try_intersection, requested_intersection, tolerance);
4497
4498 if (found_intersections.size() > 0)
4499 {
4500 for (const auto &found_intersection :
4501 found_intersections)
4502 {
4503 answer.emplace_back(request_index,
4504 found_intersection);
4505
4506 send_components.emplace_back(
4507 std::make_pair(cell->level(), cell->index()),
4508 other_rank,
4509 request_index,
4510 found_intersection);
4511 }
4512 }
4513#else
4514 (void)other_rank;
4515 (void)box_cell;
4516 Assert(false, ExcNeedsCGAL());
4517#endif
4518 }
4519 }
4520 }
4521
4522 return answer;
4523 };
4524
4525 const auto process_answer = [&](const unsigned int other_rank,
4526 const AnswerType &answer) {
4527 for (unsigned int i = 0; i < answer.size(); ++i)
4528 recv_components.emplace_back(other_rank,
4529 answer[i].first,
4530 answer[i].second);
4531 };
4532
4533 Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
4534 potential_owners_ranks,
4535 create_request,
4536 answer_request,
4537 process_answer,
4538 comm);
4539
4540 // sort according to 1) intersection index and 2) rank (keeping the order
4541 // of recv components with same indices and ranks)
4542 std::stable_sort(recv_components.begin(),
4543 recv_components.end(),
4544 [&](const auto &a, const auto &b) {
4545 // intersection index
4546 if (std::get<1>(a) != std::get<1>(b))
4547 return std::get<1>(a) < std::get<1>(b);
4548
4549 // rank
4550 return std::get<0>(a) < std::get<0>(b);
4551 });
4552
4553 // sort according to 1) rank and 2) intersection index (keeping the
4554 // order of recv components with same indices and ranks)
4555 std::stable_sort(send_components.begin(),
4556 send_components.end(),
4557 [&](const auto &a, const auto &b) {
4558 // rank
4559 if (std::get<1>(a) != std::get<1>(b))
4560 return std::get<1>(a) < std::get<1>(b);
4561
4562 // intersection idx
4563 return std::get<2>(a) < std::get<2>(b);
4564 });
4565
4566 // construct recv_ptrs
4567 recv_ptrs.assign(intersection_requests.size() + 1, 0);
4568 for (const auto &rc : recv_components)
4569 ++recv_ptrs[std::get<1>(rc) + 1];
4570 for (unsigned int i = 0; i < intersection_requests.size(); ++i)
4571 recv_ptrs[i + 1] += recv_ptrs[i];
4572
4573 return result;
4574 }
4575
4576 } // namespace internal
4577
4578
4579
4580 template <int spacedim>
4581 unsigned int
4582 find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
4583 const Point<spacedim> &p)
4584 {
4585 auto id_and_v = std::min_element(
4586 vertices.begin(),
4587 vertices.end(),
4588 [&](const std::pair<const unsigned int, Point<spacedim>> &p1,
4589 const std::pair<const unsigned int, Point<spacedim>> &p2) -> bool {
4590 return p1.second.distance(p) < p2.second.distance(p);
4591 });
4592 return id_and_v->first;
4593 }
4594
4595
4596 template <int dim, int spacedim>
4597 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4598 Point<dim>>
4600 const Cache<dim, spacedim> &cache,
4601 const Point<spacedim> &p,
4603 &cell_hint,
4604 const std::vector<bool> &marked_vertices,
4605 const double tolerance)
4606 {
4607 const auto &mesh = cache.get_triangulation();
4608 const auto &mapping = cache.get_mapping();
4609 const auto &vertex_to_cells = cache.get_vertex_to_cell_map();
4610 const auto &vertex_to_cell_centers =
4612 const auto &used_vertices_rtree = cache.get_used_vertices_rtree();
4613
4615 mesh,
4616 p,
4617 vertex_to_cells,
4618 vertex_to_cell_centers,
4619 cell_hint,
4621 used_vertices_rtree,
4622 tolerance);
4623 }
4624
4625 template <int spacedim>
4626 std::vector<std::vector<BoundingBox<spacedim>>>
4628 const std::vector<BoundingBox<spacedim>> &local_bboxes,
4629 const MPI_Comm mpi_communicator)
4630 {
4631#ifndef DEAL_II_WITH_MPI
4632 (void)local_bboxes;
4633 (void)mpi_communicator;
4634 Assert(false,
4635 ExcMessage(
4636 "GridTools::exchange_local_bounding_boxes() requires MPI."));
4637 return {};
4638#else
4639 // Step 1: preparing data to be sent
4640 unsigned int n_bboxes = local_bboxes.size();
4641 // Dimension of the array to be exchanged (number of double)
4642 int n_local_data = 2 * spacedim * n_bboxes;
4643 // data array stores each entry of each point describing the bounding
4644 // boxes
4645 std::vector<double> loc_data_array(n_local_data);
4646 for (unsigned int i = 0; i < n_bboxes; ++i)
4647 for (unsigned int d = 0; d < spacedim; ++d)
4648 {
4649 // Extracting the coordinates of each boundary point
4650 loc_data_array[2 * i * spacedim + d] =
4651 local_bboxes[i].get_boundary_points().first[d];
4652 loc_data_array[2 * i * spacedim + spacedim + d] =
4653 local_bboxes[i].get_boundary_points().second[d];
4654 }
4655
4656 // Step 2: exchanging the size of local data
4657 unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_communicator);
4658
4659 // Vector to store the size of loc_data_array for every process
4660 std::vector<int> size_all_data(n_procs);
4661
4662 // Exchanging the number of bboxes
4663 int ierr = MPI_Allgather(&n_local_data,
4664 1,
4665 MPI_INT,
4666 size_all_data.data(),
4667 1,
4668 MPI_INT,
4669 mpi_communicator);
4670 AssertThrowMPI(ierr);
4671
4672 // Now computing the displacement, relative to recvbuf,
4673 // at which to store the incoming data
4674 std::vector<int> rdispls(n_procs);
4675 rdispls[0] = 0;
4676 for (unsigned int i = 1; i < n_procs; ++i)
4677 rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
4678
4679 // Step 3: exchange the data and bounding boxes:
4680 // Allocating a vector to contain all the received data
4681 std::vector<double> data_array(rdispls.back() + size_all_data.back());
4682
4683 ierr = MPI_Allgatherv(loc_data_array.data(),
4684 n_local_data,
4685 MPI_DOUBLE,
4686 data_array.data(),
4687 size_all_data.data(),
4688 rdispls.data(),
4689 MPI_DOUBLE,
4690 mpi_communicator);
4691 AssertThrowMPI(ierr);
4692
4693 // Step 4: create the array of bboxes for output
4694 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
4695 unsigned int begin_idx = 0;
4696 for (unsigned int i = 0; i < n_procs; ++i)
4697 {
4698 // Number of local bounding boxes
4699 unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
4700 global_bboxes[i].resize(n_bbox_i);
4701 for (unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
4702 {
4703 Point<spacedim> p1, p2; // boundary points for bbox
4704 for (unsigned int d = 0; d < spacedim; ++d)
4705 {
4706 p1[d] = data_array[begin_idx + 2 * bbox * spacedim + d];
4707 p2[d] =
4708 data_array[begin_idx + 2 * bbox * spacedim + spacedim + d];
4709 }
4710 BoundingBox<spacedim> loc_bbox(std::make_pair(p1, p2));
4711 global_bboxes[i][bbox] = loc_bbox;
4712 }
4713 // Shifting the first index to the start of the next vector
4714 begin_idx += size_all_data[i];
4715 }
4716 return global_bboxes;
4717#endif // DEAL_II_WITH_MPI
4718 }
4719
4720
4721
4722 template <int spacedim>
4725 const std::vector<BoundingBox<spacedim>> &local_description,
4726 const MPI_Comm mpi_communicator)
4727 {
4728#ifndef DEAL_II_WITH_MPI
4729 (void)mpi_communicator;
4730 // Building a tree with the only boxes available without MPI
4731 std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> boxes_index(
4732 local_description.size());
4733 // Adding to each box the rank of the process owning it
4734 for (unsigned int i = 0; i < local_description.size(); ++i)
4735 boxes_index[i] = std::make_pair(local_description[i], 0u);
4736 return pack_rtree(boxes_index);
4737#else
4738 // Exchanging local bounding boxes
4739 const std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes =
4740 Utilities::MPI::all_gather(mpi_communicator, local_description);
4741
4742 // Preparing to flatten the vector
4743 const unsigned int n_procs =
4744 Utilities::MPI::n_mpi_processes(mpi_communicator);
4745 // The i'th element of the following vector contains the index of the
4746 // first local bounding box from the process of rank i
4747 std::vector<unsigned int> bboxes_position(n_procs);
4748
4749 unsigned int tot_bboxes = 0;
4750 for (const auto &process_bboxes : global_bboxes)
4751 tot_bboxes += process_bboxes.size();
4752
4753 // Now flattening the vector
4754 std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
4755 flat_global_bboxes;
4756 flat_global_bboxes.reserve(tot_bboxes);
4757 unsigned int process_index = 0;
4758 for (const auto &process_bboxes : global_bboxes)
4759 {
4760 // Initialize a vector containing bounding boxes and rank of a process
4761 std::vector<std::pair<BoundingBox<spacedim>, unsigned int>>
4762 boxes_and_indices(process_bboxes.size());
4763
4764 // Adding to each box the rank of the process owning it
4765 for (unsigned int i = 0; i < process_bboxes.size(); ++i)
4766 boxes_and_indices[i] =
4767 std::make_pair(process_bboxes[i], process_index);
4768
4769 flat_global_bboxes.insert(flat_global_bboxes.end(),
4770 boxes_and_indices.begin(),
4771 boxes_and_indices.end());
4772
4773 ++process_index;
4774 }
4775
4776 // Build a tree out of the bounding boxes. We avoid using the
4777 // insert method so that boost uses the packing algorithm
4778 return RTree<std::pair<BoundingBox<spacedim>, unsigned int>>(
4779 flat_global_bboxes.begin(), flat_global_bboxes.end());
4780#endif // DEAL_II_WITH_MPI
4781 }
4782
4783
4784
4785 template <int dim, int spacedim>
4786 void
4789 std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
4790 std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group)
4791 {
4792 // 1) determine for each vertex a vertex it coincides with and
4793 // put it into a map
4794 {
4795 // loop over all periodic face pairs
4796 for (const auto &pair : tria.get_periodic_face_map())
4797 {
4798 if (pair.first.first->level() != pair.second.first.first->level())
4799 continue;
4800
4801 const auto face_a = pair.first.first->face(pair.first.second);
4802 const auto face_b =
4803 pair.second.first.first->face(pair.second.first.second);
4804 const auto reference_cell = pair.first.first->reference_cell();
4805 const auto face_reference_cell = face_a->reference_cell();
4806 const unsigned char combined_orientation = pair.second.second;
4807 const unsigned char inverse_combined_orientation =
4808 face_reference_cell.get_inverse_combined_orientation(
4809 combined_orientation);
4810
4811 AssertDimension(face_a->n_vertices(), face_b->n_vertices());
4812
4813 // loop over all vertices on face
4814 for (unsigned int i = 0; i < face_a->n_vertices(); ++i)
4815 {
4816 // find the right local vertex index for the second face
4817 const unsigned int j =
4818 reference_cell.standard_to_real_face_vertex(
4819 i, pair.first.second, inverse_combined_orientation);
4820
4821 // get vertex indices and store in map
4822 const auto vertex_a = face_a->vertex_index(i);
4823 const auto vertex_b = face_b->vertex_index(j);
4824 unsigned int temp = std::min(vertex_a, vertex_b);
4825
4826 auto it_a = vertex_to_coinciding_vertex_group.find(vertex_a);
4827 if (it_a != vertex_to_coinciding_vertex_group.end())
4828 temp = std::min(temp, it_a->second);
4829
4830 auto it_b = vertex_to_coinciding_vertex_group.find(vertex_b);
4831 if (it_b != vertex_to_coinciding_vertex_group.end())
4832 temp = std::min(temp, it_b->second);
4833
4834 if (it_a != vertex_to_coinciding_vertex_group.end())
4835 it_a->second = temp;
4836 else
4837 vertex_to_coinciding_vertex_group[vertex_a] = temp;
4838
4839 if (it_b != vertex_to_coinciding_vertex_group.end())
4840 it_b->second = temp;
4841 else
4842 vertex_to_coinciding_vertex_group[vertex_b] = temp;
4843 }
4844 }
4845
4846 // 2) compress map: let vertices point to the coinciding vertex with
4847 // the smallest index
4848 for (auto &p : vertex_to_coinciding_vertex_group)
4849 {
4850 if (p.first == p.second)
4851 continue;
4852 unsigned int temp = p.second;
4853 while (temp != vertex_to_coinciding_vertex_group[temp])
4854 temp = vertex_to_coinciding_vertex_group[temp];
4855 p.second = temp;
4856 }
4857
4858 // 3) create a map: smallest index of coinciding index -> all
4859 // coinciding indices
4860 for (auto p : vertex_to_coinciding_vertex_group)
4861 coinciding_vertex_groups[p.second] = {};
4862
4863 for (auto p : vertex_to_coinciding_vertex_group)
4864 coinciding_vertex_groups[p.second].push_back(p.first);
4865 }
4866 }
4867
4868
4869
4870 template <int dim, int spacedim>
4871 std::map<unsigned int, std::set<::types::subdomain_id>>
4873 const Triangulation<dim, spacedim> &tria)
4874 {
4875 if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4876 &tria) == nullptr) // nothing to do for a serial triangulation
4877 return {};
4878
4879 // 1) collect for each vertex on periodic faces all vertices it coincides
4880 // with
4881 std::map<unsigned int, std::vector<unsigned int>> coinciding_vertex_groups;
4882 std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
4883
4885 coinciding_vertex_groups,
4886 vertex_to_coinciding_vertex_group);
4887
4888 // 2) collect vertices belonging to local cells
4889 std::vector<bool> vertex_of_own_cell(tria.n_vertices(), false);
4890 for (const auto &cell :
4892 for (const unsigned int v : cell->vertex_indices())
4893 vertex_of_own_cell[cell->vertex_index(v)] = true;
4894
4895 // 3) for each vertex belonging to a locally owned cell, find all ghost
4896 // neighbors (including the periodic own)
4897 std::map<unsigned int, std::set<types::subdomain_id>> result;
4898
4899 // loop over all active ghost cells
4900 for (const auto &cell : tria.active_cell_iterators())
4901 if (cell->is_ghost())
4902 {
4903 const types::subdomain_id owner = cell->subdomain_id();
4904
4905 // loop over all its vertices
4906 for (const unsigned int v : cell->vertex_indices())
4907 {
4908 // set owner if vertex belongs to a local cell
4909 if (vertex_of_own_cell[cell->vertex_index(v)])
4910 result[cell->vertex_index(v)].insert(owner);
4911
4912 // mark also nodes coinciding due to periodicity
4913 auto coinciding_vertex_group =
4914 vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
4915 if (coinciding_vertex_group !=
4916 vertex_to_coinciding_vertex_group.end())
4917 for (auto coinciding_vertex :
4918 coinciding_vertex_groups[coinciding_vertex_group->second])
4919 if (vertex_of_own_cell[coinciding_vertex])
4920 result[coinciding_vertex].insert(owner);
4921 }
4922 }
4923
4924 return result;
4925 }
4926
4927
4928
4929 namespace internal
4930 {
4931 template <int dim,
4932 unsigned int n_vertices,
4933 unsigned int n_sub_vertices,
4934 unsigned int n_configurations,
4935 unsigned int n_lines,
4936 unsigned int n_cols,
4937 typename value_type>
4938 void
4940 const std::array<unsigned int, n_configurations> &cut_line_table,
4942 const ndarray<unsigned int, n_lines, 2> &line_to_vertex_table,
4943 const std::vector<value_type> &ls_values,
4944 const std::vector<Point<dim>> &points,
4945 const std::vector<unsigned int> &mask,
4946 const double iso_level,
4947 const double tolerance,
4948 std::vector<Point<dim>> &vertices,
4949 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
4950 const bool write_back_cell_data)
4951 {
4952 // inspired by https://graphics.stanford.edu/~mdfisher/MarchingCubes.html
4953
4954 constexpr unsigned int X = static_cast<unsigned int>(-1);
4955
4956 // determine configuration
4957 unsigned int configuration = 0;
4958 for (unsigned int v = 0; v < n_vertices; ++v)
4959 if (ls_values[mask[v]] < iso_level)
4960 configuration |= (1 << v);
4961
4962 // cell is not cut (nothing to do)
4963 if (cut_line_table[configuration] == 0)
4964 return;
4965
4966 // helper function to determine where an edge (between index i and j) is
4967 // cut - see also: http://paulbourke.net/geometry/polygonise/
4968 const auto interpolate = [&](const unsigned int i, const unsigned int j) {
4969 if (std::abs(iso_level - ls_values[mask[i]]) < tolerance)
4970 return points[mask[i]];
4971 if (std::abs(iso_level - ls_values[mask[j]]) < tolerance)
4972 return points[mask[j]];
4973 if (std::abs(ls_values[mask[i]] - ls_values[mask[j]]) < tolerance)
4974 return points[mask[i]];
4975
4976 const double mu = (iso_level - ls_values[mask[i]]) /
4977 (ls_values[mask[j]] - ls_values[mask[i]]);
4978
4979 return Point<dim>(points[mask[i]] +
4980 mu * (points[mask[j]] - points[mask[i]]));
4981 };
4982
4983 // determine the position where edges are cut (if they are cut)
4984 std::array<Point<dim>, n_lines> vertex_list_all;
4985 for (unsigned int l = 0; l < n_lines; ++l)
4986 if (cut_line_table[configuration] & (1 << l))
4987 vertex_list_all[l] =
4988 interpolate(line_to_vertex_table[l][0], line_to_vertex_table[l][1]);
4989
4990 // merge duplicate vertices if possible
4991 unsigned int local_vertex_count = 0;
4992 std::array<Point<dim>, n_lines> vertex_list_reduced;
4993 std::array<unsigned int, n_lines> local_remap;
4994 std::fill(local_remap.begin(), local_remap.end(), X);
4995 for (int i = 0; new_line_table[configuration][i] != X; ++i)
4996 if (local_remap[new_line_table[configuration][i]] == X)
4997 {
4998 vertex_list_reduced[local_vertex_count] =
4999 vertex_list_all[new_line_table[configuration][i]];
5000 local_remap[new_line_table[configuration][i]] = local_vertex_count;
5001 ++local_vertex_count;
5002 }
5003
5004 // write back vertices
5005 const unsigned int n_vertices_old = vertices.size();
5006 for (unsigned int i = 0; i < local_vertex_count; ++i)
5007 vertices.push_back(vertex_list_reduced[i]);
5008
5009 // write back cells
5010 if (write_back_cell_data && dim > 1)
5011 {
5012 for (unsigned int i = 0; new_line_table[configuration][i] != X;
5013 i += n_sub_vertices)
5014 {
5015 cells.resize(cells.size() + 1);
5016 cells.back().vertices.resize(n_sub_vertices);
5017
5018 for (unsigned int v = 0; v < n_sub_vertices; ++v)
5019 cells.back().vertices[v] =
5020 local_remap[new_line_table[configuration][i + v]] +
5021 n_vertices_old;
5022 }
5023 }
5024 }
5025 } // namespace internal
5026
5027
5028
5029 template <int dim, typename VectorType>
5031 const Mapping<dim, dim> &mapping,
5032 const FiniteElement<dim, dim> &fe,
5033 const unsigned int n_subdivisions,
5034 const double tolerance)
5035 : n_subdivisions(n_subdivisions)
5036 , tolerance(tolerance)
5037 , fe_values(mapping,
5038 fe,
5039 create_quadrature_rule(n_subdivisions),
5041 {}
5042
5043
5044
5045 template <int dim, typename VectorType>
5048 const unsigned int n_subdivisions)
5049 {
5050 std::vector<Point<dim>> quadrature_points;
5051
5052 if (dim == 1)
5053 {
5054 for (unsigned int i = 0; i <= n_subdivisions; ++i)
5055 quadrature_points.emplace_back(1.0 / n_subdivisions * i);
5056 }
5057 else if (dim == 2)
5058 {
5059 for (unsigned int j = 0; j <= n_subdivisions; ++j)
5060 for (unsigned int i = 0; i <= n_subdivisions; ++i)
5061 quadrature_points.emplace_back(1.0 / n_subdivisions * i,
5062 1.0 / n_subdivisions * j);
5063 }
5064 else
5065 {
5066 for (unsigned int k = 0; k <= n_subdivisions; ++k)
5067 for (unsigned int j = 0; j <= n_subdivisions; ++j)
5068 for (unsigned int i = 0; i <= n_subdivisions; ++i)
5069 quadrature_points.emplace_back(1.0 / n_subdivisions * i,
5070 1.0 / n_subdivisions * j,
5071 1.0 / n_subdivisions * k);
5072 }
5073
5074
5075 return {quadrature_points};
5076 }
5077
5078
5079
5080 template <int dim, typename VectorType>
5081 void
5083 const DoFHandler<dim> &background_dof_handler,
5084 const VectorType &ls_vector,
5085 const double iso_level,
5086 std::vector<Point<dim>> &vertices,
5087 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const
5088 {
5090 dim > 1,
5091 ExcMessage(
5092 "Not implemented for dim==1. Use the alternative process()-function "
5093 "not returning a vector of CellData objects."));
5094
5095 for (const auto &cell : background_dof_handler.active_cell_iterators() |
5097 process_cell(cell, ls_vector, iso_level, vertices, cells);
5098 }
5099
5100 template <int dim, typename VectorType>
5101 void
5103 const DoFHandler<dim> &background_dof_handler,
5104 const VectorType &ls_vector,
5105 const double iso_level,
5106 std::vector<Point<dim>> &vertices) const
5107 {
5108 for (const auto &cell : background_dof_handler.active_cell_iterators() |
5110 process_cell(cell, ls_vector, iso_level, vertices);
5111
5112 delete_duplicated_vertices(vertices, 1e-10 /*tol*/);
5113 }
5114
5115
5116 template <int dim, typename VectorType>
5117 void
5119 const typename DoFHandler<dim>::active_cell_iterator &cell,
5120 const VectorType &ls_vector,
5121 const double iso_level,
5122 std::vector<Point<dim>> &vertices,
5123 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const
5124 {
5126 dim > 1,
5127 ExcMessage(
5128 "Not implemented for dim==1. Use the alternative process_cell()-function "
5129 "not returning a vector of CellData objects."));
5130
5131 std::vector<value_type> ls_values;
5132
5133 fe_values.reinit(cell);
5134 ls_values.resize(fe_values.n_quadrature_points);
5135 fe_values.get_function_values(ls_vector, ls_values);
5136 process_cell(
5137 ls_values, fe_values.get_quadrature_points(), iso_level, vertices, cells);
5138 }
5139
5140 template <int dim, typename VectorType>
5141 void
5143 const typename DoFHandler<dim>::active_cell_iterator &cell,
5144 const VectorType &ls_vector,
5145 const double iso_level,
5146 std::vector<Point<dim>> &vertices) const
5147 {
5148 // This vector is just a placeholder to reuse the process_cell function.
5149 std::vector<CellData<dim == 1 ? 1 : dim - 1>> dummy_cells;
5150
5151 std::vector<value_type> ls_values;
5152
5153 fe_values.reinit(cell);
5154 ls_values.resize(fe_values.n_quadrature_points);
5155 fe_values.get_function_values(ls_vector, ls_values);
5156
5157 process_cell(ls_values,
5158 fe_values.get_quadrature_points(),
5159 iso_level,
5160 vertices,
5161 dummy_cells,
5162 false /*don't write back cell data*/);
5163 }
5164
5165
5166 template <int dim, typename VectorType>
5167 void
5169 std::vector<value_type> &ls_values,
5170 const std::vector<Point<dim>> &points,
5171 const double iso_level,
5172 std::vector<Point<dim>> &vertices,
5173 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
5174 const bool write_back_cell_data) const
5175 {
5176 const unsigned p = n_subdivisions + 1;
5177
5178 if (dim == 1)
5179 {
5180 for (unsigned int i = 0; i < n_subdivisions; ++i)
5181 {
5182 std::vector<unsigned int> mask{i + 0, i + 1};
5183
5184 // check if a corner node is cut
5185 if (std::abs(iso_level - ls_values[mask[0]]) < tolerance)
5186 vertices.emplace_back(points[mask[0]]);
5187 else if (std::abs(iso_level - ls_values[mask[1]]) < tolerance)
5188 {
5189 if (i + 1 == n_subdivisions)
5190 vertices.emplace_back(points[mask[1]]);
5191 }
5192 // check if the edge is cut
5193 else if (((ls_values[mask[0]] > iso_level) &&
5194 (ls_values[mask[1]] < iso_level)) ||
5195 ((ls_values[mask[0]] < iso_level) &&
5196 (ls_values[mask[1]] > iso_level)))
5197 {
5198 // determine the interpolation weight (0<mu<1)
5199 const double mu = (iso_level - ls_values[mask[0]]) /
5200 (ls_values[mask[1]] - ls_values[mask[0]]);
5201
5202 // interpolate
5203 vertices.emplace_back(points[mask[0]] +
5204 mu * (points[mask[1]] - points[mask[0]]));
5205 }
5206 }
5207 }
5208 else if (dim == 2)
5209 {
5210 for (unsigned int j = 0; j < n_subdivisions; ++j)
5211 for (unsigned int i = 0; i < n_subdivisions; ++i)
5212 {
5213 std::vector<unsigned int> mask{p * (j + 0) + (i + 0),
5214 p * (j + 0) + (i + 1),
5215 p * (j + 1) + (i + 1),
5216 p * (j + 1) + (i + 0)};
5217
5218 process_sub_cell(ls_values,
5219 points,
5220 mask,
5221 iso_level,
5222 vertices,
5223 cells,
5224 write_back_cell_data);
5225 }
5226 }
5227 else if (dim == 3)
5228 {
5229 for (unsigned int k = 0; k < n_subdivisions; ++k)
5230 for (unsigned int j = 0; j < n_subdivisions; ++j)
5231 for (unsigned int i = 0; i < n_subdivisions; ++i)
5232 {
5233 std::vector<unsigned int> mask{
5234 p * p * (k + 0) + p * (j + 0) + (i + 0),
5235 p * p * (k + 0) + p * (j + 0) + (i + 1),
5236 p * p * (k + 0) + p * (j + 1) + (i + 1),
5237 p * p * (k + 0) + p * (j + 1) + (i + 0),
5238 p * p * (k + 1) + p * (j + 0) + (i + 0),
5239 p * p * (k + 1) + p * (j + 0) + (i + 1),
5240 p * p * (k + 1) + p * (j + 1) + (i + 1),
5241 p * p * (k + 1) + p * (j + 1) + (i + 0)};
5242
5243 process_sub_cell(ls_values,
5244 points,
5245 mask,
5246 iso_level,
5247 vertices,
5248 cells,
5249 write_back_cell_data);
5250 }
5251 }
5252 }
5253
5254
5255
5256 template <int dim, typename VectorType>
5257 void
5259 const std::vector<value_type> &ls_values,
5260 const std::vector<Point<2>> &points,
5261 const std::vector<unsigned int> &mask,
5262 const double iso_level,
5263 std::vector<Point<2>> &vertices,
5264 std::vector<CellData<1>> &cells,
5265 const bool write_back_cell_data) const
5266 {
5267 // set up dimension-dependent sizes and tables
5268 constexpr unsigned int n_vertices = 4;
5269 constexpr unsigned int n_sub_vertices = 2;
5270 constexpr unsigned int n_lines = 4;
5271 constexpr unsigned int n_configurations = Utilities::pow(2, n_vertices);
5272 constexpr unsigned int X = static_cast<unsigned int>(-1);
5273
5274 // table that indicates if an edge is cut (if the i-th bit is set the i-th
5275 // line is cut)
5276 constexpr std::array<unsigned int, n_configurations> cut_line_table = {
5277 {0b0000,
5278 0b0101,
5279 0b0110,
5280 0b0011,
5281 0b1010,
5282 0b0000,
5283 0b1100,
5284 0b1001,
5285 0b1001,
5286 0b1100,
5287 0b0000,
5288 0b1010,
5289 0b0011,
5290 0b0110,
5291 0b0101,
5292 0b0000}};
5293
5294 // list of the definition of the newly created lines (each line is defined
5295 // by two edges it cuts)
5296 constexpr ndarray<unsigned int, n_configurations, 5> new_line_table = {
5297 {{{X, X, X, X, X}},
5298 {{0, 2, X, X, X}},
5299 {{1, 2, X, X, X}},
5300 {{0, 1, X, X, X}},
5301 {{1, 3, X, X, X}},
5302 {{X, X, X, X, X}},
5303 {{2, 3, X, X, X}},
5304 {{0, 3, X, X, X}},
5305 {{0, 3, X, X, X}},
5306 {{2, 3, X, X, X}},
5307 {{X, X, X, X, X}},
5308 {{1, 3, X, X, X}},
5309 {{0, 1, X, X, X}},
5310 {{2, 1, X, X, X}},
5311 {{0, 2, X, X, X}},
5312 {{X, X, X, X, X}}}};
5313
5314 // vertices of each line
5315 constexpr ndarray<unsigned int, n_lines, 2> line_to_vertex_table = {
5316 {{{0, 3}}, {{1, 2}}, {{0, 1}}, {{3, 2}}}};
5317
5318 // run dimension-independent code
5320 n_vertices,
5321 n_sub_vertices,
5322 n_configurations,
5323 n_lines,
5324 5>(cut_line_table,
5325 new_line_table,
5326 line_to_vertex_table,
5327 ls_values,
5328 points,
5329 mask,
5330 iso_level,
5331 tolerance,
5332 vertices,
5333 cells,
5334 write_back_cell_data);
5335 }
5336
5337
5338
5339 template <int dim, typename VectorType>
5340 void
5342 const std::vector<value_type> &ls_values,
5343 const std::vector<Point<3>> &points,
5344 const std::vector<unsigned int> &mask,
5345 const double iso_level,
5346 std::vector<Point<3>> &vertices,
5347 std::vector<CellData<2>> &cells,
5348 const bool write_back_cell_data) const
5349 {
5350 // set up dimension-dependent sizes and tables
5351 constexpr unsigned int n_vertices = 8;
5352 constexpr unsigned int n_sub_vertices = 3;
5353 constexpr unsigned int n_lines = 12;
5354 constexpr unsigned int n_configurations = Utilities::pow(2, n_vertices);
5355 constexpr unsigned int X = static_cast<unsigned int>(-1);
5356
5357 // clang-format off
5358 // table that indicates if an edge is cut (if the i-th bit is set the i-th
5359 // line is cut)
5360 constexpr std::array<unsigned int, n_configurations> cut_line_table = {{
5361 0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905,
5362 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 0x190, 0x99, 0x393, 0x29a,
5363 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93,
5364 0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, 0x53c,
5365 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9,
5366 0x1a3, 0xaa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6,
5367 0xfaa, 0xea3, 0xda9, 0xca0, 0x460, 0x569, 0x663, 0x76a, 0x66, 0x16f,
5368 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
5369 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5,
5370 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a,
5371 0x256, 0x35f, 0x55, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53,
5372 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc,
5373 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9,
5374 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6,
5375 0x4ca, 0x5c3, 0x6c9, 0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f,
5376 0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
5377 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5,
5378 0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a,
5379 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663,
5380 0x569, 0x460, 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
5381 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, 0xd30, 0xc39,
5382 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636,
5383 0x13a, 0x33, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f,
5384 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190,
5385 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605,
5386 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0}};
5387 // clang-format on
5388
5389 // list of the definition of the newly created triangles (each triangles is
5390 // defined by two edges it cuts)
5391 constexpr ndarray<unsigned int, n_configurations, 16> new_line_table = {
5392 {{{X, X, X, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5393 {{0, 8, 3, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5394 {{0, 1, 9, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5395 {{1, 8, 3, 9, 8, 1, X, X, X, X, X, X, X, X, X, X}},
5396 {{1, 2, 10, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5397 {{0, 8, 3, 1, 2, 10, X, X, X, X, X, X, X, X, X, X}},
5398 {{9, 2, 10, 0, 2, 9, X, X, X, X, X, X, X, X, X, X}},
5399 {{2, 8, 3, 2, 10, 8, 10, 9, 8, X, X, X, X, X, X, X}},
5400 {{3, 11, 2, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5401 {{0, 11, 2, 8, 11, 0, X, X, X, X, X, X, X, X, X, X}},
5402 {{1, 9, 0, 2, 3, 11, X, X, X, X, X, X, X, X, X, X}},
5403 {{1, 11, 2, 1, 9, 11, 9, 8, 11, X, X, X, X, X, X, X}},
5404 {{3, 10, 1, 11, 10, 3, X, X, X, X, X, X, X, X, X, X}},
5405 {{0, 10, 1, 0, 8, 10, 8, 11, 10, X, X, X, X, X, X, X}},
5406 {{3, 9, 0, 3, 11, 9, 11, 10, 9, X, X, X, X, X, X, X}},
5407 {{9, 8, 10, 10, 8, 11, X, X, X, X, X, X, X, X, X, X}},
5408 {{4, 7, 8, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5409 {{4, 3, 0, 7, 3, 4, X, X, X, X, X, X, X, X, X, X}},
5410 {{0, 1, 9, 8, 4, 7, X, X, X, X, X, X, X, X, X, X}},
5411 {{4, 1, 9, 4, 7, 1, 7, 3, 1, X, X, X, X, X, X, X}},
5412 {{1, 2, 10, 8, 4, 7, X, X, X, X, X, X, X, X, X, X}},
5413 {{3, 4, 7, 3, 0, 4, 1, 2, 10, X, X, X, X, X, X, X}},
5414 {{9, 2, 10, 9, 0, 2, 8, 4, 7, X, X, X, X, X, X, X}},
5415 {{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, X, X, X, X}},
5416 {{8, 4, 7, 3, 11, 2, X, X, X, X, X, X, X, X, X, X}},
5417 {{11, 4, 7, 11, 2, 4, 2, 0, 4, X, X, X, X, X, X, X}},
5418 {{9, 0, 1, 8, 4, 7, 2, 3, 11, X, X, X, X, X, X, X}},
5419 {{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, X, X, X, X}},
5420 {{3, 10, 1, 3, 11, 10, 7, 8, 4, X, X, X, X, X, X, X}},
5421 {{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, X, X, X, X}},
5422 {{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, X, X, X, X}},
5423 {{4, 7, 11, 4, 11, 9, 9, 11, 10, X, X, X, X, X, X, X}},
5424 {{9, 5, 4, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5425 {{9, 5, 4, 0, 8, 3, X, X, X, X, X, X, X, X, X, X}},
5426 {{0, 5, 4, 1, 5, 0, X, X, X, X, X, X, X, X, X, X}},
5427 {{8, 5, 4, 8, 3, 5, 3, 1, 5, X, X, X, X, X, X, X}},
5428 {{1, 2, 10, 9, 5, 4, X, X, X, X, X, X, X, X, X, X}},
5429 {{3, 0, 8, 1, 2, 10, 4, 9, 5, X, X, X, X, X, X, X}},
5430 {{5, 2, 10, 5, 4, 2, 4, 0, 2, X, X, X, X, X, X, X}},
5431 {{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, X, X, X, X}},
5432 {{9, 5, 4, 2, 3, 11, X, X, X, X, X, X, X, X, X, X}},
5433 {{0, 11, 2, 0, 8, 11, 4, 9, 5, X, X, X, X, X, X, X}},
5434 {{0, 5, 4, 0, 1, 5, 2, 3, 11, X, X, X, X, X, X, X}},
5435 {{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, X, X, X, X}},
5436 {{10, 3, 11, 10, 1, 3, 9, 5, 4, X, X, X, X, X, X, X}},
5437 {{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, X, X, X, X}},
5438 {{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, X, X, X, X}},
5439 {{5, 4, 8, 5, 8, 10, 10, 8, 11, X, X, X, X, X, X, X}},
5440 {{9, 7, 8, 5, 7, 9, X, X, X, X, X, X, X, X, X, X}},
5441 {{9, 3, 0, 9, 5, 3, 5, 7, 3, X, X, X, X, X, X, X}},
5442 {{0, 7, 8, 0, 1, 7, 1, 5, 7, X, X, X, X, X, X, X}},
5443 {{1, 5, 3, 3, 5, 7, X, X, X, X, X, X, X, X, X, X}},
5444 {{9, 7, 8, 9, 5, 7, 10, 1, 2, X, X, X, X, X, X, X}},
5445 {{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, X, X, X, X}},
5446 {{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, X, X, X, X}},
5447 {{2, 10, 5, 2, 5, 3, 3, 5, 7, X, X, X, X, X, X, X}},
5448 {{7, 9, 5, 7, 8, 9, 3, 11, 2, X, X, X, X, X, X, X}},
5449 {{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, X, X, X, X}},
5450 {{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, X, X, X, X}},
5451 {{11, 2, 1, 11, 1, 7, 7, 1, 5, X, X, X, X, X, X, X}},
5452 {{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, X, X, X, X}},
5453 {{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, X}},
5454 {{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, X}},
5455 {{11, 10, 5, 7, 11, 5, X, X, X, X, X, X, X, X, X, X}},
5456 {{10, 6, 5, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5457 {{0, 8, 3, 5, 10, 6, X, X, X, X, X, X, X, X, X, X}},
5458 {{9, 0, 1, 5, 10, 6, X, X, X, X, X, X, X, X, X, X}},
5459 {{1, 8, 3, 1, 9, 8, 5, 10, 6, X, X, X, X, X, X, X}},
5460 {{1, 6, 5, 2, 6, 1, X, X, X, X, X, X, X, X, X, X}},
5461 {{1, 6, 5, 1, 2, 6, 3, 0, 8, X, X, X, X, X, X, X}},
5462 {{9, 6, 5, 9, 0, 6, 0, 2, 6, X, X, X, X, X, X, X}},
5463 {{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, X, X, X, X}},
5464 {{2, 3, 11, 10, 6, 5, X, X, X, X, X, X, X, X, X, X}},
5465 {{11, 0, 8, 11, 2, 0, 10, 6, 5, X, X, X, X, X, X, X}},
5466 {{0, 1, 9, 2, 3, 11, 5, 10, 6, X, X, X, X, X, X, X}},
5467 {{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, X, X, X, X}},
5468 {{6, 3, 11, 6, 5, 3, 5, 1, 3, X, X, X, X, X, X, X}},
5469 {{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, X, X, X, X}},
5470 {{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, X, X, X, X}},
5471 {{6, 5, 9, 6, 9, 11, 11, 9, 8, X, X, X, X, X, X, X}},
5472 {{5, 10, 6, 4, 7, 8, X, X, X, X, X, X, X, X, X, X}},
5473 {{4, 3, 0, 4, 7, 3, 6, 5, 10, X, X, X, X, X, X, X}},
5474 {{1, 9, 0, 5, 10, 6, 8, 4, 7, X, X, X, X, X, X, X}},
5475 {{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, X, X, X, X}},
5476 {{6, 1, 2, 6, 5, 1, 4, 7, 8, X, X, X, X, X, X, X}},
5477 {{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, X, X, X, X}},
5478 {{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, X, X, X, X}},
5479 {{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, X}},
5480 {{3, 11, 2, 7, 8, 4, 10, 6, 5, X, X, X, X, X, X, X}},
5481 {{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, X, X, X, X}},
5482 {{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, X, X, X, X}},
5483 {{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, X}},
5484 {{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, X, X, X, X}},
5485 {{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, X}},
5486 {{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, X}},
5487 {{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, X, X, X, X}},
5488 {{10, 4, 9, 6, 4, 10, X, X, X, X, X, X, X, X, X, X}},
5489 {{4, 10, 6, 4, 9, 10, 0, 8, 3, X, X, X, X, X, X, X}},
5490 {{10, 0, 1, 10, 6, 0, 6, 4, 0, X, X, X, X, X, X, X}},
5491 {{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, X, X, X, X}},
5492 {{1, 4, 9, 1, 2, 4, 2, 6, 4, X, X, X, X, X, X, X}},
5493 {{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, X, X, X, X}},
5494 {{0, 2, 4, 4, 2, 6, X, X, X, X, X, X, X, X, X, X}},
5495 {{8, 3, 2, 8, 2, 4, 4, 2, 6, X, X, X, X, X, X, X}},
5496 {{10, 4, 9, 10, 6, 4, 11, 2, 3, X, X, X, X, X, X, X}},
5497 {{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, X, X, X, X}},
5498 {{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, X, X, X, X}},
5499 {{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, X}},
5500 {{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, X, X, X, X}},
5501 {{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, X}},
5502 {{3, 11, 6, 3, 6, 0, 0, 6, 4, X, X, X, X, X, X, X}},
5503 {{6, 4, 8, 11, 6, 8, X, X, X, X, X, X, X, X, X, X}},
5504 {{7, 10, 6, 7, 8, 10, 8, 9, 10, X, X, X, X, X, X, X}},
5505 {{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, X, X, X, X}},
5506 {{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, X, X, X, X}},
5507 {{10, 6, 7, 10, 7, 1, 1, 7, 3, X, X, X, X, X, X, X}},
5508 {{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, X, X, X, X}},
5509 {{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, X}},
5510 {{7, 8, 0, 7, 0, 6, 6, 0, 2, X, X, X, X, X, X, X}},
5511 {{7, 3, 2, 6, 7, 2, X, X, X, X, X, X, X, X, X, X}},
5512 {{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, X, X, X, X}},
5513 {{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, X}},
5514 {{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, X}},
5515 {{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, X, X, X, X}},
5516 {{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, X}},
5517 {{0, 9, 1, 11, 6, 7, X, X, X, X, X, X, X, X, X, X}},
5518 {{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, X, X, X, X}},
5519 {{7, 11, 6, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5520 {{7, 6, 11, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5521 {{3, 0, 8, 11, 7, 6, X, X, X, X, X, X, X, X, X, X}},
5522 {{0, 1, 9, 11, 7, 6, X, X, X, X, X, X, X, X, X, X}},
5523 {{8, 1, 9, 8, 3, 1, 11, 7, 6, X, X, X, X, X, X, X}},
5524 {{10, 1, 2, 6, 11, 7, X, X, X, X, X, X, X, X, X, X}},
5525 {{1, 2, 10, 3, 0, 8, 6, 11, 7, X, X, X, X, X, X, X}},
5526 {{2, 9, 0, 2, 10, 9, 6, 11, 7, X, X, X, X, X, X, X}},
5527 {{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, X, X, X, X}},
5528 {{7, 2, 3, 6, 2, 7, X, X, X, X, X, X, X, X, X, X}},
5529 {{7, 0, 8, 7, 6, 0, 6, 2, 0, X, X, X, X, X, X, X}},
5530 {{2, 7, 6, 2, 3, 7, 0, 1, 9, X, X, X, X, X, X, X}},
5531 {{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, X, X, X, X}},
5532 {{10, 7, 6, 10, 1, 7, 1, 3, 7, X, X, X, X, X, X, X}},
5533 {{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, X, X, X, X}},
5534 {{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, X, X, X, X}},
5535 {{7, 6, 10, 7, 10, 8, 8, 10, 9, X, X, X, X, X, X, X}},
5536 {{6, 8, 4, 11, 8, 6, X, X, X, X, X, X, X, X, X, X}},
5537 {{3, 6, 11, 3, 0, 6, 0, 4, 6, X, X, X, X, X, X, X}},
5538 {{8, 6, 11, 8, 4, 6, 9, 0, 1, X, X, X, X, X, X, X}},
5539 {{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, X, X, X, X}},
5540 {{6, 8, 4, 6, 11, 8, 2, 10, 1, X, X, X, X, X, X, X}},
5541 {{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, X, X, X, X}},
5542 {{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, X, X, X, X}},
5543 {{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, X}},
5544 {{8, 2, 3, 8, 4, 2, 4, 6, 2, X, X, X, X, X, X, X}},
5545 {{0, 4, 2, 4, 6, 2, X, X, X, X, X, X, X, X, X, X}},
5546 {{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, X, X, X, X}},
5547 {{1, 9, 4, 1, 4, 2, 2, 4, 6, X, X, X, X, X, X, X}},
5548 {{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, X, X, X, X}},
5549 {{10, 1, 0, 10, 0, 6, 6, 0, 4, X, X, X, X, X, X, X}},
5550 {{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, X}},
5551 {{10, 9, 4, 6, 10, 4, X, X, X, X, X, X, X, X, X, X}},
5552 {{4, 9, 5, 7, 6, 11, X, X, X, X, X, X, X, X, X, X}},
5553 {{0, 8, 3, 4, 9, 5, 11, 7, 6, X, X, X, X, X, X, X}},
5554 {{5, 0, 1, 5, 4, 0, 7, 6, 11, X, X, X, X, X, X, X}},
5555 {{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, X, X, X, X}},
5556 {{9, 5, 4, 10, 1, 2, 7, 6, 11, X, X, X, X, X, X, X}},
5557 {{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, X, X, X, X}},
5558 {{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, X, X, X, X}},
5559 {{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, X}},
5560 {{7, 2, 3, 7, 6, 2, 5, 4, 9, X, X, X, X, X, X, X}},
5561 {{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, X, X, X, X}},
5562 {{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, X, X, X, X}},
5563 {{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, X}},
5564 {{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, X, X, X, X}},
5565 {{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, X}},
5566 {{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, X}},
5567 {{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, X, X, X, X}},
5568 {{6, 9, 5, 6, 11, 9, 11, 8, 9, X, X, X, X, X, X, X}},
5569 {{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, X, X, X, X}},
5570 {{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, X, X, X, X}},
5571 {{6, 11, 3, 6, 3, 5, 5, 3, 1, X, X, X, X, X, X, X}},
5572 {{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, X, X, X, X}},
5573 {{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, X}},
5574 {{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, X}},
5575 {{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, X, X, X, X}},
5576 {{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, X, X, X, X}},
5577 {{9, 5, 6, 9, 6, 0, 0, 6, 2, X, X, X, X, X, X, X}},
5578 {{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, X}},
5579 {{1, 5, 6, 2, 1, 6, X, X, X, X, X, X, X, X, X, X}},
5580 {{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, X}},
5581 {{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, X, X, X, X}},
5582 {{0, 3, 8, 5, 6, 10, X, X, X, X, X, X, X, X, X, X}},
5583 {{10, 5, 6, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5584 {{11, 5, 10, 7, 5, 11, X, X, X, X, X, X, X, X, X, X}},
5585 {{11, 5, 10, 11, 7, 5, 8, 3, 0, X, X, X, X, X, X, X}},
5586 {{5, 11, 7, 5, 10, 11, 1, 9, 0, X, X, X, X, X, X, X}},
5587 {{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, X, X, X, X}},
5588 {{11, 1, 2, 11, 7, 1, 7, 5, 1, X, X, X, X, X, X, X}},
5589 {{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, X, X, X, X}},
5590 {{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, X, X, X, X}},
5591 {{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, X}},
5592 {{2, 5, 10, 2, 3, 5, 3, 7, 5, X, X, X, X, X, X, X}},
5593 {{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, X, X, X, X}},
5594 {{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, X, X, X, X}},
5595 {{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, X}},
5596 {{1, 3, 5, 3, 7, 5, X, X, X, X, X, X, X, X, X, X}},
5597 {{0, 8, 7, 0, 7, 1, 1, 7, 5, X, X, X, X, X, X, X}},
5598 {{9, 0, 3, 9, 3, 5, 5, 3, 7, X, X, X, X, X, X, X}},
5599 {{9, 8, 7, 5, 9, 7, X, X, X, X, X, X, X, X, X, X}},
5600 {{5, 8, 4, 5, 10, 8, 10, 11, 8, X, X, X, X, X, X, X}},
5601 {{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, X, X, X, X}},
5602 {{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, X, X, X, X}},
5603 {{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, X}},
5604 {{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, X, X, X, X}},
5605 {{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, X}},
5606 {{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, X}},
5607 {{9, 4, 5, 2, 11, 3, X, X, X, X, X, X, X, X, X, X}},
5608 {{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, X, X, X, X}},
5609 {{5, 10, 2, 5, 2, 4, 4, 2, 0, X, X, X, X, X, X, X}},
5610 {{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, X}},
5611 {{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, X, X, X, X}},
5612 {{8, 4, 5, 8, 5, 3, 3, 5, 1, X, X, X, X, X, X, X}},
5613 {{0, 4, 5, 1, 0, 5, X, X, X, X, X, X, X, X, X, X}},
5614 {{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, X, X, X, X}},
5615 {{9, 4, 5, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5616 {{4, 11, 7, 4, 9, 11, 9, 10, 11, X, X, X, X, X, X, X}},
5617 {{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, X, X, X, X}},
5618 {{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, X, X, X, X}},
5619 {{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, X}},
5620 {{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, X, X, X, X}},
5621 {{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, X}},
5622 {{11, 7, 4, 11, 4, 2, 2, 4, 0, X, X, X, X, X, X, X}},
5623 {{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, X, X, X, X}},
5624 {{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, X, X, X, X}},
5625 {{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, X}},
5626 {{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, X}},
5627 {{1, 10, 2, 8, 7, 4, X, X, X, X, X, X, X, X, X, X}},
5628 {{4, 9, 1, 4, 1, 7, 7, 1, 3, X, X, X, X, X, X, X}},
5629 {{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, X, X, X, X}},
5630 {{4, 0, 3, 7, 4, 3, X, X, X, X, X, X, X, X, X, X}},
5631 {{4, 8, 7, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5632 {{9, 10, 8, 10, 11, 8, X, X, X, X, X, X, X, X, X, X}},
5633 {{3, 0, 9, 3, 9, 11, 11, 9, 10, X, X, X, X, X, X, X}},
5634 {{0, 1, 10, 0, 10, 8, 8, 10, 11, X, X, X, X, X, X, X}},
5635 {{3, 1, 10, 11, 3, 10, X, X, X, X, X, X, X, X, X, X}},
5636 {{1, 2, 11, 1, 11, 9, 9, 11, 8, X, X, X, X, X, X, X}},
5637 {{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, X, X, X, X}},
5638 {{0, 2, 11, 8, 0, 11, X, X, X, X, X, X, X, X, X, X}},
5639 {{3, 2, 11, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5640 {{2, 3, 8, 2, 8, 10, 10, 8, 9, X, X, X, X, X, X, X}},
5641 {{9, 10, 2, 0, 9, 2, X, X, X, X, X, X, X, X, X, X}},
5642 {{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, X, X, X, X}},
5643 {{1, 10, 2, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5644 {{1, 3, 8, 9, 1, 8, X, X, X, X, X, X, X, X, X, X}},
5645 {{0, 9, 1, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5646 {{0, 3, 8, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5647 {{X, X, X, X, X, X, X, X, X, X, X, X, X, X, X, X}}}};
5648
5649 // vertices of each line
5650 static constexpr ndarray<unsigned int, n_lines, 2> line_to_vertex_table = {
5651 {{{0, 1}},
5652 {{1, 2}},
5653 {{2, 3}},
5654 {{3, 0}},
5655 {{4, 5}},
5656 {{5, 6}},
5657 {{6, 7}},
5658 {{7, 4}},
5659 {{0, 4}},
5660 {{1, 5}},
5661 {{2, 6}},
5662 {{3, 7}}}};
5663
5664 // run dimension-independent code
5666 n_vertices,
5667 n_sub_vertices,
5668 n_configurations,
5669 n_lines,
5670 16>(cut_line_table,
5671 new_line_table,
5672 line_to_vertex_table,
5673 ls_values,
5674 points,
5675 mask,
5676 iso_level,
5677 tolerance,
5678 vertices,
5679 cells,
5680 write_back_cell_data);
5681 }
5682
5683} /* namespace GridTools */
5684
5685
5686// explicit instantiations
5687#include "grid_tools.inst"
5688
void distribute(VectorType &vec) const
std::pair< std::vector< std::pair< int, int > >, std::vector< int > > query(const QueryType &queries)
DistributedTree(const MPI_Comm comm, const std::vector< BoundingBox< dim, Number > > &bounding_boxes)
BoundingBox< spacedim, Number > create_extended_relative(const Number relative_amount) const
BoundingBox< spacedim, Number > create_extended(const Number amount) const
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
types::global_dof_index n_dofs() const
Definition fe_q.h:550
const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_vertex_to_cell_map() const
const std::vector< std::vector< Tensor< 1, spacedim > > > & get_vertex_to_cell_centers_directions() const
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_locally_owned_cell_bounding_boxes_rtree() const
const Mapping< dim, spacedim > & get_mapping() const
const Triangulation< dim, spacedim > & get_triangulation() const
const RTree< std::pair< Point< spacedim >, unsigned int > > & get_used_vertices_rtree() const
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_cell_bounding_boxes_rtree() const
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
void process_cell(const typename DoFHandler< dim >::active_cell_iterator &cell, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim > > &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 > > &cells) const
void process_sub_cell(const std::vector< value_type > &, const std::vector< Point< 1 > > &, const std::vector< unsigned int > &, const double, std::vector< Point< 1 > > &, std::vector< CellData< 1 > > &, const bool) const
void process(const DoFHandler< dim > &background_dof_handler, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim > > &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 > > &cells) const
const Tensor< 2, 2, double > rotation_matrix
Point< 2 > operator()(const Point< 2 > &p) const
Rotate2d(const double angle)
Point< 3 > operator()(const Point< 3 > &p) const
Rotate3d(const Tensor< 1, 3, double > &axis, const double angle)
const Tensor< 2, 3, double > rotation_matrix
Scale(const double factor)
Point< spacedim > operator()(const Point< spacedim > p) const
Shift(const Tensor< 1, spacedim > &shift)
Point< spacedim > operator()(const Point< spacedim > p) const
const Tensor< 1, spacedim > shift
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const =0
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
constexpr numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr numbers::NumberTraits< Number >::real_type square() const
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
Quadrature< spacedim > compute_affine_transformation(const std::array< Point< spacedim >, dim+1 > &vertices) const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
size_type n() const
size_type n_rows() const
size_type n_cols() const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
numbers::NumberTraits< Number >::real_type norm() const
IteratorState::IteratorStates state() const
virtual MPI_Comm get_communicator() const
unsigned int n_quads() const
void load_user_indices(const std::vector< unsigned int > &v)
virtual void clear()
bool all_reference_cells_are_hyper_cube() const
void clear_user_data()
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, unsigned char > > & get_periodic_face_map() const
face_iterator end_face() const
cell_iterator create_cell_iterator(const CellId &cell_id) const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_lines() const
unsigned int n_raw_faces() const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int n_levels() const
cell_iterator end() const
virtual bool has_hanging_nodes() const
virtual void execute_coarsening_and_refinement()
unsigned int n_cells() const
const std::vector< bool > & get_used_vertices() const
Triangulation< dim, spacedim > & get_triangulation()
Signals signals
Definition tria.h:2509
unsigned int n_vertices() const
void save_user_indices(std::vector< unsigned int > &v) const
virtual std::vector< types::boundary_id > get_boundary_ids() const
active_face_iterator begin_active_face() const
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]