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