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