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_mpi_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_mpi_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 if (!agglomeration_handler.used_fe_collection())
1178 {
1179 // Original version: handle case without hp::FECollection
1180 const Triangulation<dim, spacedim> &tria =
1181 agglomeration_handler.get_triangulation();
1182 const Mapping<dim> &mapping = agglomeration_handler.get_mapping();
1183 const FiniteElement<dim, spacedim> &original_fe =
1184 agglomeration_handler.get_fe();
1185
1186 // We use DGQ (on tensor-product meshes) or DGP (on simplex meshes)
1187 // nodal elements of the same degree as the ones in the
1188 // agglomeration handler to interpolate the solution onto the finer
1189 // grid.
1190 std::unique_ptr<FiniteElement<dim>> output_fe;
1191 if (tria.all_reference_cells_are_hyper_cube())
1192 output_fe = std::make_unique<FE_DGQ<dim>>(original_fe.degree);
1193 else if (tria.all_reference_cells_are_simplex())
1194 output_fe =
1195 std::make_unique<FE_SimplexDGP<dim>>(original_fe.degree);
1196 else
1197 AssertThrow(false, ExcNotImplemented());
1198
1199 DoFHandler<dim> &output_dh =
1200 const_cast<DoFHandler<dim> &>(agglomeration_handler.output_dh);
1201 output_dh.reinit(tria);
1202 output_dh.distribute_dofs(*output_fe);
1203
1204 if constexpr (std::is_same_v<VectorType,
1205 TrilinosWrappers::MPI::Vector>)
1206 {
1207 const IndexSet &locally_owned_dofs =
1208 output_dh.locally_owned_dofs();
1209 dst.reinit(locally_owned_dofs);
1210 }
1211 else if constexpr (std::is_same_v<VectorType, Vector<NumberType>>)
1212 {
1213 dst.reinit(output_dh.n_dofs());
1214 }
1215 else
1216 {
1217 // PETSc, LA::d::v options not implemented.
1218 (void)agglomeration_handler;
1219 (void)dst;
1220 (void)src;
1221 AssertThrow(false, ExcNotImplemented());
1222 }
1223
1224
1225
1226 const unsigned int dofs_per_cell =
1227 agglomeration_handler.n_dofs_per_cell();
1228 const unsigned int output_dofs_per_cell =
1229 output_fe->n_dofs_per_cell();
1230 Quadrature<dim> quad(output_fe->get_unit_support_points());
1231 FEValues<dim> output_fe_values(mapping,
1232 *output_fe,
1233 quad,
1234 update_quadrature_points);
1235
1236 std::vector<types::global_dof_index> local_dof_indices(
1237 dofs_per_cell);
1238 std::vector<types::global_dof_index> local_dof_indices_output(
1239 output_dofs_per_cell);
1240
1241 const auto &bboxes = agglomeration_handler.get_local_bboxes();
1242 for (const auto &polytope :
1243 agglomeration_handler.polytope_iterators())
1244 {
1245 if (polytope->is_locally_owned())
1246 {
1247 polytope->get_dof_indices(local_dof_indices);
1248 const BoundingBox<dim> &box = bboxes[polytope->index()];
1249
1250 const auto &deal_cells =
1251 polytope->get_agglomerate(); // fine deal.II cells
1252 for (const auto &cell : deal_cells)
1253 {
1254 const auto slave_output = cell->as_dof_handler_iterator(
1255 agglomeration_handler.output_dh);
1256 slave_output->get_dof_indices(local_dof_indices_output);
1257 output_fe_values.reinit(slave_output);
1258
1259 const auto &qpoints =
1260 output_fe_values.get_quadrature_points();
1261
1262 for (unsigned int j = 0; j < output_dofs_per_cell; ++j)
1263 {
1264 const auto &ref_qpoint =
1265 box.real_to_unit(qpoints[j]);
1266 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1267 dst(local_dof_indices_output[j]) +=
1268 src(local_dof_indices[i]) *
1269 original_fe.shape_value(i, ref_qpoint);
1270 }
1271 }
1272 }
1273 }
1274 }
1275 else
1276 {
1277 // Handle the hp::FECollection case
1278 const Triangulation<dim, spacedim> &tria =
1279 agglomeration_handler.get_triangulation();
1280 const Mapping<dim> &mapping = agglomeration_handler.get_mapping();
1281 const hp::FECollection<dim, spacedim> &original_fe_collection =
1282 agglomeration_handler.get_fe_collection();
1283
1284 // We use DGQ (on tensor-product meshes) or DGP (on simplex meshes)
1285 // nodal elements of the same degree as the ones in the
1286 // agglomeration handler to interpolate the solution onto the finer
1287 // grid.
1288 hp::FECollection<dim, spacedim> output_fe_collection;
1289
1290 Assert(original_fe_collection[0].n_components() >= 1,
1291 ExcMessage("Invalid FE: must have at least one component."));
1292 if (original_fe_collection[0].n_components() == 1)
1293 {
1294 // Scalar case
1295 for (unsigned int i = 0; i < original_fe_collection.size(); ++i)
1296 {
1297 std::unique_ptr<FiniteElement<dim>> output_fe;
1298 if (tria.all_reference_cells_are_hyper_cube())
1299 output_fe = std::make_unique<FE_DGQ<dim>>(
1300 original_fe_collection[i].degree);
1301 else if (tria.all_reference_cells_are_simplex())
1302 output_fe = std::make_unique<FE_SimplexDGP<dim>>(
1303 original_fe_collection[i].degree);
1304 else
1305 AssertThrow(false, ExcNotImplemented());
1306 output_fe_collection.push_back(*output_fe);
1307 }
1308 }
1309 else if (original_fe_collection[0].n_components() > 1)
1310 {
1311 // System case
1312 for (unsigned int i = 0; i < original_fe_collection.size(); ++i)
1313 {
1314 std::vector<const FiniteElement<dim, spacedim> *>
1315 base_elements;
1316 std::vector<unsigned int> multiplicities;
1317 for (unsigned int b = 0;
1318 b < original_fe_collection[i].n_base_elements();
1319 ++b)
1320 {
1321 if (dynamic_cast<const FE_Nothing<dim> *>(
1322 &original_fe_collection[i].base_element(b)))
1323 base_elements.push_back(
1324 new FE_Nothing<dim, spacedim>());
1325 else
1326 {
1327 if (tria.all_reference_cells_are_hyper_cube())
1328 base_elements.push_back(new FE_DGQ<dim, spacedim>(
1329 original_fe_collection[i]
1330 .base_element(b)
1331 .degree));
1332 else if (tria.all_reference_cells_are_simplex())
1333 base_elements.push_back(
1334 new FE_SimplexDGP<dim, spacedim>(
1335 original_fe_collection[i]
1336 .base_element(b)
1337 .degree));
1338 else
1339 AssertThrow(false, ExcNotImplemented());
1340 }
1341 multiplicities.push_back(
1342 original_fe_collection[i].element_multiplicity(b));
1343 }
1344
1345 FESystem<dim, spacedim> output_fe_system(base_elements,
1346 multiplicities);
1347 for (const auto *ptr : base_elements)
1348 delete ptr;
1349 output_fe_collection.push_back(output_fe_system);
1350 }
1351 }
1352
1353
1354 DoFHandler<dim> &output_dh =
1355 const_cast<DoFHandler<dim> &>(agglomeration_handler.output_dh);
1356 output_dh.reinit(tria);
1357 for (const auto &polytope :
1358 agglomeration_handler.polytope_iterators())
1359 {
1360 if (polytope->is_locally_owned())
1361 {
1362 const auto &deal_cells =
1363 polytope->get_agglomerate(); // fine deal.II cells
1364 const unsigned int active_fe_idx =
1365 polytope->active_fe_index();
1366
1367 for (const auto &cell : deal_cells)
1368 {
1369 const typename DoFHandler<dim>::active_cell_iterator
1370 slave_cell_dh_iterator =
1371 cell->as_dof_handler_iterator(output_dh);
1372 slave_cell_dh_iterator->set_active_fe_index(
1373 active_fe_idx);
1374 }
1375 }
1376 }
1377 output_dh.distribute_dofs(output_fe_collection);
1378
1379 if constexpr (std::is_same_v<VectorType,
1380 TrilinosWrappers::MPI::Vector>)
1381 {
1382 const IndexSet &locally_owned_dofs =
1383 output_dh.locally_owned_dofs();
1384 dst.reinit(locally_owned_dofs);
1385 }
1386 else if constexpr (std::is_same_v<VectorType, Vector<NumberType>>)
1387 {
1388 dst.reinit(output_dh.n_dofs());
1389 }
1390 else
1391 {
1392 // PETSc, LA::d::v options not implemented.
1393 (void)agglomeration_handler;
1394 (void)dst;
1395 (void)src;
1396 AssertThrow(false, ExcNotImplemented());
1397 }
1398
1399 const auto &bboxes = agglomeration_handler.get_local_bboxes();
1400 for (const auto &polytope :
1401 agglomeration_handler.polytope_iterators())
1402 {
1403 if (polytope->is_locally_owned())
1404 {
1405 const unsigned int active_fe_idx =
1406 polytope->active_fe_index();
1407 const unsigned int dofs_per_cell =
1408 polytope->get_fe().dofs_per_cell;
1409 const unsigned int output_dofs_per_cell =
1410 output_fe_collection[active_fe_idx].n_dofs_per_cell();
1411 Quadrature<dim> quad(output_fe_collection[active_fe_idx]
1412 .get_unit_support_points());
1413 FEValues<dim> output_fe_values(
1414 mapping,
1415 output_fe_collection[active_fe_idx],
1416 quad,
1417 update_quadrature_points);
1418 std::vector<types::global_dof_index> local_dof_indices(
1419 dofs_per_cell);
1420 std::vector<types::global_dof_index>
1421 local_dof_indices_output(output_dofs_per_cell);
1422
1423 polytope->get_dof_indices(local_dof_indices);
1424 const BoundingBox<dim> &box = bboxes[polytope->index()];
1425
1426 const auto &deal_cells =
1427 polytope->get_agglomerate(); // fine deal.II cells
1428 for (const auto &cell : deal_cells)
1429 {
1430 const auto slave_output = cell->as_dof_handler_iterator(
1431 agglomeration_handler.output_dh);
1432 slave_output->get_dof_indices(local_dof_indices_output);
1433 output_fe_values.reinit(slave_output);
1434
1435 const auto &qpoints =
1436 output_fe_values.get_quadrature_points();
1437
1438 for (unsigned int j = 0; j < output_dofs_per_cell; ++j)
1439 {
1440 const unsigned int component_idx_of_this_dof =
1441 slave_output->get_fe()
1442 .system_to_component_index(j)
1443 .first;
1444 const auto &ref_qpoint =
1445 box.real_to_unit(qpoints[j]);
1446 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1447 dst(local_dof_indices_output[j]) +=
1448 src(local_dof_indices[i]) *
1449 original_fe_collection[active_fe_idx]
1450 .shape_value_component(
1451 i, ref_qpoint, component_idx_of_this_dof);
1452 }
1453 }
1454 }
1455 }
1456 }
1457 }
1458 }
1459
1460
1461
1469 template <int dim, int spacedim, typename MatrixType>
1470 void
1471 fill_interpolation_matrix(
1472 const AgglomerationHandler<dim, spacedim> &agglomeration_handler,
1473 MatrixType &interpolation_matrix)
1474 {
1475 Assert((dim == spacedim), ExcNotImplemented());
1476
1477 using NumberType = typename MatrixType::value_type;
1478 constexpr bool is_trilinos_matrix =
1479 std::is_same_v<MatrixType, TrilinosWrappers::MPI::Vector>;
1480
1481 [[maybe_unused]]
1482 typename std::conditional_t<!is_trilinos_matrix, SparsityPattern, void *>
1483 sp;
1484
1485 // Get some info from the handler
1486 const DoFHandler<dim, spacedim> &agglo_dh = agglomeration_handler.agglo_dh;
1487
1488 DoFHandler<dim, spacedim> *output_dh =
1489 const_cast<DoFHandler<dim, spacedim> *>(&agglomeration_handler.output_dh);
1490 const Mapping<dim, spacedim> &mapping = agglomeration_handler.get_mapping();
1491 const FiniteElement<dim, spacedim> &fe = agglomeration_handler.get_fe();
1492 const Triangulation<dim, spacedim> &tria =
1493 agglomeration_handler.get_triangulation();
1494 const auto &bboxes = agglomeration_handler.get_local_bboxes();
1495
1496 // Setup an auxiliary DoFHandler for output purposes
1497 output_dh->reinit(tria);
1498 output_dh->distribute_dofs(fe);
1499
1500 const IndexSet &locally_owned_dofs = output_dh->locally_owned_dofs();
1501 const IndexSet locally_relevant_dofs =
1502 DoFTools::extract_locally_relevant_dofs(*output_dh);
1503
1504 const IndexSet &locally_owned_dofs_agglo = agglo_dh.locally_owned_dofs();
1505
1506
1507 DynamicSparsityPattern dsp(output_dh->n_dofs(),
1508 agglo_dh.n_dofs(),
1509 locally_relevant_dofs);
1510
1511 std::vector<types::global_dof_index> agglo_dof_indices(fe.dofs_per_cell);
1512 std::vector<types::global_dof_index> standard_dof_indices(fe.dofs_per_cell);
1513 std::vector<types::global_dof_index> output_dof_indices(fe.dofs_per_cell);
1514
1515 Quadrature<dim> quad(fe.get_unit_support_points());
1516 FEValues<dim, spacedim> output_fe_values(mapping,
1517 fe,
1518 quad,
1519 update_quadrature_points);
1520
1521 for (const auto &cell : agglo_dh.active_cell_iterators())
1522 if (cell->is_locally_owned())
1523 {
1524 if (agglomeration_handler.is_master_cell(cell))
1525 {
1526 auto slaves = agglomeration_handler.get_slaves_of_idx(
1527 cell->active_cell_index());
1528 slaves.emplace_back(cell);
1529
1530 cell->get_dof_indices(agglo_dof_indices);
1531
1532 for (const auto &slave : slaves)
1533 {
1534 // addd master-slave relationship
1535 const auto slave_output =
1536 slave->as_dof_handler_iterator(*output_dh);
1537 slave_output->get_dof_indices(output_dof_indices);
1538 for (const auto row : output_dof_indices)
1539 dsp.add_entries(row,
1540 agglo_dof_indices.begin(),
1541 agglo_dof_indices.end());
1542 }
1543 }
1544 }
1545
1546
1547 const auto assemble_interpolation_matrix = [&]() {
1548 FullMatrix<NumberType> local_matrix(fe.dofs_per_cell, fe.dofs_per_cell);
1549 std::vector<Point<dim>> reference_q_points(fe.dofs_per_cell);
1550
1551 // Dummy AffineConstraints, only needed for loc2glb
1552 AffineConstraints<NumberType> c;
1553 c.close();
1554
1555 for (const auto &cell : agglo_dh.active_cell_iterators())
1556 if (cell->is_locally_owned())
1557 {
1558 if (agglomeration_handler.is_master_cell(cell))
1559 {
1560 auto slaves = agglomeration_handler.get_slaves_of_idx(
1561 cell->active_cell_index());
1562 slaves.emplace_back(cell);
1563
1564 cell->get_dof_indices(agglo_dof_indices);
1565
1566 const types::global_cell_index polytope_index =
1567 agglomeration_handler.cell_to_polytope_index(cell);
1568
1569 // Get the box of this agglomerate.
1570 const BoundingBox<dim> &box = bboxes[polytope_index];
1571
1572 for (const auto &slave : slaves)
1573 {
1574 // add master-slave relationship
1575 const auto slave_output =
1576 slave->as_dof_handler_iterator(*output_dh);
1577
1578 slave_output->get_dof_indices(output_dof_indices);
1579 output_fe_values.reinit(slave_output);
1580
1581 local_matrix = 0.;
1582
1583 const auto &q_points =
1584 output_fe_values.get_quadrature_points();
1585 for (const auto i : output_fe_values.dof_indices())
1586 {
1587 const auto &p = box.real_to_unit(q_points[i]);
1588 for (const auto j : output_fe_values.dof_indices())
1589 {
1590 local_matrix(i, j) = fe.shape_value(j, p);
1591 }
1592 }
1593 c.distribute_local_to_global(local_matrix,
1594 output_dof_indices,
1595 agglo_dof_indices,
1596 interpolation_matrix);
1597 }
1598 }
1599 }
1600 };
1601
1602
1603 if constexpr (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix>)
1604 {
1605 const MPI_Comm &communicator = tria.get_mpi_communicator();
1606 SparsityTools::distribute_sparsity_pattern(dsp,
1607 locally_owned_dofs,
1608 communicator,
1609 locally_relevant_dofs);
1610
1611 interpolation_matrix.reinit(locally_owned_dofs,
1612 locally_owned_dofs_agglo,
1613 dsp,
1614 communicator);
1615 assemble_interpolation_matrix();
1616 }
1617 else if constexpr (std::is_same_v<MatrixType, SparseMatrix<NumberType>>)
1618 {
1619 sp.copy_from(dsp);
1620 interpolation_matrix.reinit(sp);
1621 assemble_interpolation_matrix();
1622 }
1623 else
1624 {
1625 // PETSc, LA::d::v options not implemented.
1626 (void)agglomeration_handler;
1627 AssertThrow(false, ExcNotImplemented());
1628 }
1629
1630 // If tria is distributed
1631 if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1632 &tria) != nullptr)
1633 interpolation_matrix.compress(VectorOperation::add);
1634 }
1635
1636
1637
1647 template <int dim, typename Number, typename VectorType>
1648 void
1649 compute_global_error(const AgglomerationHandler<dim> &agglomeration_handler,
1650 const VectorType &solution,
1651 const Function<dim, Number> &exact_solution,
1652 const std::vector<VectorTools::NormType> &norms,
1653 std::vector<double> &global_errors)
1654 {
1655 Assert(solution.size() > 0,
1656 ExcNotImplemented(
1657 "Solution vector must be non-empty upon calling this function."));
1658 Assert(std::any_of(norms.cbegin(),
1659 norms.cend(),
1660 [](VectorTools::NormType norm_type) {
1661 return (norm_type ==
1662 VectorTools::NormType::H1_seminorm ||
1663 norm_type == VectorTools::NormType::L2_norm);
1664 }),
1665 ExcMessage("Norm type not supported"));
1666 global_errors.resize(norms.size());
1667 std::fill(global_errors.begin(), global_errors.end(), 0.);
1668
1669 // Vector storing errors local to the current processor.
1670 std::vector<double> local_errors(norms.size());
1671 std::fill(local_errors.begin(), local_errors.end(), 0.);
1672
1673 // Get some info from the handler
1674 const unsigned int dofs_per_cell = agglomeration_handler.n_dofs_per_cell();
1675
1676 const bool compute_semi_H1 =
1677 std::any_of(norms.cbegin(),
1678 norms.cend(),
1679 [](VectorTools::NormType norm_type) {
1680 return norm_type == VectorTools::NormType::H1_seminorm;
1681 });
1682
1683 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1684 for (const auto &polytope : agglomeration_handler.polytope_iterators())
1685 {
1686 if (polytope->is_locally_owned())
1687 {
1688 const auto &agglo_values = agglomeration_handler.reinit(polytope);
1689 polytope->get_dof_indices(local_dof_indices);
1690
1691 const auto &q_points = agglo_values.get_quadrature_points();
1692 const unsigned int n_qpoints = q_points.size();
1693 std::vector<double> analyical_sol_at_qpoints(n_qpoints);
1694 exact_solution.value_list(q_points, analyical_sol_at_qpoints);
1695 std::vector<Tensor<1, dim>> grad_analyical_sol_at_qpoints(
1696 n_qpoints);
1697
1698 if (compute_semi_H1)
1699 exact_solution.gradient_list(q_points,
1700 grad_analyical_sol_at_qpoints);
1701
1702 for (unsigned int q_index : agglo_values.quadrature_point_indices())
1703 {
1704 double solution_at_qpoint = 0.;
1705 Tensor<1, dim> grad_solution_at_qpoint;
1706 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1707 {
1708 solution_at_qpoint += solution(local_dof_indices[i]) *
1709 agglo_values.shape_value(i, q_index);
1710
1711 if (compute_semi_H1)
1712 grad_solution_at_qpoint +=
1713 solution(local_dof_indices[i]) *
1714 agglo_values.shape_grad(i, q_index);
1715 }
1716 // L2
1717 local_errors[0] += std::pow((analyical_sol_at_qpoints[q_index] -
1718 solution_at_qpoint),
1719 2) *
1720 agglo_values.JxW(q_index);
1721
1722 // H1 seminorm
1723 if (compute_semi_H1)
1724 for (unsigned int d = 0; d < dim; ++d)
1725 local_errors[1] +=
1726 std::pow((grad_analyical_sol_at_qpoints[q_index][d] -
1727 grad_solution_at_qpoint[d]),
1728 2) *
1729 agglo_values.JxW(q_index);
1730 }
1731 }
1732 }
1733
1734 // Perform reduction and take sqrt of each error
1735 global_errors[0] = Utilities::MPI::reduce<double>(
1736 local_errors[0],
1737 agglomeration_handler.get_triangulation().get_mpi_communicator(),
1738 [](const double a, const double b) { return a + b; });
1739
1740 global_errors[0] = std::sqrt(global_errors[0]);
1741
1742 if (compute_semi_H1)
1743 {
1744 global_errors[1] = Utilities::MPI::reduce<double>(
1745 local_errors[1],
1746 agglomeration_handler.get_triangulation().get_mpi_communicator(),
1747 [](const double a, const double b) { return a + b; });
1748 global_errors[1] = std::sqrt(global_errors[1]);
1749 }
1750 }
1751
1752
1753
1761 template <int dim>
1762 unsigned int
1763 construct_agglomerated_levels(
1764 const Triangulation<dim> &tria,
1765 std::vector<std::unique_ptr<AgglomerationHandler<dim>>>
1766 &agglomeration_handlers,
1767 const FE_DGQ<dim> &fe_dg,
1768 const Mapping<dim> &mapping,
1769 const unsigned int starting_tree_level)
1770 {
1771 const auto parallel_tria =
1772 dynamic_cast<const parallel::TriangulationBase<dim> *>(&tria);
1773
1774 GridTools::Cache<dim> cached_tria(tria);
1775 Assert(parallel_tria->n_active_cells() > 0, ExcInternalError());
1776
1777 const MPI_Comm comm = parallel_tria->get_mpi_communicator();
1778 ConditionalOStream pcout(std::cout,
1779 (Utilities::MPI::this_mpi_process(comm) == 0));
1780
1781 // Start building R-tree
1782 namespace bgi = boost::geometry::index;
1783 static constexpr unsigned int max_elem_per_node =
1784 constexpr_pow(2, dim); // 2^dim
1785 std::vector<std::pair<BoundingBox<dim>,
1786 typename Triangulation<dim>::active_cell_iterator>>
1787 boxes(parallel_tria->n_locally_owned_active_cells());
1788 unsigned int i = 0;
1789 for (const auto &cell : parallel_tria->active_cell_iterators())
1790 if (cell->is_locally_owned())
1791 boxes[i++] = std::make_pair(mapping.get_bounding_box(cell), cell);
1792
1793 auto tree = pack_rtree<bgi::rstar<max_elem_per_node>>(boxes);
1794 Assert(n_levels(tree) >= 2, ExcMessage("At least two levels are needed."));
1795 pcout << "Total number of available levels: " << n_levels(tree)
1796 << std::endl;
1797
1798 pcout << "Starting level: " << starting_tree_level << std::endl;
1799 const unsigned int total_tree_levels =
1800 n_levels(tree) - starting_tree_level + 1;
1801
1802
1803 // Resize the agglomeration handlers to the right size
1804
1805 agglomeration_handlers.resize(total_tree_levels);
1806 // Loop through the available levels and set AgglomerationHandlers up.
1807 for (unsigned int extraction_level = starting_tree_level;
1808 extraction_level <= n_levels(tree);
1809 ++extraction_level)
1810 {
1811 agglomeration_handlers[extraction_level - starting_tree_level] =
1812 std::make_unique<AgglomerationHandler<dim>>(cached_tria);
1813 CellsAgglomerator<dim, decltype(tree)> agglomerator{tree,
1814 extraction_level};
1815 const auto agglomerates = agglomerator.extract_agglomerates();
1816 agglomeration_handlers[extraction_level - starting_tree_level]
1817 ->connect_hierarchy(agglomerator);
1818
1819 // Flag elements for agglomeration
1820 unsigned int agglo_index = 0;
1821 for (unsigned int i = 0; i < agglomerates.size(); ++i)
1822 {
1823 const auto &agglo = agglomerates[i]; // i-th agglomerate
1824 for (const auto &el : agglo)
1825 {
1826 el->set_material_id(agglo_index);
1827 }
1828 ++agglo_index;
1829 }
1830
1831 const unsigned int n_local_agglomerates = agglo_index;
1832 unsigned int total_agglomerates =
1833 Utilities::MPI::sum(n_local_agglomerates, comm);
1834 pcout << "Total agglomerates per (tree) level: " << extraction_level
1835 << ": " << total_agglomerates << std::endl;
1836
1837
1838 // Now, perform agglomeration within each locally owned partition
1839 std::vector<
1840 std::vector<typename Triangulation<dim>::active_cell_iterator>>
1841 cells_per_subdomain(n_local_agglomerates);
1842 for (const auto &cell : parallel_tria->active_cell_iterators())
1843 if (cell->is_locally_owned())
1844 cells_per_subdomain[cell->material_id()].push_back(cell);
1845
1846 // For every subdomain, agglomerate elements together
1847 for (std::size_t i = 0; i < cells_per_subdomain.size(); ++i)
1848 agglomeration_handlers[extraction_level - starting_tree_level]
1849 ->define_agglomerate(cells_per_subdomain[i]);
1850
1851 agglomeration_handlers[extraction_level - starting_tree_level]
1852 ->initialize_fe_values(QGauss<dim>(fe_dg.degree + 1),
1853 update_values | update_gradients |
1854 update_JxW_values | update_quadrature_points,
1855 QGauss<dim - 1>(fe_dg.degree + 1),
1856 update_JxW_values);
1857 agglomeration_handlers[extraction_level - starting_tree_level]
1858 ->distribute_agglomerated_dofs(fe_dg);
1859 }
1860
1861 return total_tree_levels;
1862 }
1863
1864
1865
1870 template <int dim>
1871 void
1872 assemble_local_jumps_and_averages(FullMatrix<double> &M11,
1873 FullMatrix<double> &M12,
1874 FullMatrix<double> &M21,
1875 FullMatrix<double> &M22,
1876 const FEValuesBase<dim> &fe_faces0,
1877 const FEValuesBase<dim> &fe_faces1,
1878 const double penalty_constant,
1879 const double h_f)
1880 {
1881 const std::vector<Tensor<1, dim>> &normals = fe_faces0.get_normal_vectors();
1882 const unsigned int dofs_per_cell =
1883 M11.m(); // size of local matrices equals the #DoFs
1884 for (unsigned int q_index : fe_faces0.quadrature_point_indices())
1885 {
1886 const Tensor<1, dim> &normal = normals[q_index];
1887 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1888 {
1889 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1890 {
1891 M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normal *
1892 fe_faces0.shape_value(j, q_index) -
1893 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
1894 fe_faces0.shape_value(i, q_index) +
1895 (penalty_constant / h_f) *
1896 fe_faces0.shape_value(i, q_index) *
1897 fe_faces0.shape_value(j, q_index)) *
1898 fe_faces0.JxW(q_index);
1899 M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normal *
1900 fe_faces1.shape_value(j, q_index) -
1901 0.5 * fe_faces1.shape_grad(j, q_index) * normal *
1902 fe_faces0.shape_value(i, q_index) -
1903 (penalty_constant / h_f) *
1904 fe_faces0.shape_value(i, q_index) *
1905 fe_faces1.shape_value(j, q_index)) *
1906 fe_faces1.JxW(q_index);
1907 M21(i, j) += (-0.5 * fe_faces1.shape_grad(i, q_index) * normal *
1908 fe_faces0.shape_value(j, q_index) +
1909 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
1910 fe_faces1.shape_value(i, q_index) -
1911 (penalty_constant / h_f) *
1912 fe_faces1.shape_value(i, q_index) *
1913 fe_faces0.shape_value(j, q_index)) *
1914 fe_faces1.JxW(q_index);
1915 M22(i, j) += (0.5 * fe_faces1.shape_grad(i, q_index) * normal *
1916 fe_faces1.shape_value(j, q_index) +
1917 0.5 * fe_faces1.shape_grad(j, q_index) * normal *
1918 fe_faces1.shape_value(i, q_index) +
1919 (penalty_constant / h_f) *
1920 fe_faces1.shape_value(i, q_index) *
1921 fe_faces1.shape_value(j, q_index)) *
1922 fe_faces1.JxW(q_index);
1923 }
1924 }
1925 }
1926 }
1927
1930 template <int dim>
1931 void
1932 assemble_local_jumps_and_averages_ghost(
1933 FullMatrix<double> &M11,
1934 FullMatrix<double> &M12,
1935 FullMatrix<double> &M21,
1936 FullMatrix<double> &M22,
1937 const FEValuesBase<dim> &fe_faces0,
1938 const std::vector<std::vector<double>> &recv_values,
1939 const std::vector<std::vector<Tensor<1, dim>>> &recv_gradients,
1940 const std::vector<double> &recv_jxws,
1941 const double penalty_constant,
1942 const double h_f)
1943 {
1944 Assert(
1945 (recv_values.size() > 0 && recv_gradients.size() && recv_jxws.size()),
1946 ExcMessage(
1947 "Not possible to assemble jumps and averages at a ghosted interface."));
1948 const unsigned int dofs_per_cell = M11.m();
1949 const std::vector<Tensor<1, dim>> &normals = fe_faces0.get_normal_vectors();
1950 for (unsigned int q_index : fe_faces0.quadrature_point_indices())
1951 {
1952 const Tensor<1, dim> &normal = normals[q_index];
1953 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1954 {
1955 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1956 {
1957 M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normal *
1958 fe_faces0.shape_value(j, q_index) -
1959 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
1960 fe_faces0.shape_value(i, q_index) +
1961 (penalty_constant / h_f) *
1962 fe_faces0.shape_value(i, q_index) *
1963 fe_faces0.shape_value(j, q_index)) *
1964 fe_faces0.JxW(q_index);
1965 M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normal *
1966 recv_values[j][q_index] -
1967 0.5 * recv_gradients[j][q_index] * normal *
1968 fe_faces0.shape_value(i, q_index) -
1969 (penalty_constant / h_f) *
1970 fe_faces0.shape_value(i, q_index) *
1971 recv_values[j][q_index]) *
1972 recv_jxws[q_index];
1973 M21(i, j) +=
1974 (-0.5 * recv_gradients[i][q_index] * normal *
1975 fe_faces0.shape_value(j, q_index) +
1976 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
1977 recv_values[i][q_index] -
1978 (penalty_constant / h_f) * recv_values[i][q_index] *
1979 fe_faces0.shape_value(j, q_index)) *
1980 recv_jxws[q_index];
1981 M22(i, j) +=
1982 (0.5 * recv_gradients[i][q_index] * normal *
1983 recv_values[j][q_index] +
1984 0.5 * recv_gradients[j][q_index] * normal *
1985 recv_values[i][q_index] +
1986 (penalty_constant / h_f) * recv_values[i][q_index] *
1987 recv_values[j][q_index]) *
1988 recv_jxws[q_index];
1989 }
1990 }
1991 }
1992 }
1993
1994
2000 template <int dim, typename MatrixType>
2001 void
2002 assemble_dg_matrix(MatrixType &system_matrix,
2003 const FiniteElement<dim> &fe_dg,
2004 const AgglomerationHandler<dim> &ah)
2005 {
2006 static_assert(
2007 (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix> ||
2008 std::is_same_v<MatrixType,
2009 SparseMatrix<typename MatrixType::value_type>>));
2010
2011 Assert((dynamic_cast<const FE_DGQ<dim> *>(&fe_dg) ||
2012 dynamic_cast<const FE_DGP<dim> *>(&fe_dg) ||
2013 dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_dg)),
2014 ExcMessage("FE type not supported."));
2015
2016 AffineConstraints constraints;
2017 constraints.close();
2018 const double penalty_constant =
2019 10 * (fe_dg.degree + dim) * (fe_dg.degree + 1);
2020 TrilinosWrappers::SparsityPattern dsp;
2021 const_cast<AgglomerationHandler<dim> &>(ah)
2022 .create_agglomeration_sparsity_pattern(dsp);
2023 system_matrix.reinit(dsp);
2024 const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell();
2025 FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
2026 FullMatrix<double> M11(dofs_per_cell, dofs_per_cell);
2027 FullMatrix<double> M12(dofs_per_cell, dofs_per_cell);
2028 FullMatrix<double> M21(dofs_per_cell, dofs_per_cell);
2029 FullMatrix<double> M22(dofs_per_cell, dofs_per_cell);
2030 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2031 std::vector<types::global_dof_index> local_dof_indices_neighbor(
2032 dofs_per_cell);
2033
2034 for (const auto &polytope : ah.polytope_iterators())
2035 {
2036 if (polytope->is_locally_owned())
2037 {
2038 cell_matrix = 0.;
2039 const auto &agglo_values = ah.reinit(polytope);
2040 for (unsigned int q_index : agglo_values.quadrature_point_indices())
2041 {
2042 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2043 {
2044 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2045 {
2046 cell_matrix(i, j) +=
2047 agglo_values.shape_grad(i, q_index) *
2048 agglo_values.shape_grad(j, q_index) *
2049 agglo_values.JxW(q_index);
2050 }
2051 }
2052 }
2053 // get volumetric DoFs
2054 polytope->get_dof_indices(local_dof_indices);
2055 // Assemble face terms
2056 unsigned int n_faces = polytope->n_faces();
2057 const double h_f = polytope->diameter();
2058 for (unsigned int f = 0; f < n_faces; ++f)
2059 {
2060 if (polytope->at_boundary(f))
2061 {
2062 // Get normal vectors seen from each agglomeration.
2063 const auto &fe_face = ah.reinit(polytope, f);
2064 const auto &normals = fe_face.get_normal_vectors();
2065 for (unsigned int q_index :
2066 fe_face.quadrature_point_indices())
2067 {
2068 const Tensor<1, dim> &normal = normals[q_index];
2069 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2070 {
2071 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2072 {
2073 cell_matrix(i, j) +=
2074 (-fe_face.shape_value(i, q_index) *
2075 fe_face.shape_grad(j, q_index) * normal -
2076 fe_face.shape_grad(i, q_index) * normal *
2077 fe_face.shape_value(j, q_index) +
2078 (penalty_constant / h_f) *
2079 fe_face.shape_value(i, q_index) *
2080 fe_face.shape_value(j, q_index)) *
2081 fe_face.JxW(q_index);
2082 }
2083 }
2084 }
2085 }
2086 else
2087 {
2088 const auto &neigh_polytope = polytope->neighbor(f);
2089 if (polytope->id() < neigh_polytope->id())
2090 {
2091 unsigned int nofn =
2092 polytope->neighbor_of_agglomerated_neighbor(f);
2093 Assert(neigh_polytope->neighbor(nofn)->id() ==
2094 polytope->id(),
2095 ExcMessage("Mismatch."));
2096 const auto &fe_faces = ah.reinit_interface(
2097 polytope, neigh_polytope, f, nofn);
2098 const auto &fe_faces0 = fe_faces.first;
2099 if (neigh_polytope->is_locally_owned())
2100 {
2101 // use both fevalues
2102 const auto &fe_faces1 = fe_faces.second;
2103 M11 = 0.;
2104 M12 = 0.;
2105 M21 = 0.;
2106 M22 = 0.;
2107 assemble_local_jumps_and_averages(M11,
2108 M12,
2109 M21,
2110 M22,
2111 fe_faces0,
2112 fe_faces1,
2113 penalty_constant,
2114 h_f);
2115 // distribute DoFs accordingly
2116 // fluxes
2117 neigh_polytope->get_dof_indices(
2118 local_dof_indices_neighbor);
2119 constraints.distribute_local_to_global(
2120 M11, local_dof_indices, system_matrix);
2121 constraints.distribute_local_to_global(
2122 M12,
2123 local_dof_indices,
2124 local_dof_indices_neighbor,
2125 system_matrix);
2126 constraints.distribute_local_to_global(
2127 M21,
2128 local_dof_indices_neighbor,
2129 local_dof_indices,
2130 system_matrix);
2131 constraints.distribute_local_to_global(
2132 M22, local_dof_indices_neighbor, system_matrix);
2133 }
2134 else
2135 {
2136 // neigh polytope is ghosted, so retrieve necessary
2137 // metadata.
2138 types::subdomain_id neigh_rank =
2139 neigh_polytope->subdomain_id();
2140 const auto &recv_jxws =
2141 ah.recv_jxws.at(neigh_rank)
2142 .at({neigh_polytope->id(), nofn});
2143 const auto &recv_values =
2144 ah.recv_values.at(neigh_rank)
2145 .at({neigh_polytope->id(), nofn});
2146 const auto &recv_gradients =
2147 ah.recv_gradients.at(neigh_rank)
2148 .at({neigh_polytope->id(), nofn});
2149 M11 = 0.;
2150 M12 = 0.;
2151 M21 = 0.;
2152 M22 = 0.;
2153 // there's no FEFaceValues on the other side (it's
2154 // ghosted), so we just pass the actual data we have
2155 // recevied from the neighboring ghosted polytope
2156 assemble_local_jumps_and_averages_ghost(
2157 M11,
2158 M12,
2159 M21,
2160 M22,
2161 fe_faces0,
2162 recv_values,
2163 recv_gradients,
2164 recv_jxws,
2165 penalty_constant,
2166 h_f);
2167 // distribute DoFs accordingly
2168 // fluxes
2169 neigh_polytope->get_dof_indices(
2170 local_dof_indices_neighbor);
2171 constraints.distribute_local_to_global(
2172 M11, local_dof_indices, system_matrix);
2173 constraints.distribute_local_to_global(
2174 M12,
2175 local_dof_indices,
2176 local_dof_indices_neighbor,
2177 system_matrix);
2178 constraints.distribute_local_to_global(
2179 M21,
2180 local_dof_indices_neighbor,
2181 local_dof_indices,
2182 system_matrix);
2183 constraints.distribute_local_to_global(
2184 M22, local_dof_indices_neighbor, system_matrix);
2185 } // ghosted polytope case
2186 } // only once
2187 } // internal face
2188 } // face loop
2189 constraints.distribute_local_to_global(cell_matrix,
2190 local_dof_indices,
2191 system_matrix);
2192 } // locally owned polytopes
2193 }
2194 system_matrix.compress(VectorOperation::add);
2195 }
2196
2197
2198
2204 template <int dim, typename MatrixType, typename VectorType>
2205 void
2206 assemble_dg_matrix_on_standard_mesh(MatrixType &system_matrix,
2207 VectorType &system_rhs,
2208 const Mapping<dim> &mapping,
2209 const FiniteElement<dim> &fe_dg,
2210 const DoFHandler<dim> &dof_handler)
2211 {
2212 static_assert(
2213 (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix> ||
2214 std::is_same_v<MatrixType,
2215 SparseMatrix<typename MatrixType::value_type>>));
2216
2217 Assert((dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_dg) != nullptr),
2218 ExcNotImplemented(
2219 "Implemented only for simplex meshes for the time being."));
2220
2221 Assert(dof_handler.get_triangulation().all_reference_cells_are_simplex(),
2222 ExcNotImplemented());
2223
2224 const double penalty_constant = .5 * fe_dg.degree * (fe_dg.degree + 1);
2225 AffineConstraints<typename MatrixType::value_type> constraints;
2226 constraints.close();
2227
2228 const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
2229 const IndexSet locally_relevant_dofs =
2230 DoFTools::extract_locally_relevant_dofs(dof_handler);
2231
2232 DynamicSparsityPattern dsp(locally_relevant_dofs);
2233 DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
2234 SparsityTools::distribute_sparsity_pattern(
2235 dsp,
2236 dof_handler.locally_owned_dofs(),
2237 dof_handler.get_mpi_communicator(),
2238 locally_relevant_dofs);
2239
2240 system_matrix.reinit(locally_owned_dofs,
2241 locally_owned_dofs,
2242 dsp,
2243 dof_handler.get_mpi_communicator());
2244
2245 system_rhs.reinit(locally_owned_dofs, dof_handler.get_mpi_communicator());
2246
2247 const unsigned int quadrature_degree = fe_dg.degree + 1;
2248 FEFaceValues<dim> fe_faces0(mapping,
2249 fe_dg,
2250 QGaussSimplex<dim - 1>(quadrature_degree),
2251 update_values | update_JxW_values |
2252 update_gradients | update_quadrature_points |
2253 update_normal_vectors);
2254
2255
2256 FEValues<dim> fe_values(mapping,
2257 fe_dg,
2258 QGaussSimplex<dim>(quadrature_degree),
2259 update_values | update_JxW_values |
2260 update_gradients | update_quadrature_points);
2261
2262 FEFaceValues<dim> fe_faces1(mapping,
2263 fe_dg,
2264 QGaussSimplex<dim - 1>(quadrature_degree),
2265 update_values | update_JxW_values |
2266 update_gradients | update_quadrature_points |
2267 update_normal_vectors);
2268 const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell();
2269
2270 FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
2271 Vector<double> cell_rhs(dofs_per_cell);
2272
2273 FullMatrix<double> M11(dofs_per_cell, dofs_per_cell);
2274 FullMatrix<double> M12(dofs_per_cell, dofs_per_cell);
2275 FullMatrix<double> M21(dofs_per_cell, dofs_per_cell);
2276 FullMatrix<double> M22(dofs_per_cell, dofs_per_cell);
2277
2278 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2279
2280 // Loop over standard deal.II cells
2281 for (const auto &cell : dof_handler.active_cell_iterators())
2282 {
2283 if (cell->is_locally_owned())
2284 {
2285 cell_matrix = 0.;
2286 cell_rhs = 0.;
2287
2288 fe_values.reinit(cell);
2289
2290 // const auto &q_points = fe_values.get_quadrature_points();
2291 // const unsigned int n_qpoints = q_points.size();
2292
2293 for (unsigned int q_index : fe_values.quadrature_point_indices())
2294 {
2295 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2296 {
2297 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2298 {
2299 cell_matrix(i, j) += fe_values.shape_grad(i, q_index) *
2300 fe_values.shape_grad(j, q_index) *
2301 fe_values.JxW(q_index);
2302 }
2303 cell_rhs(i) +=
2304 fe_values.shape_value(i, q_index) * 1. *
2305 fe_values.JxW(q_index); // TODO: pass functional
2306 }
2307 }
2308
2309 // distribute volumetric DoFs
2310 cell->get_dof_indices(local_dof_indices);
2311 double hf = 0.;
2312 for (const auto f : cell->face_indices())
2313 {
2314 const double extent1 =
2315 cell->measure() / cell->face(f)->measure();
2316
2317 if (cell->face(f)->at_boundary())
2318 {
2319 hf = (1. / extent1 + 1. / extent1);
2320 fe_faces0.reinit(cell, f);
2321
2322 const auto &normals = fe_faces0.get_normal_vectors();
2323 for (unsigned int q_index :
2324 fe_faces0.quadrature_point_indices())
2325 {
2326 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2327 {
2328 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2329 {
2330 cell_matrix(i, j) +=
2331 (-fe_faces0.shape_value(i, q_index) *
2332 fe_faces0.shape_grad(j, q_index) *
2333 normals[q_index] -
2334 fe_faces0.shape_grad(i, q_index) *
2335 normals[q_index] *
2336 fe_faces0.shape_value(j, q_index) +
2337 (penalty_constant * hf) *
2338 fe_faces0.shape_value(i, q_index) *
2339 fe_faces0.shape_value(j, q_index)) *
2340 fe_faces0.JxW(q_index);
2341 }
2342 cell_rhs(i) +=
2343 0.; // TODO: add bdary conditions functional
2344 }
2345 }
2346 }
2347 else
2348 {
2349 const auto &neigh_cell = cell->neighbor(f);
2350 if (cell->global_active_cell_index() <
2351 neigh_cell->global_active_cell_index())
2352 {
2353 const double extent2 =
2354 neigh_cell->measure() /
2355 neigh_cell->face(cell->neighbor_of_neighbor(f))
2356 ->measure();
2357 hf = (1. / extent1 + 1. / extent2);
2358 fe_faces0.reinit(cell, f);
2359 fe_faces1.reinit(neigh_cell,
2360 cell->neighbor_of_neighbor(f));
2361
2362 std::vector<types::global_dof_index>
2363 local_dof_indices_neighbor(dofs_per_cell);
2364
2365 M11 = 0.;
2366 M12 = 0.;
2367 M21 = 0.;
2368 M22 = 0.;
2369
2370 const auto &normals = fe_faces0.get_normal_vectors();
2371 // M11
2372 for (unsigned int q_index :
2373 fe_faces0.quadrature_point_indices())
2374 {
2375 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2376 {
2377 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2378 {
2379 M11(i, j) +=
2380 (-0.5 * fe_faces0.shape_grad(i, q_index) *
2381 normals[q_index] *
2382 fe_faces0.shape_value(j, q_index) -
2383 0.5 * fe_faces0.shape_grad(j, q_index) *
2384 normals[q_index] *
2385 fe_faces0.shape_value(i, q_index) +
2386 (penalty_constant * hf) *
2387 fe_faces0.shape_value(i, q_index) *
2388 fe_faces0.shape_value(j, q_index)) *
2389 fe_faces0.JxW(q_index);
2390
2391 M12(i, j) +=
2392 (0.5 * fe_faces0.shape_grad(i, q_index) *
2393 normals[q_index] *
2394 fe_faces1.shape_value(j, q_index) -
2395 0.5 * fe_faces1.shape_grad(j, q_index) *
2396 normals[q_index] *
2397 fe_faces0.shape_value(i, q_index) -
2398 (penalty_constant * hf) *
2399 fe_faces0.shape_value(i, q_index) *
2400 fe_faces1.shape_value(j, q_index)) *
2401 fe_faces1.JxW(q_index);
2402
2403 // A10
2404 M21(i, j) +=
2405 (-0.5 * fe_faces1.shape_grad(i, q_index) *
2406 normals[q_index] *
2407 fe_faces0.shape_value(j, q_index) +
2408 0.5 * fe_faces0.shape_grad(j, q_index) *
2409 normals[q_index] *
2410 fe_faces1.shape_value(i, q_index) -
2411 (penalty_constant * hf) *
2412 fe_faces1.shape_value(i, q_index) *
2413 fe_faces0.shape_value(j, q_index)) *
2414 fe_faces1.JxW(q_index);
2415
2416 // A11
2417 M22(i, j) +=
2418 (0.5 * fe_faces1.shape_grad(i, q_index) *
2419 normals[q_index] *
2420 fe_faces1.shape_value(j, q_index) +
2421 0.5 * fe_faces1.shape_grad(j, q_index) *
2422 normals[q_index] *
2423 fe_faces1.shape_value(i, q_index) +
2424 (penalty_constant * hf) *
2425 fe_faces1.shape_value(i, q_index) *
2426 fe_faces1.shape_value(j, q_index)) *
2427 fe_faces1.JxW(q_index);
2428 }
2429 }
2430 }
2431
2432 // distribute DoFs accordingly
2433
2434 neigh_cell->get_dof_indices(local_dof_indices_neighbor);
2435
2436 constraints.distribute_local_to_global(
2437 M11, local_dof_indices, system_matrix);
2438 constraints.distribute_local_to_global(
2439 M12,
2440 local_dof_indices,
2441 local_dof_indices_neighbor,
2442 system_matrix);
2443 constraints.distribute_local_to_global(
2444 M21,
2445 local_dof_indices_neighbor,
2446 local_dof_indices,
2447 system_matrix);
2448 constraints.distribute_local_to_global(
2449 M22, local_dof_indices_neighbor, system_matrix);
2450
2451 } // check idx neighbors
2452 } // over faces
2453 }
2454 constraints.distribute_local_to_global(cell_matrix,
2455 cell_rhs,
2456 local_dof_indices,
2457 system_matrix,
2458 system_rhs);
2459 }
2460 }
2461 system_matrix.compress(VectorOperation::add);
2462 system_rhs.compress(VectorOperation::add);
2463 }
2464
2465
2466} // namespace dealii::PolyUtils
2467
2468#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
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_Nothing
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::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:1099
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:996
AgglomerationHandler::get_dof_handler
const DoFHandler< dim, spacedim > & get_dof_handler() const
Definition agglomeration_handler.h:1038
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:1048
AgglomerationHandler::reinit
const FEValues< dim, spacedim > & reinit(const AgglomerationIterator< dim, spacedim > &polytope) const
Definition agglomeration_handler.cc:731
AgglomerationHandler::recv_jxws
std::map< types::subdomain_id, std::map< std::pair< CellId, unsigned int >, std::vector< double > > > recv_jxws
Definition agglomeration_handler.h:533
AgglomerationHandler::output_dh
DoFHandler< dim, spacedim > output_dh
Definition agglomeration_handler.h:506
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:808
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:548
AgglomerationHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
Definition agglomeration_handler.h:979
AgglomerationHandler::get_mapping
const Mapping< dim > & get_mapping() const
Definition agglomeration_handler.h:961
AgglomerationHandler::n_agglomerates
unsigned int n_agglomerates() const
Definition agglomeration_handler.h:1167
AgglomerationHandler::get_local_bboxes
const std::vector< BoundingBox< dim > > & get_local_bboxes() const
Definition agglomeration_handler.h:987
AgglomerationHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe() const
Definition agglomeration_handler.h:952
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:543
AgglomerationHandler::is_master_cell
bool is_master_cell(const CellIterator &cell) const
Definition agglomeration_handler.h:1059
AgglomerationHandler::polytope_iterators
IteratorRange< agglomeration_iterator > polytope_iterators() const
Definition agglomeration_handler.h:1243
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:501
AgglomerationHandler::get_fe_collection
const hp::FECollection< dim, spacedim > & get_fe_collection() const
Definition agglomeration_handler.h:1262
dealii::CellsAgglomerator::extract_agglomerates
const std::vector< std::vector< typename Triangulation< dim >::active_cell_iterator > > & extract_agglomerates()
Definition agglomerator.h:392
hp::FECollection::size
unsigned int size() const
hp::FECollection::push_back
void push_back(const FiniteElement< dim, spacedim > &new_fe)
conditional_ostream.h
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:1471
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:1872
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:1649
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:1763
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:2002
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:1932
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:2206
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
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
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
pack_rtree
RTree< typename LeafTypeIterator::value_type, IndexType, IndexableGetter > pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
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