PolyDEAL
 
Loading...
Searching...
No Matches
poly_utils.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) XXXX - YYYY by the polyDEAL authors
5//
6// This file is part of the polyDEAL library.
7//
8// Detailed license information governing the source code
9// can be found in LICENSE.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef poly_utils_h
15#define poly_utils_h
16
17
18#include <deal.II/base/config.h>
19
20#include <deal.II/base/conditional_ostream.h>
21#include <deal.II/base/point.h>
22#include <deal.II/base/quadrature.h>
23#include <deal.II/base/std_cxx20/iota_view.h>
24
25#include <deal.II/boost_adaptors/bounding_box.h>
26#include <deal.II/boost_adaptors/point.h>
27#include <deal.II/boost_adaptors/segment.h>
28
29#include <deal.II/distributed/tria.h>
30
31#include <deal.II/dofs/dof_handler.h>
32
33#include <deal.II/fe/fe_dgq.h>
34#include <deal.II/fe/fe_values.h>
35
36#include <deal.II/grid/grid_tools.h>
37
38#include <deal.II/lac/dynamic_sparsity_pattern.h>
39#include <deal.II/lac/sparse_matrix.h>
40#include <deal.II/lac/sparsity_pattern.h>
41#include <deal.II/lac/sparsity_tools.h>
42#include <deal.II/lac/trilinos_sparse_matrix.h>
43
44#include <deal.II/numerics/vector_tools_common.h>
45
46#include <boost/geometry/algorithms/distance.hpp>
47#include <boost/geometry/index/detail/rtree/utilities/print.hpp>
48#include <boost/geometry/index/rtree.hpp>
49#include <boost/geometry/strategies/strategies.hpp>
50
51#include <deal.II/cgal/point_conversion.h>
52
53#ifdef DEAL_II_WITH_TRILINOS
54# include <EpetraExt_RowMatrixOut.h>
55#endif
56
57#ifdef DEAL_II_WITH_CGAL
58
59# include <CGAL/Constrained_Delaunay_triangulation_2.h>
60# include <CGAL/Constrained_triangulation_plus_2.h>
61# include <CGAL/Exact_predicates_exact_constructions_kernel.h>
62# include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
63# include <CGAL/Polygon_2.h>
64# include <CGAL/Polygon_with_holes_2.h>
65# include <CGAL/Segment_Delaunay_graph_2.h>
66# include <CGAL/Segment_Delaunay_graph_traits_2.h>
67# include <CGAL/intersections.h>
68# include <CGAL/squared_distance_2.h>
69# include <CGAL/squared_distance_3.h>
70
71
72#endif
73
74#include <memory>
75
76
77namespace dealii::PolyUtils::internal
78{
82 inline types::global_cell_index
83 get_index(const std::vector<types::global_cell_index> &v,
84 const types::global_cell_index index)
85 {
86 return std::distance(v.begin(), std::find(v.begin(), v.end(), index));
87 }
88
89
90
95 template <int dim, int spacedim>
96 void
97 get_face_connectivity_of_cells(
98 const parallel::fullydistributed::Triangulation<dim, spacedim>
99 &triangulation,
100 DynamicSparsityPattern &cell_connectivity,
101 const std::vector<types::global_cell_index> locally_owned_cells)
102 {
103 cell_connectivity.reinit(triangulation.n_locally_owned_active_cells(),
104 triangulation.n_locally_owned_active_cells());
105
106
107 // loop over all cells and their neighbors to build the sparsity
108 // pattern. note that it's a bit hard to enter all the connections when
109 // a neighbor has children since we would need to find out which of its
110 // children is adjacent to the current cell. this problem can be omitted
111 // if we only do something if the neighbor has no children -- in that
112 // case it is either on the same or a coarser level than we are. in
113 // return, we have to add entries in both directions for both cells
114 for (const auto &cell : triangulation.active_cell_iterators())
115 {
116 if (cell->is_locally_owned())
117 {
118 const unsigned int index = cell->active_cell_index();
119 cell_connectivity.add(get_index(locally_owned_cells, index),
120 get_index(locally_owned_cells, index));
121 for (auto f : cell->face_indices())
122 if ((cell->at_boundary(f) == false) &&
123 (cell->neighbor(f)->has_children() == false) &&
124 cell->neighbor(f)->is_locally_owned())
125 {
126 const unsigned int other_index =
127 cell->neighbor(f)->active_cell_index();
128
129 cell_connectivity.add(get_index(locally_owned_cells, index),
130 get_index(locally_owned_cells,
131 other_index));
132 cell_connectivity.add(get_index(locally_owned_cells,
133 other_index),
134 get_index(locally_owned_cells, index));
135 }
136 }
137 }
138 }
139} // namespace dealii::PolyUtils::internal
140
141
142namespace dealii::PolyUtils
143{
144 template <typename Value,
145 typename Options,
146 typename Translator,
147 typename Box,
148 typename Allocators>
149 struct Rtree_visitor : public boost::geometry::index::detail::rtree::visitor<
150 Value,
151 typename Options::parameters_type,
152 Box,
153 Allocators,
154 typename Options::node_tag,
155 true>::type
156 {
157 inline Rtree_visitor(
158 const Translator &translator,
159 unsigned int target_level,
160 std::vector<std::vector<typename Triangulation<
161 boost::geometry::dimension<Box>::value>::active_cell_iterator>> &boxes,
162 std::vector<std::vector<unsigned int>> &csr);
163
164
168 using InternalNode =
169 typename boost::geometry::index::detail::rtree::internal_node<
170 Value,
171 typename Options::parameters_type,
172 Box,
173 Allocators,
174 typename Options::node_tag>::type;
175
179 using Leaf = typename boost::geometry::index::detail::rtree::leaf<
180 Value,
181 typename Options::parameters_type,
182 Box,
183 Allocators,
184 typename Options::node_tag>::type;
185
190 inline void
191 operator()(const InternalNode &node);
192
196 inline void
197 operator()(const Leaf &);
198
202 const Translator &translator;
203
207 size_t level;
208
213 size_t node_counter;
214
215 size_t next_level_leafs_processed;
220 const size_t target_level;
221
227 std::vector<std::vector<typename Triangulation<
228 boost::geometry::dimension<Box>::value>::active_cell_iterator>>
229 &agglomerates;
230
231 std::vector<std::vector<unsigned int>> &row_ptr;
232 };
233
234
235
236 template <typename Value,
237 typename Options,
238 typename Translator,
239 typename Box,
240 typename Allocators>
241 Rtree_visitor<Value, Options, Translator, Box, Allocators>::Rtree_visitor(
242 const Translator &translator,
243 const unsigned int target_level,
244 std::vector<std::vector<typename Triangulation<
245 boost::geometry::dimension<Box>::value>::active_cell_iterator>>
246 &bb_in_boxes,
247 std::vector<std::vector<unsigned int>> &csr)
248 : translator(translator)
249 , level(0)
250 , node_counter(0)
251 , next_level_leafs_processed(0)
252 , target_level(target_level)
253 , agglomerates(bb_in_boxes)
254 , row_ptr(csr)
255 {}
256
257
258
259 template <typename Value,
260 typename Options,
261 typename Translator,
262 typename Box,
263 typename Allocators>
264 void
265 Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()(
266 const Rtree_visitor::InternalNode &node)
267 {
268 using elements_type =
269 typename boost::geometry::index::detail::rtree::elements_type<
270 InternalNode>::type; // pairs of bounding box and pointer to child node
271 const elements_type &elements =
272 boost::geometry::index::detail::rtree::elements(node);
273
274 if (level < target_level)
275 {
276 size_t level_backup = level;
277 ++level;
278
279 for (typename elements_type::const_iterator it = elements.begin();
280 it != elements.end();
281 ++it)
282 {
283 boost::geometry::index::detail::rtree::apply_visitor(*this,
284 *it->second);
285 }
286
287 level = level_backup;
288 }
289 else if (level == target_level)
290 {
291 // const unsigned int n_children = elements.size();
292 const auto offset = agglomerates.size();
293 agglomerates.resize(offset + 1);
294 row_ptr.resize(row_ptr.size() + 1);
295 next_level_leafs_processed = 0;
296 row_ptr.back().push_back(
297 next_level_leafs_processed); // convention: row_ptr[0]=0
298 size_t level_backup = level;
299
300 ++level;
301 for (const auto &child : elements)
302 {
303 boost::geometry::index::detail::rtree::apply_visitor(*this,
304 *child.second);
305 }
306 // Done with node number 'node_counter'
307
308 ++node_counter; // visited all children of an internal node
309
310 level = level_backup;
311 }
312 else if (level > target_level)
313 {
314 // Keep visiting until you go to the leafs.
315 size_t level_backup = level;
316
317 ++level;
318
319 for (const auto &child : elements)
320 {
321 boost::geometry::index::detail::rtree::apply_visitor(*this,
322 *child.second);
323 }
324 level = level_backup;
325 row_ptr[node_counter].push_back(next_level_leafs_processed);
326 }
327 }
328
329
330
331 template <typename Value,
332 typename Options,
333 typename Translator,
334 typename Box,
335 typename Allocators>
336 void
337 Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()(
338 const Rtree_visitor::Leaf &leaf)
339 {
340 using elements_type =
341 typename boost::geometry::index::detail::rtree::elements_type<
342 Leaf>::type; // pairs of bounding box and pointer to child node
343 const elements_type &elements =
344 boost::geometry::index::detail::rtree::elements(leaf);
345
346
347 for (const auto &it : elements)
348 {
349 agglomerates[node_counter].push_back(it.second);
350 }
351 next_level_leafs_processed += elements.size();
352 }
353
354 template <typename Rtree>
355 inline std::pair<
356 std::vector<std::vector<unsigned int>>,
357 std::vector<std::vector<typename Triangulation<boost::geometry::dimension<
358 typename Rtree::indexable_type>::value>::active_cell_iterator>>>
359 extract_children_of_level(const Rtree &tree, const unsigned int level)
360 {
361 using RtreeView =
362 boost::geometry::index::detail::rtree::utilities::view<Rtree>;
363 RtreeView rtv(tree);
364
365 std::vector<std::vector<unsigned int>> csrs;
366 std::vector<std::vector<typename Triangulation<boost::geometry::dimension<
367 typename Rtree::indexable_type>::value>::active_cell_iterator>>
368 agglomerates;
369
370 if (rtv.depth() == 0)
371 {
372 // The below algorithm does not work for `rtv.depth()==0`, which might
373 // happen if the number entries in the tree is too small.
374 // In this case, simply return a single bounding box.
375 agglomerates.resize(1);
376 agglomerates[0].resize(1);
377 csrs.resize(1);
378 csrs[0].resize(1);
379 }
380 else
381 {
382 const unsigned int target_level =
383 std::min<unsigned int>(level, rtv.depth());
384
385 Rtree_visitor<typename RtreeView::value_type,
386 typename RtreeView::options_type,
387 typename RtreeView::translator_type,
388 typename RtreeView::box_type,
389 typename RtreeView::allocators_type>
390 node_visitor(rtv.translator(), target_level, agglomerates, csrs);
391 rtv.apply_visitor(node_visitor);
392 }
393 AssertDimension(agglomerates.size(), csrs.size());
394
395 return {csrs, agglomerates};
396 }
397
398
399 template <int dim, typename Number = double>
400 Number
401 compute_h_orthogonal(
402 const unsigned int face_index,
403 const std::vector<typename Triangulation<dim>::active_face_iterator>
404 &polygon_boundary,
405 const Tensor<1, dim> &deal_normal)
406 {
407#ifdef DEAL_II_WITH_CGAL
408
409 using Kernel = CGAL::Exact_predicates_exact_constructions_kernel;
410 std::vector<typename Kernel::FT> candidates;
411 candidates.reserve(polygon_boundary.size() - 1);
412
413 // Initialize the range of faces to be checked for intersection: they are
414 // {0,..,n_faces-1}\setminus the current face index face_index.
415 std::vector<unsigned int> face_indices(polygon_boundary.size());
416 std::iota(face_indices.begin(), face_indices.end(), 0); // fill the range
417 face_indices.erase(face_indices.cbegin() +
418 face_index); // remove current index
419
420 if constexpr (dim == 2)
421 {
422 typename Kernel::Segment_2 face_segm(
423 {polygon_boundary[face_index]->vertex(0)[0],
424 polygon_boundary[face_index]->vertex(0)[1]},
425 {polygon_boundary[face_index]->vertex(1)[0],
426 polygon_boundary[face_index]->vertex(1)[1]});
427
428 // Shoot a ray from the midpoint of the face in the orthogonal direction
429 // given by deal.II normals
430 const auto &midpoint = CGAL::midpoint(face_segm);
431 // deal.II normal is always outward, flip the direction
432 const typename Kernel::Vector_2 orthogonal_direction{-deal_normal[0],
433 -deal_normal[1]};
434 const typename Kernel::Ray_2 ray(midpoint, orthogonal_direction);
435 for (const auto f : face_indices)
436 {
437 typename Kernel::Segment_2 segm({polygon_boundary[f]->vertex(0)[0],
438 polygon_boundary[f]->vertex(0)[1]},
439 {polygon_boundary[f]->vertex(1)[0],
440 polygon_boundary[f]->vertex(
441 1)[1]});
442
443 if (CGAL::do_intersect(ray, segm))
444 candidates.push_back(CGAL::squared_distance(midpoint, segm));
445 }
446 return std::sqrt(CGAL::to_double(
447 *std::min_element(candidates.cbegin(), candidates.cend())));
448 }
449 else if constexpr (dim == 3)
450 {
451 const typename Kernel::Point_3 &center{
452 polygon_boundary[face_index]->center()[0],
453 polygon_boundary[face_index]->center()[1],
454 polygon_boundary[face_index]->center()[2]};
455 // deal.II normal is always outward, flip the direction
456 const typename Kernel::Vector_3 orthogonal_direction{-deal_normal[0],
457 -deal_normal[1],
458 -deal_normal[2]};
459 const typename Kernel::Ray_3 ray(center, orthogonal_direction);
460
461 for (const auto f : face_indices)
462 {
463 // split the face into 2 triangles and compute distances
464 typename Kernel::Triangle_3 first_triangle(
465 {polygon_boundary[f]->vertex(0)[0],
466 polygon_boundary[f]->vertex(0)[1],
467 polygon_boundary[f]->vertex(0)[2]},
468 {polygon_boundary[f]->vertex(1)[0],
469 polygon_boundary[f]->vertex(1)[1],
470 polygon_boundary[f]->vertex(1)[2]},
471 {polygon_boundary[f]->vertex(3)[0],
472 polygon_boundary[f]->vertex(3)[1],
473 polygon_boundary[f]->vertex(3)[2]});
474 typename Kernel::Triangle_3 second_triangle(
475 {polygon_boundary[f]->vertex(0)[0],
476 polygon_boundary[f]->vertex(0)[1],
477 polygon_boundary[f]->vertex(0)[2]},
478 {polygon_boundary[f]->vertex(3)[0],
479 polygon_boundary[f]->vertex(3)[1],
480 polygon_boundary[f]->vertex(3)[2]},
481 {polygon_boundary[f]->vertex(2)[0],
482 polygon_boundary[f]->vertex(2)[1],
483 polygon_boundary[f]->vertex(2)[2]});
484
485 // compute point-triangle distance only if the orthogonal ray
486 // hits the triangle
487 if (CGAL::do_intersect(ray, first_triangle))
488 candidates.push_back(
489 CGAL::squared_distance(center, first_triangle));
490 if (CGAL::do_intersect(ray, second_triangle))
491 candidates.push_back(
492 CGAL::squared_distance(center, second_triangle));
493 }
494
495 return std::sqrt(CGAL::to_double(
496 *std::min_element(candidates.cbegin(), candidates.cend())));
497 }
498 else
499 {
500 Assert(false, ExcImpossibleInDim(dim));
501 (void)face_index;
502 (void)polygon_boundary;
503 return {};
504 }
505
506#else
507
508 Assert(false, ExcNeedsCGAL());
509 (void)face_index;
510 (void)polygon_boundary;
511 return {};
512#endif
513 }
514
515
516
521 template <int dim, int spacedim = dim>
522 void
523 collect_cells_for_agglomeration(
524 const Triangulation<dim, spacedim> &tria,
525 const std::vector<types::global_cell_index> &cell_idxs,
526 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
527 &cells_to_be_agglomerated)
528 {
529 Assert(cells_to_be_agglomerated.size() == 0,
530 ExcMessage(
531 "The vector of cells is supposed to be filled by this function."));
532 for (const auto &cell : tria.active_cell_iterators())
533 if (std::find(cell_idxs.begin(),
534 cell_idxs.end(),
535 cell->active_cell_index()) != cell_idxs.end())
536 {
537 cells_to_be_agglomerated.push_back(cell);
538 }
539 }
540
541
542
553 template <int dim, int spacedim>
554 void
555 partition_locally_owned_regions(const unsigned int n_partitions,
556 Triangulation<dim, spacedim> &triangulation,
557 const SparsityTools::Partitioner partitioner)
558 {
559 AssertDimension(dim, spacedim);
560 Assert(n_partitions > 0,
561 ExcMessage("Invalid number of partitions, you provided " +
562 std::to_string(n_partitions)));
563
564 auto parallel_triangulation =
565 dynamic_cast<parallel::fullydistributed::Triangulation<dim, spacedim> *>(
566 &triangulation);
567 Assert(
568 (parallel_triangulation != nullptr),
569 ExcMessage(
570 "Only fully distributed triangulations are supported. If you are using"
571 "a parallel::distributed::triangulation, you must convert it to a fully"
572 "distributed as explained in the documentation."));
573
574 // check for an easy return
575 if (n_partitions == 1)
576 {
577 for (const auto &cell : parallel_triangulation->active_cell_iterators())
578 if (cell->is_locally_owned())
579 cell->set_material_id(0);
580 return;
581 }
582
583 // collect all locally owned cells
584 std::vector<types::global_cell_index> locally_owned_cells;
585 for (const auto &cell : triangulation.active_cell_iterators())
586 if (cell->is_locally_owned())
587 locally_owned_cells.push_back(cell->active_cell_index());
588
589 DynamicSparsityPattern cell_connectivity;
590 internal::get_face_connectivity_of_cells(*parallel_triangulation,
591 cell_connectivity,
592 locally_owned_cells);
593
594 SparsityPattern sp_cell_connectivity;
595 sp_cell_connectivity.copy_from(cell_connectivity);
596
597 // partition each locally owned connection graph and get
598 // back a vector of indices, one per degree
599 // of freedom (which is associated with a
600 // cell)
601 std::vector<unsigned int> partition_indices(
602 parallel_triangulation->n_locally_owned_active_cells());
603 SparsityTools::partition(sp_cell_connectivity,
604 n_partitions,
605 partition_indices,
606 partitioner);
607
608
609 // finally loop over all cells and set the material ids
610 for (const auto &cell : parallel_triangulation->active_cell_iterators())
611 if (cell->is_locally_owned())
612 cell->set_material_id(
613 partition_indices[internal::get_index(locally_owned_cells,
614 cell->active_cell_index())]);
615 }
616
617
618
629 template <int dim, int spacedim>
630 void
631 partition_locally_owned_regions(
632 AgglomerationHandler<dim> &agglomeration_handler,
633 const unsigned int n_partitions,
634 Triangulation<dim, spacedim> &triangulation,
635 const SparsityTools::Partitioner partitioner)
636 {
637 AssertDimension(dim, spacedim);
638 Assert(
639 agglomeration_handler.n_agglomerates() == 0,
640 ExcMessage(
641 "The agglomerated grid must be empty upon calling this function."));
642 Assert(n_partitions > 0,
643 ExcMessage("Invalid number of partitions, you provided " +
644 std::to_string(n_partitions)));
645
646 auto parallel_triangulation =
647 dynamic_cast<parallel::fullydistributed::Triangulation<dim, spacedim> *>(
648 &triangulation);
649 Assert(
650 (parallel_triangulation != nullptr),
651 ExcMessage(
652 "Only fully distributed triangulations are supported. If you are using"
653 "a parallel::distributed::triangulation, you must convert it to a"
654 "fully distributed as explained in the documentation."));
655
656 // check for an easy return
657 if (n_partitions == 1)
658 {
659 for (const auto &cell : parallel_triangulation->active_cell_iterators())
660 if (cell->is_locally_owned())
661 agglomeration_handler.define_agglomerate({cell});
662 return;
663 }
664
665 // collect all locally owned cells
666 std::vector<types::global_cell_index> locally_owned_cells;
667 for (const auto &cell : triangulation.active_cell_iterators())
668 if (cell->is_locally_owned())
669 locally_owned_cells.push_back(cell->active_cell_index());
670
671 DynamicSparsityPattern cell_connectivity;
672 internal::get_face_connectivity_of_cells(*parallel_triangulation,
673 cell_connectivity,
674 locally_owned_cells);
675
676 SparsityPattern sp_cell_connectivity;
677 sp_cell_connectivity.copy_from(cell_connectivity);
678
679 // partition each locally owned connection graph and get
680 // back a vector of indices, one per degree
681 // of freedom (which is associated with a
682 // cell)
683 std::vector<unsigned int> partition_indices(
684 parallel_triangulation->n_locally_owned_active_cells());
685 SparsityTools::partition(sp_cell_connectivity,
686 n_partitions,
687 partition_indices,
688 partitioner);
689
690 std::vector<std::vector<typename Triangulation<dim>::active_cell_iterator>>
691 cells_per_partion_id;
692 cells_per_partion_id.resize(n_partitions); // number of agglomerates
693
694 // finally loop over all cells and store the ones with same partition index
695 for (const auto &cell : parallel_triangulation->active_cell_iterators())
696 if (cell->is_locally_owned())
697 cells_per_partion_id[partition_indices[internal::get_index(
698 locally_owned_cells, cell->active_cell_index())]]
699 .push_back(cell);
700
701 // All the cells with the same partition index will be merged together.
702 for (unsigned int i = 0; i < n_partitions; ++i)
703 agglomeration_handler.define_agglomerate(cells_per_partion_id[i]);
704 }
705
706
707
708 template <int dim>
709 std::
710 tuple<std::vector<double>, std::vector<double>, std::vector<double>, double>
711 compute_quality_metrics(const AgglomerationHandler<dim> &ah)
712 {
713 static_assert(dim == 2); // only 2D case is implemented.
714#ifdef DEAL_II_WITH_CGAL
715 using Kernel = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt;
716 using Polygon_with_holes = typename CGAL::Polygon_with_holes_2<Kernel>;
717 using Gt = typename CGAL::Segment_Delaunay_graph_traits_2<Kernel>;
718 using SDG2 = typename CGAL::Segment_Delaunay_graph_2<Gt>;
719 using CDT = typename CGAL::Constrained_Delaunay_triangulation_2<Kernel>;
720 using CDTP = typename CGAL::Constrained_triangulation_plus_2<CDT>;
721 using Point = typename CDTP::Point;
722 using Cid = typename CDTP::Constraint_id;
723 using Vertex_handle = typename CDTP::Vertex_handle;
724
725
726 const auto compute_radius_inscribed_circle =
727 [](const CGAL::Polygon_2<Kernel> &polygon) -> double {
728 SDG2 sdg;
729
730 sdg.insert_segments(polygon.edges_begin(), polygon.edges_end());
731
732 double sd = 0, sqdist = 0;
733 typename SDG2::Finite_faces_iterator fit = sdg.finite_faces_begin();
734 for (; fit != sdg.finite_faces_end(); ++fit)
735 {
736 typename Kernel::Point_2 pp = sdg.primal(fit);
737 for (int i = 0; i < 3; ++i)
738 {
739 assert(!sdg.is_infinite(fit->vertex(i)));
740 if (fit->vertex(i)->site().is_segment())
741 {
742 typename Kernel::Segment_2 s =
743 fit->vertex(i)->site().segment();
744 sqdist = CGAL::to_double(CGAL::squared_distance(pp, s));
745 }
746 else
747 {
748 typename Kernel::Point_2 p = fit->vertex(i)->site().point();
749 sqdist = CGAL::to_double(CGAL::squared_distance(pp, p));
750 }
751 }
752
753 if (polygon.bounded_side(pp) == CGAL::ON_BOUNDED_SIDE)
754 sd = std::max(sqdist, sd);
755 }
756
757 return std::sqrt(sd);
758 };
759
760 const auto mesh_size = [&ah]() -> double {
761 double hmax = 0.;
762 for (const auto &polytope : ah.polytope_iterators())
763 if (polytope->is_locally_owned())
764 {
765 const double diameter = polytope->diameter();
766 if (diameter > hmax)
767 hmax = diameter;
768 }
769 return hmax;
770 }();
771
772
773 // vectors holding quality metrics
774
775 // ration between radius of radius_inscribed_circle and circumscribed circle
776 std::vector<double> circle_ratios;
777 std::vector<double> unformity_factors; // diameter of element over mesh size
778 std::vector<double>
779 box_ratio; // ratio between measure of bbox and measure of element.
780
781 const std::vector<BoundingBox<dim>> &bboxes = ah.get_local_bboxes();
782 // Loop over all polytopes and compute metrics.
783 for (const auto &polytope : ah.polytope_iterators())
784 {
785 if (polytope->is_locally_owned())
786 {
787 const std::vector<typename Triangulation<dim>::active_face_iterator>
788 &boundary = polytope->polytope_boundary();
789
790 const double diameter = polytope->diameter();
791 const double radius_circumscribed_circle = .5 * diameter;
792
793 CDTP cdtp;
794 for (unsigned int f = 0; f < boundary.size(); f += 1)
795 {
796 // polyline
797 cdtp.insert_constraint(
798 {boundary[f]->vertex(0)[0], boundary[f]->vertex(0)[1]},
799 {boundary[f]->vertex(1)[0], boundary[f]->vertex(1)[1]});
800 }
801 cdtp.split_subconstraint_graph_into_constraints();
802
803 CGAL::Polygon_2<Kernel> outer_polygon;
804 auto it = outer_polygon.vertices_begin();
805 for (typename CDTP::Constraint_id cid : cdtp.constraints())
806 {
807 for (typename CDTP::Vertex_handle vh :
808 cdtp.vertices_in_constraint(cid))
809 {
810 it = outer_polygon.insert(outer_polygon.vertices_end(),
811 vh->point());
812 }
813 }
814 outer_polygon.erase(it); // remove duplicate final point
815
816 const double radius_inscribed_circle =
817 compute_radius_inscribed_circle(outer_polygon);
818
819 circle_ratios.push_back(radius_inscribed_circle /
820 radius_circumscribed_circle);
821 unformity_factors.push_back(diameter / mesh_size);
822
823 // box_ratio
824
825 const auto &agglo_values = ah.reinit(polytope);
826 const double measure_element =
827 std::accumulate(agglo_values.get_JxW_values().cbegin(),
828 agglo_values.get_JxW_values().cend(),
829 0.);
830 box_ratio.push_back(measure_element /
831 bboxes[polytope->index()].volume());
832 }
833 }
834
835
836
837 // Get all of the local bounding boxes
838 double covering_bboxes = 0.;
839 for (unsigned int i = 0; i < bboxes.size(); ++i)
840 covering_bboxes += bboxes[i].volume();
841
842 const double overlap_factor =
843 Utilities::MPI::sum(covering_bboxes,
844 ah.get_dof_handler().get_communicator()) /
845 GridTools::volume(ah.get_triangulation()); // assuming a linear mapping
846
847
848
849 return {unformity_factors, circle_ratios, box_ratio, overlap_factor};
850#else
851
852 (void)ah;
853 return {};
854#endif
855 }
856
857
861 template <int dim>
862 void
863 export_polygon_to_csv_file(
864 const AgglomerationHandler<dim> &agglomeration_handler,
865 const std::string &filename)
866 {
867 static_assert(dim == 2); // With 3D, Paraview is much better
868 std::ofstream myfile;
869 myfile.open(filename + ".csv");
870
871 for (const auto &polytope : agglomeration_handler.polytope_iterators())
872 if (polytope->is_locally_owned())
873 {
874 const std::vector<typename Triangulation<dim>::active_face_iterator>
875 &boundary = polytope->polytope_boundary();
876 for (unsigned int f = 0; f < boundary.size(); ++f)
877 {
878 myfile << boundary[f]->vertex(0)[0];
879 myfile << ",";
880 myfile << boundary[f]->vertex(0)[1];
881 myfile << ",";
882 myfile << boundary[f]->vertex(1)[0];
883 myfile << ",";
884 myfile << boundary[f]->vertex(1)[1];
885 myfile << "\n";
886 }
887 }
888
889
890 myfile.close();
891 }
892
893
894 template <typename T>
895 inline constexpr T
896 constexpr_pow(T num, unsigned int pow)
897 {
898 return (pow >= sizeof(unsigned int) * 8) ? 0 :
899 pow == 0 ? 1 :
900 num * constexpr_pow(num, pow - 1);
901 }
902
903
904
905 void
906 write_to_matrix_market_format(const std::string &filename,
907 const std::string &matrix_name,
908 const TrilinosWrappers::SparseMatrix &matrix)
909 {
910#ifdef DEAL_II_WITH_TRILINOS
911 const Epetra_CrsMatrix &trilinos_matrix = matrix.trilinos_matrix();
912
913 const int ierr =
914 EpetraExt::RowMatrixToMatrixMarketFile(filename.c_str(),
915 trilinos_matrix,
916 matrix_name.c_str(),
917 0 /*description field empty*/,
918 true /*write header*/);
919 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
920#else
921 (void)filename;
922 (void)matrix_name;
923 (void)matrix;
924#endif
925 }
926
927
928
929 namespace internal
930 {
936 template <int dim, int spacedim, typename VectorType>
937 void
938 interpolate_to_fine_grid(
939 const AgglomerationHandler<dim, spacedim> &agglomeration_handler,
940 VectorType &dst,
941 const VectorType &src)
942 {
943 Assert((dim == spacedim), ExcNotImplemented());
944 Assert(
945 dst.size() == 0,
946 ExcMessage(
947 "The destination vector must the empt upon calling this function."));
948
949 using NumberType = typename VectorType::value_type;
950 constexpr bool is_trilinos_vector =
951 std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>;
952 using MatrixType = std::conditional_t<is_trilinos_vector,
953 TrilinosWrappers::SparseMatrix,
954 SparseMatrix<NumberType>>;
955
956 MatrixType interpolation_matrix;
957
958 [[maybe_unused]]
959 typename std::conditional_t<!is_trilinos_vector, SparsityPattern, void *>
960 sp;
961
962 // Get some info from the handler
963 const DoFHandler<dim, spacedim> &agglo_dh =
964 agglomeration_handler.agglo_dh;
965
966 DoFHandler<dim, spacedim> *output_dh =
967 const_cast<DoFHandler<dim, spacedim> *>(
968 &agglomeration_handler.output_dh);
969 const FiniteElement<dim, spacedim> &fe = agglomeration_handler.get_fe();
970 const Mapping<dim> &mapping = agglomeration_handler.get_mapping();
971 const Triangulation<dim, spacedim> &tria =
972 agglomeration_handler.get_triangulation();
973 const auto &bboxes = agglomeration_handler.get_local_bboxes();
974
975 std::unique_ptr<FiniteElement<dim>> output_fe;
976 if (tria.all_reference_cells_are_hyper_cube())
977 output_fe = std::make_unique<FE_DGQ<dim>>(fe.degree);
978 else if (tria.all_reference_cells_are_simplex())
979 output_fe = std::make_unique<FE_SimplexDGP<dim>>(fe.degree);
980 else
981 AssertThrow(false, ExcNotImplemented());
982
983 // Setup an auxiliary DoFHandler for output purposes
984 output_dh->reinit(tria);
985 output_dh->distribute_dofs(*output_fe);
986
987 const IndexSet &locally_owned_dofs = output_dh->locally_owned_dofs();
988 const IndexSet locally_relevant_dofs =
989 DoFTools::extract_locally_relevant_dofs(*output_dh);
990
991 const IndexSet &locally_owned_dofs_agglo = agglo_dh.locally_owned_dofs();
992
993
994 DynamicSparsityPattern dsp(output_dh->n_dofs(),
995 agglo_dh.n_dofs(),
996 locally_relevant_dofs);
997
998 std::vector<types::global_dof_index> agglo_dof_indices(fe.dofs_per_cell);
999 std::vector<types::global_dof_index> standard_dof_indices(
1000 fe.dofs_per_cell);
1001 std::vector<types::global_dof_index> output_dof_indices(
1002 output_fe->dofs_per_cell);
1003
1004 Quadrature<dim> quad(output_fe->get_unit_support_points());
1005 FEValues<dim, spacedim> output_fe_values(mapping,
1006 *output_fe,
1007 quad,
1008 update_quadrature_points);
1009
1010 for (const auto &cell : agglo_dh.active_cell_iterators())
1011 if (cell->is_locally_owned())
1012 {
1013 if (agglomeration_handler.is_master_cell(cell))
1014 {
1015 auto slaves = agglomeration_handler.get_slaves_of_idx(
1016 cell->active_cell_index());
1017 slaves.emplace_back(cell);
1018
1019 cell->get_dof_indices(agglo_dof_indices);
1020
1021 for (const auto &slave : slaves)
1022 {
1023 // addd master-slave relationship
1024 const auto slave_output =
1025 slave->as_dof_handler_iterator(*output_dh);
1026 slave_output->get_dof_indices(output_dof_indices);
1027 for (const auto row : output_dof_indices)
1028 dsp.add_entries(row,
1029 agglo_dof_indices.begin(),
1030 agglo_dof_indices.end());
1031 }
1032 }
1033 }
1034
1035
1036 const auto assemble_interpolation_matrix = [&]() {
1037 FullMatrix<NumberType> local_matrix(fe.dofs_per_cell, fe.dofs_per_cell);
1038 std::vector<Point<dim>> reference_q_points(fe.dofs_per_cell);
1039
1040 // Dummy AffineConstraints, only needed for loc2glb
1041 AffineConstraints<NumberType> c;
1042 c.close();
1043
1044 for (const auto &cell : agglo_dh.active_cell_iterators())
1045 if (cell->is_locally_owned())
1046 {
1047 if (agglomeration_handler.is_master_cell(cell))
1048 {
1049 auto slaves = agglomeration_handler.get_slaves_of_idx(
1050 cell->active_cell_index());
1051 slaves.emplace_back(cell);
1052
1053 cell->get_dof_indices(agglo_dof_indices);
1054
1055 const types::global_cell_index polytope_index =
1056 agglomeration_handler.cell_to_polytope_index(cell);
1057
1058 // Get the box of this agglomerate.
1059 const BoundingBox<dim> &box = bboxes[polytope_index];
1060
1061 for (const auto &slave : slaves)
1062 {
1063 // add master-slave relationship
1064 const auto slave_output =
1065 slave->as_dof_handler_iterator(*output_dh);
1066
1067 slave_output->get_dof_indices(output_dof_indices);
1068 output_fe_values.reinit(slave_output);
1069
1070 local_matrix = 0.;
1071
1072 const auto &q_points =
1073 output_fe_values.get_quadrature_points();
1074 for (const auto i : output_fe_values.dof_indices())
1075 {
1076 const auto &p = box.real_to_unit(q_points[i]);
1077 for (const auto j : output_fe_values.dof_indices())
1078 {
1079 local_matrix(i, j) = fe.shape_value(j, p);
1080 }
1081 }
1082 c.distribute_local_to_global(local_matrix,
1083 output_dof_indices,
1084 agglo_dof_indices,
1085 interpolation_matrix);
1086 }
1087 }
1088 }
1089 };
1090
1091
1092 if constexpr (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix>)
1093 {
1094 const MPI_Comm &communicator = tria.get_communicator();
1095 SparsityTools::distribute_sparsity_pattern(dsp,
1096 locally_owned_dofs,
1097 communicator,
1098 locally_relevant_dofs);
1099
1100 interpolation_matrix.reinit(locally_owned_dofs,
1101 locally_owned_dofs_agglo,
1102 dsp,
1103 communicator);
1104 dst.reinit(locally_owned_dofs);
1105 assemble_interpolation_matrix();
1106 }
1107 else if constexpr (std::is_same_v<MatrixType, SparseMatrix<NumberType>>)
1108 {
1109 sp.copy_from(dsp);
1110 interpolation_matrix.reinit(sp);
1111 dst.reinit(output_dh->n_dofs());
1112 assemble_interpolation_matrix();
1113 }
1114 else
1115 {
1116 // PETSc, LA::d::v options not implemented.
1117 (void)agglomeration_handler;
1118 (void)dst;
1119 (void)src;
1120 AssertThrow(false, ExcNotImplemented());
1121 }
1122
1123 // If tria is distributed
1124 if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1125 &tria) != nullptr)
1126 interpolation_matrix.compress(VectorOperation::add);
1127
1128 // Finally, perform the interpolation.
1129 interpolation_matrix.vmult(dst, src);
1130 }
1131 } // namespace internal
1132
1133
1134
1145 template <int dim, int spacedim, typename VectorType>
1146 void
1147 interpolate_to_fine_grid(
1148 const AgglomerationHandler<dim, spacedim> &agglomeration_handler,
1149 VectorType &dst,
1150 const VectorType &src,
1151 const bool on_the_fly = true)
1152 {
1153 Assert((dim == spacedim), ExcNotImplemented());
1154 Assert(
1155 dst.size() == 0,
1156 ExcMessage(
1157 "The destination vector must the empt upon calling this function."));
1158
1159 using NumberType = typename VectorType::value_type;
1160 static constexpr bool is_trilinos_vector =
1161 std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>;
1162
1163 static constexpr bool is_supported_vector =
1164 std::is_same_v<VectorType, Vector<NumberType>> || is_trilinos_vector;
1165 static_assert(is_supported_vector);
1166
1167 // First, check for an easy return
1168 if (on_the_fly == false)
1169 {
1170 return internal::interpolate_to_fine_grid(agglomeration_handler,
1171 dst,
1172 src);
1173 }
1174 else
1175 {
1176 // otherwise, do not create any matrix
1177 const Triangulation<dim, spacedim> &tria =
1178 agglomeration_handler.get_triangulation();
1179 const Mapping<dim> &mapping = agglomeration_handler.get_mapping();
1180 const FiniteElement<dim, spacedim> &original_fe =
1181 agglomeration_handler.get_fe();
1182
1183 // We use DGQ (on tensor-product meshes) or DGP (on simplex meshes)
1184 // nodal elements of the same degree as the ones in the agglomeration
1185 // handler to interpolate the solution onto the finer grid.
1186 std::unique_ptr<FiniteElement<dim>> output_fe;
1187 if (tria.all_reference_cells_are_hyper_cube())
1188 output_fe = std::make_unique<FE_DGQ<dim>>(original_fe.degree);
1189 else if (tria.all_reference_cells_are_simplex())
1190 output_fe = std::make_unique<FE_SimplexDGP<dim>>(original_fe.degree);
1191 else
1192 AssertThrow(false, ExcNotImplemented());
1193
1194 DoFHandler<dim> &output_dh =
1195 const_cast<DoFHandler<dim> &>(agglomeration_handler.output_dh);
1196 output_dh.reinit(tria);
1197 output_dh.distribute_dofs(*output_fe);
1198
1199 if constexpr (std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>)
1200 {
1201 const IndexSet &locally_owned_dofs = output_dh.locally_owned_dofs();
1202 dst.reinit(locally_owned_dofs);
1203 }
1204 else if constexpr (std::is_same_v<VectorType, Vector<NumberType>>)
1205 {
1206 dst.reinit(output_dh.n_dofs());
1207 }
1208 else
1209 {
1210 // PETSc, LA::d::v options not implemented.
1211 (void)agglomeration_handler;
1212 (void)dst;
1213 (void)src;
1214 AssertThrow(false, ExcNotImplemented());
1215 }
1216
1217
1218
1219 const unsigned int dofs_per_cell =
1220 agglomeration_handler.n_dofs_per_cell();
1221 const unsigned int output_dofs_per_cell = output_fe->n_dofs_per_cell();
1222 Quadrature<dim> quad(output_fe->get_unit_support_points());
1223 FEValues<dim> output_fe_values(mapping,
1224 *output_fe,
1225 quad,
1226 update_quadrature_points);
1227
1228 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1229 std::vector<types::global_dof_index> local_dof_indices_output(
1230 output_dofs_per_cell);
1231
1232 const auto &bboxes = agglomeration_handler.get_local_bboxes();
1233 for (const auto &polytope : agglomeration_handler.polytope_iterators())
1234 {
1235 if (polytope->is_locally_owned())
1236 {
1237 polytope->get_dof_indices(local_dof_indices);
1238 const BoundingBox<dim> &box = bboxes[polytope->index()];
1239
1240 const auto &deal_cells =
1241 polytope->get_agglomerate(); // fine deal.II cells
1242 for (const auto &cell : deal_cells)
1243 {
1244 const auto slave_output = cell->as_dof_handler_iterator(
1245 agglomeration_handler.output_dh);
1246 slave_output->get_dof_indices(local_dof_indices_output);
1247 output_fe_values.reinit(slave_output);
1248
1249 const auto &qpoints =
1250 output_fe_values.get_quadrature_points();
1251
1252 for (unsigned int j = 0; j < output_dofs_per_cell; ++j)
1253 {
1254 const auto &ref_qpoint = box.real_to_unit(qpoints[j]);
1255 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1256 dst(local_dof_indices_output[j]) +=
1257 src(local_dof_indices[i]) *
1258 original_fe.shape_value(i, ref_qpoint);
1259 }
1260 }
1261 }
1262 }
1263 }
1264 }
1265
1266
1267
1275 template <int dim, int spacedim, typename MatrixType>
1276 void
1277 fill_interpolation_matrix(
1278 const AgglomerationHandler<dim, spacedim> &agglomeration_handler,
1279 MatrixType &interpolation_matrix)
1280 {
1281 Assert((dim == spacedim), ExcNotImplemented());
1282
1283 using NumberType = typename MatrixType::value_type;
1284 constexpr bool is_trilinos_matrix =
1285 std::is_same_v<MatrixType, TrilinosWrappers::MPI::Vector>;
1286
1287 [[maybe_unused]]
1288 typename std::conditional_t<!is_trilinos_matrix, SparsityPattern, void *>
1289 sp;
1290
1291 // Get some info from the handler
1292 const DoFHandler<dim, spacedim> &agglo_dh = agglomeration_handler.agglo_dh;
1293
1294 DoFHandler<dim, spacedim> *output_dh =
1295 const_cast<DoFHandler<dim, spacedim> *>(&agglomeration_handler.output_dh);
1296 const Mapping<dim, spacedim> &mapping = agglomeration_handler.get_mapping();
1297 const FiniteElement<dim, spacedim> &fe = agglomeration_handler.get_fe();
1298 const Triangulation<dim, spacedim> &tria =
1299 agglomeration_handler.get_triangulation();
1300 const auto &bboxes = agglomeration_handler.get_local_bboxes();
1301
1302 // Setup an auxiliary DoFHandler for output purposes
1303 output_dh->reinit(tria);
1304 output_dh->distribute_dofs(fe);
1305
1306 const IndexSet &locally_owned_dofs = output_dh->locally_owned_dofs();
1307 const IndexSet locally_relevant_dofs =
1308 DoFTools::extract_locally_relevant_dofs(*output_dh);
1309
1310 const IndexSet &locally_owned_dofs_agglo = agglo_dh.locally_owned_dofs();
1311
1312
1313 DynamicSparsityPattern dsp(output_dh->n_dofs(),
1314 agglo_dh.n_dofs(),
1315 locally_relevant_dofs);
1316
1317 std::vector<types::global_dof_index> agglo_dof_indices(fe.dofs_per_cell);
1318 std::vector<types::global_dof_index> standard_dof_indices(fe.dofs_per_cell);
1319 std::vector<types::global_dof_index> output_dof_indices(fe.dofs_per_cell);
1320
1321 Quadrature<dim> quad(fe.get_unit_support_points());
1322 FEValues<dim, spacedim> output_fe_values(mapping,
1323 fe,
1324 quad,
1325 update_quadrature_points);
1326
1327 for (const auto &cell : agglo_dh.active_cell_iterators())
1328 if (cell->is_locally_owned())
1329 {
1330 if (agglomeration_handler.is_master_cell(cell))
1331 {
1332 auto slaves = agglomeration_handler.get_slaves_of_idx(
1333 cell->active_cell_index());
1334 slaves.emplace_back(cell);
1335
1336 cell->get_dof_indices(agglo_dof_indices);
1337
1338 for (const auto &slave : slaves)
1339 {
1340 // addd master-slave relationship
1341 const auto slave_output =
1342 slave->as_dof_handler_iterator(*output_dh);
1343 slave_output->get_dof_indices(output_dof_indices);
1344 for (const auto row : output_dof_indices)
1345 dsp.add_entries(row,
1346 agglo_dof_indices.begin(),
1347 agglo_dof_indices.end());
1348 }
1349 }
1350 }
1351
1352
1353 const auto assemble_interpolation_matrix = [&]() {
1354 FullMatrix<NumberType> local_matrix(fe.dofs_per_cell, fe.dofs_per_cell);
1355 std::vector<Point<dim>> reference_q_points(fe.dofs_per_cell);
1356
1357 // Dummy AffineConstraints, only needed for loc2glb
1358 AffineConstraints<NumberType> c;
1359 c.close();
1360
1361 for (const auto &cell : agglo_dh.active_cell_iterators())
1362 if (cell->is_locally_owned())
1363 {
1364 if (agglomeration_handler.is_master_cell(cell))
1365 {
1366 auto slaves = agglomeration_handler.get_slaves_of_idx(
1367 cell->active_cell_index());
1368 slaves.emplace_back(cell);
1369
1370 cell->get_dof_indices(agglo_dof_indices);
1371
1372 const types::global_cell_index polytope_index =
1373 agglomeration_handler.cell_to_polytope_index(cell);
1374
1375 // Get the box of this agglomerate.
1376 const BoundingBox<dim> &box = bboxes[polytope_index];
1377
1378 for (const auto &slave : slaves)
1379 {
1380 // add master-slave relationship
1381 const auto slave_output =
1382 slave->as_dof_handler_iterator(*output_dh);
1383
1384 slave_output->get_dof_indices(output_dof_indices);
1385 output_fe_values.reinit(slave_output);
1386
1387 local_matrix = 0.;
1388
1389 const auto &q_points =
1390 output_fe_values.get_quadrature_points();
1391 for (const auto i : output_fe_values.dof_indices())
1392 {
1393 const auto &p = box.real_to_unit(q_points[i]);
1394 for (const auto j : output_fe_values.dof_indices())
1395 {
1396 local_matrix(i, j) = fe.shape_value(j, p);
1397 }
1398 }
1399 c.distribute_local_to_global(local_matrix,
1400 output_dof_indices,
1401 agglo_dof_indices,
1402 interpolation_matrix);
1403 }
1404 }
1405 }
1406 };
1407
1408
1409 if constexpr (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix>)
1410 {
1411 const MPI_Comm &communicator = tria.get_communicator();
1412 SparsityTools::distribute_sparsity_pattern(dsp,
1413 locally_owned_dofs,
1414 communicator,
1415 locally_relevant_dofs);
1416
1417 interpolation_matrix.reinit(locally_owned_dofs,
1418 locally_owned_dofs_agglo,
1419 dsp,
1420 communicator);
1421 assemble_interpolation_matrix();
1422 }
1423 else if constexpr (std::is_same_v<MatrixType, SparseMatrix<NumberType>>)
1424 {
1425 sp.copy_from(dsp);
1426 interpolation_matrix.reinit(sp);
1427 assemble_interpolation_matrix();
1428 }
1429 else
1430 {
1431 // PETSc, LA::d::v options not implemented.
1432 (void)agglomeration_handler;
1433 AssertThrow(false, ExcNotImplemented());
1434 }
1435
1436 // If tria is distributed
1437 if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1438 &tria) != nullptr)
1439 interpolation_matrix.compress(VectorOperation::add);
1440 }
1441
1442
1443
1453 template <int dim, typename Number, typename VectorType>
1454 void
1455 compute_global_error(const AgglomerationHandler<dim> &agglomeration_handler,
1456 const VectorType &solution,
1457 const Function<dim, Number> &exact_solution,
1458 const std::vector<VectorTools::NormType> &norms,
1459 std::vector<double> &global_errors)
1460 {
1461 Assert(solution.size() > 0,
1462 ExcNotImplemented(
1463 "Solution vector must be non-empty upon calling this function."));
1464 Assert(std::any_of(norms.cbegin(),
1465 norms.cend(),
1466 [](VectorTools::NormType norm_type) {
1467 return (norm_type ==
1468 VectorTools::NormType::H1_seminorm ||
1469 norm_type == VectorTools::NormType::L2_norm);
1470 }),
1471 ExcMessage("Norm type not supported"));
1472 global_errors.resize(norms.size());
1473 std::fill(global_errors.begin(), global_errors.end(), 0.);
1474
1475 // Vector storing errors local to the current processor.
1476 std::vector<double> local_errors(norms.size());
1477 std::fill(local_errors.begin(), local_errors.end(), 0.);
1478
1479 // Get some info from the handler
1480 const unsigned int dofs_per_cell = agglomeration_handler.n_dofs_per_cell();
1481
1482 const bool compute_semi_H1 =
1483 std::any_of(norms.cbegin(),
1484 norms.cend(),
1485 [](VectorTools::NormType norm_type) {
1486 return norm_type == VectorTools::NormType::H1_seminorm;
1487 });
1488
1489 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1490 for (const auto &polytope : agglomeration_handler.polytope_iterators())
1491 {
1492 if (polytope->is_locally_owned())
1493 {
1494 const auto &agglo_values = agglomeration_handler.reinit(polytope);
1495 polytope->get_dof_indices(local_dof_indices);
1496
1497 const auto &q_points = agglo_values.get_quadrature_points();
1498 const unsigned int n_qpoints = q_points.size();
1499 std::vector<double> analyical_sol_at_qpoints(n_qpoints);
1500 exact_solution.value_list(q_points, analyical_sol_at_qpoints);
1501 std::vector<Tensor<1, dim>> grad_analyical_sol_at_qpoints(
1502 n_qpoints);
1503
1504 if (compute_semi_H1)
1505 exact_solution.gradient_list(q_points,
1506 grad_analyical_sol_at_qpoints);
1507
1508 for (unsigned int q_index : agglo_values.quadrature_point_indices())
1509 {
1510 double solution_at_qpoint = 0.;
1511 Tensor<1, dim> grad_solution_at_qpoint;
1512 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1513 {
1514 solution_at_qpoint += solution(local_dof_indices[i]) *
1515 agglo_values.shape_value(i, q_index);
1516
1517 if (compute_semi_H1)
1518 grad_solution_at_qpoint +=
1519 solution(local_dof_indices[i]) *
1520 agglo_values.shape_grad(i, q_index);
1521 }
1522 // L2
1523 local_errors[0] += std::pow((analyical_sol_at_qpoints[q_index] -
1524 solution_at_qpoint),
1525 2) *
1526 agglo_values.JxW(q_index);
1527
1528 // H1 seminorm
1529 if (compute_semi_H1)
1530 for (unsigned int d = 0; d < dim; ++d)
1531 local_errors[1] +=
1532 std::pow((grad_analyical_sol_at_qpoints[q_index][d] -
1533 grad_solution_at_qpoint[d]),
1534 2) *
1535 agglo_values.JxW(q_index);
1536 }
1537 }
1538 }
1539
1540 // Perform reduction and take sqrt of each error
1541 global_errors[0] = Utilities::MPI::reduce<double>(
1542 local_errors[0],
1543 agglomeration_handler.get_triangulation().get_mpi_communicator(),
1544 [](const double a, const double b) { return a + b; });
1545
1546 global_errors[0] = std::sqrt(global_errors[0]);
1547
1548 if (compute_semi_H1)
1549 {
1550 global_errors[1] = Utilities::MPI::reduce<double>(
1551 local_errors[1],
1552 agglomeration_handler.get_triangulation().get_mpi_communicator(),
1553 [](const double a, const double b) { return a + b; });
1554 global_errors[1] = std::sqrt(global_errors[1]);
1555 }
1556 }
1557
1558
1559
1567 template <int dim>
1568 unsigned int
1569 construct_agglomerated_levels(
1570 const Triangulation<dim> &tria,
1571 std::vector<std::unique_ptr<AgglomerationHandler<dim>>>
1572 &agglomeration_handlers,
1573 const FE_DGQ<dim> &fe_dg,
1574 const Mapping<dim> &mapping,
1575 const unsigned int starting_tree_level)
1576 {
1577 const auto parallel_tria =
1578 dynamic_cast<const parallel::TriangulationBase<dim> *>(&tria);
1579
1580 GridTools::Cache<dim> cached_tria(tria);
1581 Assert(parallel_tria->n_active_cells() > 0, ExcInternalError());
1582
1583 const MPI_Comm comm = parallel_tria->get_communicator();
1584 ConditionalOStream pcout(std::cout,
1585 (Utilities::MPI::this_mpi_process(comm) == 0));
1586
1587 // Start building R-tree
1588 namespace bgi = boost::geometry::index;
1589 static constexpr unsigned int max_elem_per_node =
1590 constexpr_pow(2, dim); // 2^dim
1591 std::vector<std::pair<BoundingBox<dim>,
1592 typename Triangulation<dim>::active_cell_iterator>>
1593 boxes(parallel_tria->n_locally_owned_active_cells());
1594 unsigned int i = 0;
1595 for (const auto &cell : parallel_tria->active_cell_iterators())
1596 if (cell->is_locally_owned())
1597 boxes[i++] = std::make_pair(mapping.get_bounding_box(cell), cell);
1598
1599 auto tree = pack_rtree<bgi::rstar<max_elem_per_node>>(boxes);
1600 Assert(n_levels(tree) >= 2, ExcMessage("At least two levels are needed."));
1601 pcout << "Total number of available levels: " << n_levels(tree)
1602 << std::endl;
1603
1604 pcout << "Starting level: " << starting_tree_level << std::endl;
1605 const unsigned int total_tree_levels =
1606 n_levels(tree) - starting_tree_level + 1;
1607
1608
1609 // Resize the agglomeration handlers to the right size
1610
1611 agglomeration_handlers.resize(total_tree_levels);
1612 // Loop through the available levels and set AgglomerationHandlers up.
1613 for (unsigned int extraction_level = starting_tree_level;
1614 extraction_level <= n_levels(tree);
1615 ++extraction_level)
1616 {
1617 agglomeration_handlers[extraction_level - starting_tree_level] =
1618 std::make_unique<AgglomerationHandler<dim>>(cached_tria);
1619 CellsAgglomerator<dim, decltype(tree)> agglomerator{tree,
1620 extraction_level};
1621 const auto agglomerates = agglomerator.extract_agglomerates();
1622 agglomeration_handlers[extraction_level - starting_tree_level]
1623 ->connect_hierarchy(agglomerator);
1624
1625 // Flag elements for agglomeration
1626 unsigned int agglo_index = 0;
1627 for (unsigned int i = 0; i < agglomerates.size(); ++i)
1628 {
1629 const auto &agglo = agglomerates[i]; // i-th agglomerate
1630 for (const auto &el : agglo)
1631 {
1632 el->set_material_id(agglo_index);
1633 }
1634 ++agglo_index;
1635 }
1636
1637 const unsigned int n_local_agglomerates = agglo_index;
1638 unsigned int total_agglomerates =
1639 Utilities::MPI::sum(n_local_agglomerates, comm);
1640 pcout << "Total agglomerates per (tree) level: " << extraction_level
1641 << ": " << total_agglomerates << std::endl;
1642
1643
1644 // Now, perform agglomeration within each locally owned partition
1645 std::vector<
1646 std::vector<typename Triangulation<dim>::active_cell_iterator>>
1647 cells_per_subdomain(n_local_agglomerates);
1648 for (const auto &cell : parallel_tria->active_cell_iterators())
1649 if (cell->is_locally_owned())
1650 cells_per_subdomain[cell->material_id()].push_back(cell);
1651
1652 // For every subdomain, agglomerate elements together
1653 for (std::size_t i = 0; i < cells_per_subdomain.size(); ++i)
1654 agglomeration_handlers[extraction_level - starting_tree_level]
1655 ->define_agglomerate(cells_per_subdomain[i]);
1656
1657 agglomeration_handlers[extraction_level - starting_tree_level]
1658 ->initialize_fe_values(QGauss<dim>(fe_dg.degree + 1),
1659 update_values | update_gradients |
1660 update_JxW_values | update_quadrature_points,
1661 QGauss<dim - 1>(fe_dg.degree + 1),
1662 update_JxW_values);
1663 agglomeration_handlers[extraction_level - starting_tree_level]
1664 ->distribute_agglomerated_dofs(fe_dg);
1665 }
1666
1667 return total_tree_levels;
1668 }
1669
1670
1671
1676 template <int dim>
1677 void
1678 assemble_local_jumps_and_averages(FullMatrix<double> &M11,
1679 FullMatrix<double> &M12,
1680 FullMatrix<double> &M21,
1681 FullMatrix<double> &M22,
1682 const FEValuesBase<dim> &fe_faces0,
1683 const FEValuesBase<dim> &fe_faces1,
1684 const double penalty_constant,
1685 const double h_f)
1686 {
1687 const std::vector<Tensor<1, dim>> &normals = fe_faces0.get_normal_vectors();
1688 const unsigned int dofs_per_cell =
1689 M11.m(); // size of local matrices equals the #DoFs
1690 for (unsigned int q_index : fe_faces0.quadrature_point_indices())
1691 {
1692 const Tensor<1, dim> &normal = normals[q_index];
1693 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1694 {
1695 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1696 {
1697 M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normal *
1698 fe_faces0.shape_value(j, q_index) -
1699 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
1700 fe_faces0.shape_value(i, q_index) +
1701 (penalty_constant / h_f) *
1702 fe_faces0.shape_value(i, q_index) *
1703 fe_faces0.shape_value(j, q_index)) *
1704 fe_faces0.JxW(q_index);
1705 M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normal *
1706 fe_faces1.shape_value(j, q_index) -
1707 0.5 * fe_faces1.shape_grad(j, q_index) * normal *
1708 fe_faces0.shape_value(i, q_index) -
1709 (penalty_constant / h_f) *
1710 fe_faces0.shape_value(i, q_index) *
1711 fe_faces1.shape_value(j, q_index)) *
1712 fe_faces1.JxW(q_index);
1713 M21(i, j) += (-0.5 * fe_faces1.shape_grad(i, q_index) * normal *
1714 fe_faces0.shape_value(j, q_index) +
1715 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
1716 fe_faces1.shape_value(i, q_index) -
1717 (penalty_constant / h_f) *
1718 fe_faces1.shape_value(i, q_index) *
1719 fe_faces0.shape_value(j, q_index)) *
1720 fe_faces1.JxW(q_index);
1721 M22(i, j) += (0.5 * fe_faces1.shape_grad(i, q_index) * normal *
1722 fe_faces1.shape_value(j, q_index) +
1723 0.5 * fe_faces1.shape_grad(j, q_index) * normal *
1724 fe_faces1.shape_value(i, q_index) +
1725 (penalty_constant / h_f) *
1726 fe_faces1.shape_value(i, q_index) *
1727 fe_faces1.shape_value(j, q_index)) *
1728 fe_faces1.JxW(q_index);
1729 }
1730 }
1731 }
1732 }
1736 template <int dim>
1737 void
1738 assemble_local_jumps_and_averages_ghost(
1739 FullMatrix<double> &M11,
1740 FullMatrix<double> &M12,
1741 FullMatrix<double> &M21,
1742 FullMatrix<double> &M22,
1743 const FEValuesBase<dim> &fe_faces0,
1744 const std::vector<std::vector<double>> &recv_values,
1745 const std::vector<std::vector<Tensor<1, dim>>> &recv_gradients,
1746 const std::vector<double> &recv_jxws,
1747 const double penalty_constant,
1748 const double h_f)
1749 {
1750 Assert(
1751 (recv_values.size() > 0 && recv_gradients.size() && recv_jxws.size()),
1752 ExcMessage(
1753 "Not possible to assemble jumps and averages at a ghosted interface."));
1754 const unsigned int dofs_per_cell = M11.m();
1755 const std::vector<Tensor<1, dim>> &normals = fe_faces0.get_normal_vectors();
1756 for (unsigned int q_index : fe_faces0.quadrature_point_indices())
1757 {
1758 const Tensor<1, dim> &normal = normals[q_index];
1759 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1760 {
1761 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1762 {
1763 M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normal *
1764 fe_faces0.shape_value(j, q_index) -
1765 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
1766 fe_faces0.shape_value(i, q_index) +
1767 (penalty_constant / h_f) *
1768 fe_faces0.shape_value(i, q_index) *
1769 fe_faces0.shape_value(j, q_index)) *
1770 fe_faces0.JxW(q_index);
1771 M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normal *
1772 recv_values[j][q_index] -
1773 0.5 * recv_gradients[j][q_index] * normal *
1774 fe_faces0.shape_value(i, q_index) -
1775 (penalty_constant / h_f) *
1776 fe_faces0.shape_value(i, q_index) *
1777 recv_values[j][q_index]) *
1778 recv_jxws[q_index];
1779 M21(i, j) +=
1780 (-0.5 * recv_gradients[i][q_index] * normal *
1781 fe_faces0.shape_value(j, q_index) +
1782 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
1783 recv_values[i][q_index] -
1784 (penalty_constant / h_f) * recv_values[i][q_index] *
1785 fe_faces0.shape_value(j, q_index)) *
1786 recv_jxws[q_index];
1787 M22(i, j) +=
1788 (0.5 * recv_gradients[i][q_index] * normal *
1789 recv_values[j][q_index] +
1790 0.5 * recv_gradients[j][q_index] * normal *
1791 recv_values[i][q_index] +
1792 (penalty_constant / h_f) * recv_values[i][q_index] *
1793 recv_values[j][q_index]) *
1794 recv_jxws[q_index];
1795 }
1796 }
1797 }
1798 }
1799
1800
1806 template <int dim, typename MatrixType>
1807 void
1808 assemble_dg_matrix(MatrixType &system_matrix,
1809 const FiniteElement<dim> &fe_dg,
1810 const AgglomerationHandler<dim> &ah)
1811 {
1812 static_assert(
1813 (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix> ||
1814 std::is_same_v<MatrixType,
1815 SparseMatrix<typename MatrixType::value_type>>));
1816
1817 Assert((dynamic_cast<const FE_DGQ<dim> *>(&fe_dg) ||
1818 dynamic_cast<const FE_DGP<dim> *>(&fe_dg) ||
1819 dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_dg)),
1820 ExcMessage("FE type not supported."));
1821
1822 AffineConstraints constraints;
1823 constraints.close();
1824 const double penalty_constant =
1825 10 * (fe_dg.degree + dim) * (fe_dg.degree + 1);
1826 TrilinosWrappers::SparsityPattern dsp;
1827 const_cast<AgglomerationHandler<dim> &>(ah)
1828 .create_agglomeration_sparsity_pattern(dsp);
1829 system_matrix.reinit(dsp);
1830 const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell();
1831 FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1832 FullMatrix<double> M11(dofs_per_cell, dofs_per_cell);
1833 FullMatrix<double> M12(dofs_per_cell, dofs_per_cell);
1834 FullMatrix<double> M21(dofs_per_cell, dofs_per_cell);
1835 FullMatrix<double> M22(dofs_per_cell, dofs_per_cell);
1836 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1837 std::vector<types::global_dof_index> local_dof_indices_neighbor(
1838 dofs_per_cell);
1839
1840 for (const auto &polytope : ah.polytope_iterators())
1841 {
1842 if (polytope->is_locally_owned())
1843 {
1844 cell_matrix = 0.;
1845 const auto &agglo_values = ah.reinit(polytope);
1846 for (unsigned int q_index : agglo_values.quadrature_point_indices())
1847 {
1848 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1849 {
1850 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1851 {
1852 cell_matrix(i, j) +=
1853 agglo_values.shape_grad(i, q_index) *
1854 agglo_values.shape_grad(j, q_index) *
1855 agglo_values.JxW(q_index);
1856 }
1857 }
1858 }
1859 // get volumetric DoFs
1860 polytope->get_dof_indices(local_dof_indices);
1861 // Assemble face terms
1862 unsigned int n_faces = polytope->n_faces();
1863 const double h_f = polytope->diameter();
1864 for (unsigned int f = 0; f < n_faces; ++f)
1865 {
1866 if (polytope->at_boundary(f))
1867 {
1868 // Get normal vectors seen from each agglomeration.
1869 const auto &fe_face = ah.reinit(polytope, f);
1870 const auto &normals = fe_face.get_normal_vectors();
1871 for (unsigned int q_index :
1872 fe_face.quadrature_point_indices())
1873 {
1874 const Tensor<1, dim> &normal = normals[q_index];
1875 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1876 {
1877 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1878 {
1879 cell_matrix(i, j) +=
1880 (-fe_face.shape_value(i, q_index) *
1881 fe_face.shape_grad(j, q_index) * normal -
1882 fe_face.shape_grad(i, q_index) * normal *
1883 fe_face.shape_value(j, q_index) +
1884 (penalty_constant / h_f) *
1885 fe_face.shape_value(i, q_index) *
1886 fe_face.shape_value(j, q_index)) *
1887 fe_face.JxW(q_index);
1888 }
1889 }
1890 }
1891 }
1892 else
1893 {
1894 const auto &neigh_polytope = polytope->neighbor(f);
1895 if (polytope->id() < neigh_polytope->id())
1896 {
1897 unsigned int nofn =
1898 polytope->neighbor_of_agglomerated_neighbor(f);
1899 Assert(neigh_polytope->neighbor(nofn)->id() ==
1900 polytope->id(),
1901 ExcMessage("Mismatch."));
1902 const auto &fe_faces = ah.reinit_interface(
1903 polytope, neigh_polytope, f, nofn);
1904 const auto &fe_faces0 = fe_faces.first;
1905 if (neigh_polytope->is_locally_owned())
1906 {
1907 // use both fevalues
1908 const auto &fe_faces1 = fe_faces.second;
1909 M11 = 0.;
1910 M12 = 0.;
1911 M21 = 0.;
1912 M22 = 0.;
1913 assemble_local_jumps_and_averages(M11,
1914 M12,
1915 M21,
1916 M22,
1917 fe_faces0,
1918 fe_faces1,
1919 penalty_constant,
1920 h_f);
1921 // distribute DoFs accordingly
1922 // fluxes
1923 neigh_polytope->get_dof_indices(
1924 local_dof_indices_neighbor);
1925 constraints.distribute_local_to_global(
1926 M11, local_dof_indices, system_matrix);
1927 constraints.distribute_local_to_global(
1928 M12,
1929 local_dof_indices,
1930 local_dof_indices_neighbor,
1931 system_matrix);
1932 constraints.distribute_local_to_global(
1933 M21,
1934 local_dof_indices_neighbor,
1935 local_dof_indices,
1936 system_matrix);
1937 constraints.distribute_local_to_global(
1938 M22, local_dof_indices_neighbor, system_matrix);
1939 }
1940 else
1941 {
1942 // neigh polytope is ghosted, so retrieve necessary
1943 // metadata.
1944 types::subdomain_id neigh_rank =
1945 neigh_polytope->subdomain_id();
1946 const auto &recv_jxws =
1947 ah.recv_jxws.at(neigh_rank)
1948 .at({neigh_polytope->id(), nofn});
1949 const auto &recv_values =
1950 ah.recv_values.at(neigh_rank)
1951 .at({neigh_polytope->id(), nofn});
1952 const auto &recv_gradients =
1953 ah.recv_gradients.at(neigh_rank)
1954 .at({neigh_polytope->id(), nofn});
1955 M11 = 0.;
1956 M12 = 0.;
1957 M21 = 0.;
1958 M22 = 0.;
1959 // there's no FEFaceValues on the other side (it's
1960 // ghosted), so we just pass the actual data we have
1961 // recevied from the neighboring ghosted polytope
1962 assemble_local_jumps_and_averages_ghost(
1963 M11,
1964 M12,
1965 M21,
1966 M22,
1967 fe_faces0,
1968 recv_values,
1969 recv_gradients,
1970 recv_jxws,
1971 penalty_constant,
1972 h_f);
1973 // distribute DoFs accordingly
1974 // fluxes
1975 neigh_polytope->get_dof_indices(
1976 local_dof_indices_neighbor);
1977 constraints.distribute_local_to_global(
1978 M11, local_dof_indices, system_matrix);
1979 constraints.distribute_local_to_global(
1980 M12,
1981 local_dof_indices,
1982 local_dof_indices_neighbor,
1983 system_matrix);
1984 constraints.distribute_local_to_global(
1985 M21,
1986 local_dof_indices_neighbor,
1987 local_dof_indices,
1988 system_matrix);
1989 constraints.distribute_local_to_global(
1990 M22, local_dof_indices_neighbor, system_matrix);
1991 } // ghosted polytope case
1992 } // only once
1993 } // internal face
1994 } // face loop
1995 constraints.distribute_local_to_global(cell_matrix,
1996 local_dof_indices,
1997 system_matrix);
1998 } // locally owned polytopes
1999 }
2000 system_matrix.compress(VectorOperation::add);
2001 }
2002
2003
2004
2010 template <int dim, typename MatrixType, typename VectorType>
2011 void
2012 assemble_dg_matrix_on_standard_mesh(MatrixType &system_matrix,
2013 VectorType &system_rhs,
2014 const Mapping<dim> &mapping,
2015 const FiniteElement<dim> &fe_dg,
2016 const DoFHandler<dim> &dof_handler)
2017 {
2018 static_assert(
2019 (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix> ||
2020 std::is_same_v<MatrixType,
2021 SparseMatrix<typename MatrixType::value_type>>));
2022
2023 Assert((dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_dg) != nullptr),
2024 ExcNotImplemented(
2025 "Implemented only for simplex meshes for the time being."));
2026
2027 Assert(dof_handler.get_triangulation().all_reference_cells_are_simplex(),
2028 ExcNotImplemented());
2029
2030 const double penalty_constant = .5 * fe_dg.degree * (fe_dg.degree + 1);
2031 AffineConstraints<typename MatrixType::value_type> constraints;
2032 constraints.close();
2033
2034 const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
2035 const IndexSet locally_relevant_dofs =
2036 DoFTools::extract_locally_relevant_dofs(dof_handler);
2037
2038 DynamicSparsityPattern dsp(locally_relevant_dofs);
2039 DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
2040 SparsityTools::distribute_sparsity_pattern(
2041 dsp,
2042 dof_handler.locally_owned_dofs(),
2043 dof_handler.get_mpi_communicator(),
2044 locally_relevant_dofs);
2045
2046 system_matrix.reinit(locally_owned_dofs,
2047 locally_owned_dofs,
2048 dsp,
2049 dof_handler.get_communicator());
2050
2051 system_rhs.reinit(locally_owned_dofs, dof_handler.get_communicator());
2052
2053 const unsigned int quadrature_degree = fe_dg.degree + 1;
2054 FEFaceValues<dim> fe_faces0(mapping,
2055 fe_dg,
2056 QGaussSimplex<dim - 1>(quadrature_degree),
2057 update_values | update_JxW_values |
2058 update_gradients | update_quadrature_points |
2059 update_normal_vectors);
2060
2061
2062 FEValues<dim> fe_values(mapping,
2063 fe_dg,
2064 QGaussSimplex<dim>(quadrature_degree),
2065 update_values | update_JxW_values |
2066 update_gradients | update_quadrature_points);
2067
2068 FEFaceValues<dim> fe_faces1(mapping,
2069 fe_dg,
2070 QGaussSimplex<dim - 1>(quadrature_degree),
2071 update_values | update_JxW_values |
2072 update_gradients | update_quadrature_points |
2073 update_normal_vectors);
2074 const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell();
2075
2076 FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
2077 Vector<double> cell_rhs(dofs_per_cell);
2078
2079 FullMatrix<double> M11(dofs_per_cell, dofs_per_cell);
2080 FullMatrix<double> M12(dofs_per_cell, dofs_per_cell);
2081 FullMatrix<double> M21(dofs_per_cell, dofs_per_cell);
2082 FullMatrix<double> M22(dofs_per_cell, dofs_per_cell);
2083
2084 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2085
2086 // Loop over standard deal.II cells
2087 for (const auto &cell : dof_handler.active_cell_iterators())
2088 {
2089 if (cell->is_locally_owned())
2090 {
2091 cell_matrix = 0.;
2092 cell_rhs = 0.;
2093
2094 fe_values.reinit(cell);
2095
2096 // const auto &q_points = fe_values.get_quadrature_points();
2097 // const unsigned int n_qpoints = q_points.size();
2098
2099 for (unsigned int q_index : fe_values.quadrature_point_indices())
2100 {
2101 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2102 {
2103 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2104 {
2105 cell_matrix(i, j) += fe_values.shape_grad(i, q_index) *
2106 fe_values.shape_grad(j, q_index) *
2107 fe_values.JxW(q_index);
2108 }
2109 cell_rhs(i) +=
2110 fe_values.shape_value(i, q_index) * 1. *
2111 fe_values.JxW(q_index); // TODO: pass functional
2112 }
2113 }
2114
2115 // distribute volumetric DoFs
2116 cell->get_dof_indices(local_dof_indices);
2117 double hf = 0.;
2118 for (const auto f : cell->face_indices())
2119 {
2120 const double extent1 =
2121 cell->measure() / cell->face(f)->measure();
2122
2123 if (cell->face(f)->at_boundary())
2124 {
2125 hf = (1. / extent1 + 1. / extent1);
2126 fe_faces0.reinit(cell, f);
2127
2128 const auto &normals = fe_faces0.get_normal_vectors();
2129 for (unsigned int q_index :
2130 fe_faces0.quadrature_point_indices())
2131 {
2132 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2133 {
2134 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2135 {
2136 cell_matrix(i, j) +=
2137 (-fe_faces0.shape_value(i, q_index) *
2138 fe_faces0.shape_grad(j, q_index) *
2139 normals[q_index] -
2140 fe_faces0.shape_grad(i, q_index) *
2141 normals[q_index] *
2142 fe_faces0.shape_value(j, q_index) +
2143 (penalty_constant * hf) *
2144 fe_faces0.shape_value(i, q_index) *
2145 fe_faces0.shape_value(j, q_index)) *
2146 fe_faces0.JxW(q_index);
2147 }
2148 cell_rhs(i) +=
2149 0.; // TODO: add bdary conditions functional
2150 }
2151 }
2152 }
2153 else
2154 {
2155 const auto &neigh_cell = cell->neighbor(f);
2156 if (cell->global_active_cell_index() <
2157 neigh_cell->global_active_cell_index())
2158 {
2159 const double extent2 =
2160 neigh_cell->measure() /
2161 neigh_cell->face(cell->neighbor_of_neighbor(f))
2162 ->measure();
2163 hf = (1. / extent1 + 1. / extent2);
2164 fe_faces0.reinit(cell, f);
2165 fe_faces1.reinit(neigh_cell,
2166 cell->neighbor_of_neighbor(f));
2167
2168 std::vector<types::global_dof_index>
2169 local_dof_indices_neighbor(dofs_per_cell);
2170
2171 M11 = 0.;
2172 M12 = 0.;
2173 M21 = 0.;
2174 M22 = 0.;
2175
2176 const auto &normals = fe_faces0.get_normal_vectors();
2177 // M11
2178 for (unsigned int q_index :
2179 fe_faces0.quadrature_point_indices())
2180 {
2181 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2182 {
2183 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2184 {
2185 M11(i, j) +=
2186 (-0.5 * fe_faces0.shape_grad(i, q_index) *
2187 normals[q_index] *
2188 fe_faces0.shape_value(j, q_index) -
2189 0.5 * fe_faces0.shape_grad(j, q_index) *
2190 normals[q_index] *
2191 fe_faces0.shape_value(i, q_index) +
2192 (penalty_constant * hf) *
2193 fe_faces0.shape_value(i, q_index) *
2194 fe_faces0.shape_value(j, q_index)) *
2195 fe_faces0.JxW(q_index);
2196
2197 M12(i, j) +=
2198 (0.5 * fe_faces0.shape_grad(i, q_index) *
2199 normals[q_index] *
2200 fe_faces1.shape_value(j, q_index) -
2201 0.5 * fe_faces1.shape_grad(j, q_index) *
2202 normals[q_index] *
2203 fe_faces0.shape_value(i, q_index) -
2204 (penalty_constant * hf) *
2205 fe_faces0.shape_value(i, q_index) *
2206 fe_faces1.shape_value(j, q_index)) *
2207 fe_faces1.JxW(q_index);
2208
2209 // A10
2210 M21(i, j) +=
2211 (-0.5 * fe_faces1.shape_grad(i, q_index) *
2212 normals[q_index] *
2213 fe_faces0.shape_value(j, q_index) +
2214 0.5 * fe_faces0.shape_grad(j, q_index) *
2215 normals[q_index] *
2216 fe_faces1.shape_value(i, q_index) -
2217 (penalty_constant * hf) *
2218 fe_faces1.shape_value(i, q_index) *
2219 fe_faces0.shape_value(j, q_index)) *
2220 fe_faces1.JxW(q_index);
2221
2222 // A11
2223 M22(i, j) +=
2224 (0.5 * fe_faces1.shape_grad(i, q_index) *
2225 normals[q_index] *
2226 fe_faces1.shape_value(j, q_index) +
2227 0.5 * fe_faces1.shape_grad(j, q_index) *
2228 normals[q_index] *
2229 fe_faces1.shape_value(i, q_index) +
2230 (penalty_constant * hf) *
2231 fe_faces1.shape_value(i, q_index) *
2232 fe_faces1.shape_value(j, q_index)) *
2233 fe_faces1.JxW(q_index);
2234 }
2235 }
2236 }
2237
2238 // distribute DoFs accordingly
2239
2240 neigh_cell->get_dof_indices(local_dof_indices_neighbor);
2241
2242 constraints.distribute_local_to_global(
2243 M11, local_dof_indices, system_matrix);
2244 constraints.distribute_local_to_global(
2245 M12,
2246 local_dof_indices,
2247 local_dof_indices_neighbor,
2248 system_matrix);
2249 constraints.distribute_local_to_global(
2250 M21,
2251 local_dof_indices_neighbor,
2252 local_dof_indices,
2253 system_matrix);
2254 constraints.distribute_local_to_global(
2255 M22, local_dof_indices_neighbor, system_matrix);
2256
2257 } // check idx neighbors
2258 } // over faces
2259 }
2260 constraints.distribute_local_to_global(cell_matrix,
2261 cell_rhs,
2262 local_dof_indices,
2263 system_matrix,
2264 system_rhs);
2265 }
2266 }
2267 system_matrix.compress(VectorOperation::add);
2268 system_rhs.compress(VectorOperation::add);
2269 }
2270
2271
2272} // namespace dealii::PolyUtils
2273
2274#endif
bounding_box.h
point.h
AffineConstraints::distribute_local_to_global
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
BoundingBox::real_to_unit
Point< spacedim, Number > real_to_unit(const Point< spacedim, Number > &point) const
DoFHandler::distribute_dofs
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DoFHandler::locally_owned_dofs
const IndexSet & locally_owned_dofs() const
DoFHandler::reinit
void reinit(const Triangulation< dim, spacedim > &tria)
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
DoFHandler::get_mpi_communicator
MPI_Comm get_mpi_communicator() const
DoFHandler::get_communicator
MPI_Comm get_communicator() const
DynamicSparsityPattern::add_entries
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
DynamicSparsityPattern::reinit
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
DynamicSparsityPattern::add
void add(const size_type i, const size_type j)
FEFaceValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell, const unsigned int face_no)
FEValues::get_quadrature_points
const std::vector< Point< spacedim > > & get_quadrature_points() const
FEValuesBase::get_normal_vectors
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
FEValues::dof_indices
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
FEValuesBase::quadrature_point_indices
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
FEValuesBase::shape_grad
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
FEValuesBase::JxW
double JxW(const unsigned int q_point) const
FEValuesBase::shape_value
const double & shape_value(const unsigned int i, const unsigned int q_point) const
FEValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
FE_DGP
FE_DGQ
FE_SimplexDGP
FiniteElement::degree
const unsigned int degree
FiniteElement::n_dofs_per_cell
unsigned int n_dofs_per_cell() const
FiniteElement::dofs_per_cell
const unsigned int dofs_per_cell
FiniteElement::get_unit_support_points
const std::vector< Point< dim > > & get_unit_support_points() const
FiniteElement::shape_value
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
FullMatrix::m
size_type m() const
Function::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim, RangeNumberType > > &gradients, const unsigned int component=0) const
Function::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
Mapping::get_bounding_box
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
SparsityPattern::copy_from
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
Triangulation::get_communicator
MPI_Comm get_communicator() const
Triangulation::all_reference_cells_are_hyper_cube
bool all_reference_cells_are_hyper_cube() const
Triangulation::get_mpi_communicator
virtual MPI_Comm get_mpi_communicator() const
Triangulation::all_reference_cells_are_simplex
bool all_reference_cells_are_simplex() const
AgglomerationHandler::n_dofs_per_cell
unsigned int n_dofs_per_cell() const noexcept
Definition agglomeration_handler.h:1029
AgglomerationHandler::cell_to_polytope_index
types::global_cell_index cell_to_polytope_index(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell) const
Definition agglomeration_handler.h:926
AgglomerationHandler::get_dof_handler
const DoFHandler< dim, spacedim > & get_dof_handler() const
Definition agglomeration_handler.h:968
AgglomerationHandler::get_slaves_of_idx
const std::vector< typename Triangulation< dim, spacedim >::active_cell_iterator > & get_slaves_of_idx(types::global_cell_index idx) const
Definition agglomeration_handler.h:978
AgglomerationHandler::reinit
const FEValues< dim, spacedim > & reinit(const AgglomerationIterator< dim, spacedim > &polytope) const
Definition agglomeration_handler.cc:500
AgglomerationHandler::recv_jxws
std::map< types::subdomain_id, std::map< std::pair< CellId, unsigned int >, std::vector< double > > > recv_jxws
Definition agglomeration_handler.h:499
AgglomerationHandler::output_dh
DoFHandler< dim, spacedim > output_dh
Definition agglomeration_handler.h:472
AgglomerationHandler::reinit_interface
std::pair< const FEValuesBase< dim, spacedim > &, const FEValuesBase< dim, spacedim > & > reinit_interface(const AgglomerationIterator< dim, spacedim > &polytope_in, const AgglomerationIterator< dim, spacedim > &neigh_polytope, const unsigned int local_in, const unsigned int local_outside) const
Definition agglomeration_handler.cc:565
AgglomerationHandler::recv_gradients
std::map< types::subdomain_id, std::map< std::pair< CellId, unsigned int >, std::vector< std::vector< Tensor< 1, spacedim > > > > > recv_gradients
Definition agglomeration_handler.h:514
AgglomerationHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
Definition agglomeration_handler.h:909
AgglomerationHandler::get_mapping
const Mapping< dim > & get_mapping() const
Definition agglomeration_handler.h:891
AgglomerationHandler::n_agglomerates
unsigned int n_agglomerates() const
Definition agglomeration_handler.h:1097
AgglomerationHandler::get_local_bboxes
const std::vector< BoundingBox< dim > > & get_local_bboxes() const
Definition agglomeration_handler.h:917
AgglomerationHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe() const
Definition agglomeration_handler.h:882
AgglomerationHandler::recv_values
std::map< types::subdomain_id, std::map< std::pair< CellId, unsigned int >, std::vector< std::vector< double > > > > recv_values
Definition agglomeration_handler.h:509
AgglomerationHandler::is_master_cell
bool is_master_cell(const CellIterator &cell) const
Definition agglomeration_handler.h:989
AgglomerationHandler::polytope_iterators
IteratorRange< agglomeration_iterator > polytope_iterators() const
Definition agglomeration_handler.h:1173
AgglomerationHandler::define_agglomerate
agglomeration_iterator define_agglomerate(const AgglomerationContainer &cells)
Definition agglomeration_handler.cc:46
AgglomerationHandler::agglo_dh
DoFHandler< dim, spacedim > agglo_dh
Definition agglomeration_handler.h:467
dealii::CellsAgglomerator::extract_agglomerates
const std::vector< std::vector< typename Triangulation< dim >::active_cell_iterator > > & extract_agglomerates()
Definition agglomerator.h:392
conditional_ostream.h
center
Point< 3 > center
dof_handler.h
dynamic_sparsity_pattern.h
fe_values.h
fe_dgq.h
level
unsigned int level
grid_tools.h
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Triangulation::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
ExcTrilinosError
static ::ExceptionBase & ExcTrilinosError(int arg1)
Assert
#define Assert(cond, exc)
AssertDimension
#define AssertDimension(dim1, dim2)
AssertThrow
#define AssertThrow(cond, exc)
DoFTools::make_flux_sparsity_pattern
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
update_values
update_values
update_normal_vectors
update_normal_vectors
update_JxW_values
update_JxW_values
update_gradients
update_gradients
update_quadrature_points
update_quadrature_points
iota_view.h
pow
Expression pow(const Expression &base, const Expression &exponent)
DoFTools::extract_locally_relevant_dofs
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
volume
double volume(const Triangulation< dim, spacedim > &tria)
diameter
double diameter(const Triangulation< dim, spacedim > &tria)
b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparsityTools::partition
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
SparsityTools::distribute_sparsity_pattern
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm mpi_communicator)
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Utilities::MPI::reduce
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
dealii::PolyUtils::internal::get_face_connectivity_of_cells
void get_face_connectivity_of_cells(const parallel::fullydistributed::Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &cell_connectivity, const std::vector< types::global_cell_index > locally_owned_cells)
Definition poly_utils.h:97
dealii::PolyUtils::internal::interpolate_to_fine_grid
void interpolate_to_fine_grid(const AgglomerationHandler< dim, spacedim > &agglomeration_handler, VectorType &dst, const VectorType &src)
Definition poly_utils.h:938
dealii::PolyUtils::internal::get_index
types::global_cell_index get_index(const std::vector< types::global_cell_index > &v, const types::global_cell_index index)
Definition poly_utils.h:83
dealii::PolyUtils::partition_locally_owned_regions
void partition_locally_owned_regions(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner)
Definition poly_utils.h:555
dealii::PolyUtils::fill_interpolation_matrix
void fill_interpolation_matrix(const AgglomerationHandler< dim, spacedim > &agglomeration_handler, MatrixType &interpolation_matrix)
Definition poly_utils.h:1277
dealii::PolyUtils::assemble_local_jumps_and_averages
void assemble_local_jumps_and_averages(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe_faces0, const FEValuesBase< dim > &fe_faces1, const double penalty_constant, const double h_f)
Definition poly_utils.h:1678
dealii::PolyUtils::extract_children_of_level
std::pair< std::vector< std::vector< unsigned int > >, std::vector< std::vector< typename Triangulation< boost::geometry::dimension< typename Rtree::indexable_type >::value >::active_cell_iterator > > > extract_children_of_level(const Rtree &tree, const unsigned int level)
Definition poly_utils.h:359
dealii::PolyUtils::interpolate_to_fine_grid
void interpolate_to_fine_grid(const AgglomerationHandler< dim, spacedim > &agglomeration_handler, VectorType &dst, const VectorType &src, const bool on_the_fly=true)
Definition poly_utils.h:1147
dealii::PolyUtils::compute_global_error
void compute_global_error(const AgglomerationHandler< dim > &agglomeration_handler, const VectorType &solution, const Function< dim, Number > &exact_solution, const std::vector< VectorTools::NormType > &norms, std::vector< double > &global_errors)
Definition poly_utils.h:1455
dealii::PolyUtils::compute_h_orthogonal
Number compute_h_orthogonal(const unsigned int face_index, const std::vector< typename Triangulation< dim >::active_face_iterator > &polygon_boundary, const Tensor< 1, dim > &deal_normal)
Definition poly_utils.h:401
dealii::PolyUtils::construct_agglomerated_levels
unsigned int construct_agglomerated_levels(const Triangulation< dim > &tria, std::vector< std::unique_ptr< AgglomerationHandler< dim > > > &agglomeration_handlers, const FE_DGQ< dim > &fe_dg, const Mapping< dim > &mapping, const unsigned int starting_tree_level)
Definition poly_utils.h:1569
dealii::PolyUtils::constexpr_pow
constexpr T constexpr_pow(T num, unsigned int pow)
Definition poly_utils.h:896
dealii::PolyUtils::assemble_dg_matrix
void assemble_dg_matrix(MatrixType &system_matrix, const FiniteElement< dim > &fe_dg, const AgglomerationHandler< dim > &ah)
Definition poly_utils.h:1808
dealii::PolyUtils::collect_cells_for_agglomeration
void collect_cells_for_agglomeration(const Triangulation< dim, spacedim > &tria, const std::vector< types::global_cell_index > &cell_idxs, std::vector< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_be_agglomerated)
Definition poly_utils.h:523
dealii::PolyUtils::assemble_local_jumps_and_averages_ghost
void assemble_local_jumps_and_averages_ghost(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe_faces0, const std::vector< std::vector< double > > &recv_values, const std::vector< std::vector< Tensor< 1, dim > > > &recv_gradients, const std::vector< double > &recv_jxws, const double penalty_constant, const double h_f)
Definition poly_utils.h:1738
dealii::PolyUtils::compute_quality_metrics
std::tuple< std::vector< double >, std::vector< double >, std::vector< double >, double > compute_quality_metrics(const AgglomerationHandler< dim > &ah)
Definition poly_utils.h:711
dealii::PolyUtils::assemble_dg_matrix_on_standard_mesh
void assemble_dg_matrix_on_standard_mesh(MatrixType &system_matrix, VectorType &system_rhs, const Mapping< dim > &mapping, const FiniteElement< dim > &fe_dg, const DoFHandler< dim > &dof_handler)
Definition poly_utils.h:2012
dealii::PolyUtils::write_to_matrix_market_format
void write_to_matrix_market_format(const std::string &filename, const std::string &matrix_name, const TrilinosWrappers::SparseMatrix &matrix)
Definition poly_utils.h:906
dealii::PolyUtils::export_polygon_to_csv_file
void export_polygon_to_csv_file(const AgglomerationHandler< dim > &agglomeration_handler, const std::string &filename)
Definition poly_utils.h:863
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
types::subdomain_id
unsigned int subdomain_id
types::global_cell_index
unsigned int global_cell_index
triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
comm
*braid_SplitCommworld & comm
point_conversion.h
sparse_matrix.h
sparsity_pattern.h
sparsity_tools.h
dealii::PolyUtils::Rtree_visitor
Definition poly_utils.h:156
dealii::PolyUtils::Rtree_visitor::translator
const Translator & translator
Definition poly_utils.h:202
dealii::PolyUtils::Rtree_visitor::level
size_t level
Definition poly_utils.h:207
dealii::PolyUtils::Rtree_visitor::operator()
void operator()(const InternalNode &node)
Definition poly_utils.h:265
dealii::PolyUtils::Rtree_visitor::agglomerates
std::vector< std::vector< typename Triangulation< boost::geometry::dimension< Box >::value >::active_cell_iterator > > & agglomerates
Definition poly_utils.h:229
dealii::PolyUtils::Rtree_visitor::Rtree_visitor
Rtree_visitor(const Translator &translator, unsigned int target_level, std::vector< std::vector< typename Triangulation< boost::geometry::dimension< Box >::value >::active_cell_iterator > > &boxes, std::vector< std::vector< unsigned int > > &csr)
Definition poly_utils.h:241
dealii::PolyUtils::Rtree_visitor::node_counter
size_t node_counter
Definition poly_utils.h:213
dealii::PolyUtils::Rtree_visitor::target_level
const size_t target_level
Definition poly_utils.h:220
dealii::PolyUtils::Rtree_visitor::row_ptr
std::vector< std::vector< unsigned int > > & row_ptr
Definition poly_utils.h:231
dealii::PolyUtils::Rtree_visitor::Leaf
typename boost::geometry::index::detail::rtree::leaf< Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag >::type Leaf
Definition poly_utils.h:179
dealii::PolyUtils::Rtree_visitor::next_level_leafs_processed
size_t next_level_leafs_processed
Definition poly_utils.h:215
dealii::PolyUtils::Rtree_visitor::InternalNode
typename boost::geometry::index::detail::rtree::internal_node< Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag >::type InternalNode
Definition poly_utils.h:168
trilinos_sparse_matrix.h
vector_tools_common.h