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