Reference documentation for deal.II version GIT relicensing-480-geae235577e 2024-04-25 01:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
quadrature_generator.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
17#include "deal.II/fe/fe_q.h"
19
21
31#include <deal.II/lac/vector.h>
32
34
36
37#include <boost/math/special_functions/relative_difference.hpp>
38#include <boost/math/special_functions/sign.hpp>
39#include <boost/math/tools/roots.hpp>
40
41#include <algorithm>
42#include <vector>
43
45
46namespace NonMatching
47{
48 namespace internal
49 {
50 namespace QuadratureGeneratorImplementation
51 {
52 template <int dim>
53 void
55 const double weight,
56 const Quadrature<1> &quadrature1D,
57 const double start,
58 const double end,
59 const unsigned int component_in_dim,
60 ExtendableQuadrature<dim> &quadrature)
61 {
62 Assert(start < end,
63 ExcMessage("Interval start must be less than interval end."));
64
65 const double length = end - start;
66 for (unsigned int j = 0; j < quadrature1D.size(); ++j)
67 {
68 const double x = start + length * quadrature1D.point(j)[0];
70 point, component_in_dim, x),
71 length * weight * quadrature1D.weight(j));
72 }
73 }
74
75
76
81 template <int dim>
82 void
84 const Quadrature<1> &quadrature1D,
85 const double start,
86 const double end,
87 const unsigned int component_in_dim,
88 ExtendableQuadrature<dim> &quadrature)
89 {
90 for (unsigned int j = 0; j < lower.size(); ++j)
91 {
93 lower.weight(j),
94 quadrature1D,
95 start,
96 end,
97 component_in_dim,
98 quadrature);
99 }
100 }
101
102
103
104 template <int dim>
107 const std::vector<std::reference_wrapper<const Function<dim>>>
108 &functions,
109 const Point<dim> &point)
110 {
111 Assert(functions.size() > 0,
113 "The incoming vector must contain at least one function."));
114
115 const int sign_of_first =
116 boost::math::sign(functions[0].get().value(point));
117
118 if (sign_of_first == 0)
120
121 for (unsigned int j = 1; j < functions.size(); ++j)
122 {
123 const int sign = boost::math::sign(functions[j].get().value(point));
124
125 if (sign != sign_of_first)
127 }
128 // If we got here all functions have the same sign.
129 if (sign_of_first < 0)
131 else
133 }
134
135
136
149 template <int dim>
150 void
152 const BoundingBox<dim> &box,
153 std::pair<double, double> &value_bounds)
154 {
155 const ReferenceCell &cube = ReferenceCells::get_hypercube<dim>();
156 for (unsigned int i = 0; i < cube.n_vertices(); ++i)
157 {
158 const double vertex_value = function.value(box.vertex(i));
159
160 value_bounds.first = std::min(value_bounds.first, vertex_value);
161 value_bounds.second = std::max(value_bounds.second, vertex_value);
162 }
163 }
164
165
166
177 template <int dim>
178 void
180 const std::vector<std::reference_wrapper<const Function<dim>>>
181 &functions,
182 const BoundingBox<dim> &box,
183 std::vector<FunctionBounds<dim>> &all_function_bounds)
184 {
185 all_function_bounds.clear();
186 all_function_bounds.reserve(functions.size());
187 for (const Function<dim> &function : functions)
188 {
189 FunctionBounds<dim> bounds;
190 FunctionTools::taylor_estimate_function_bounds<dim>(
191 function, box, bounds.value, bounds.gradient);
192 take_min_max_at_vertices(function, box, bounds.value);
193
194 all_function_bounds.push_back(bounds);
195 }
196 }
197
198
200 template <int dim>
201 std::pair<double, double>
202 find_extreme_values(const std::vector<FunctionBounds<dim>> &bounds)
203 {
204 Assert(bounds.size() > 0, ExcMessage("The incoming vector is empty."));
205
206 std::pair<double, double> extremes = bounds[0].value;
207 for (unsigned int i = 1; i < bounds.size(); ++i)
208 {
209 extremes.first = std::min(extremes.first, bounds[i].value.first);
210 extremes.second = std::max(extremes.second, bounds[i].value.second);
211 }
212
213 return extremes;
214 }
215
216
217
222 inline bool
223 is_indefinite(const std::pair<double, double> &function_bounds)
224 {
225 if (function_bounds.first > 0)
226 return false;
227 if (function_bounds.second < 0)
228 return false;
229 return true;
230 }
231
232
233
254 inline double
255 lower_bound_on_abs(const std::pair<double, double> &function_bounds)
256 {
257 Assert(function_bounds.first <= function_bounds.second,
258 ExcMessage("Function bounds reversed, max < min."));
260 return 0.5 * (std::abs(function_bounds.second + function_bounds.first) -
261 (function_bounds.second - function_bounds.first));
262 }
263
264
265
271
272
273
274 template <int dim>
275 std::optional<HeightDirectionData>
277 const std::vector<FunctionBounds<dim>> &all_function_bounds)
278 {
279 // Minimum (taken over the indefinite functions) on the lower bound on
280 // each component of the gradient.
281 std::optional<std::array<double, dim>> min_lower_abs_grad;
282
283 for (const FunctionBounds<dim> &bounds : all_function_bounds)
284 {
285 if (is_indefinite(bounds.value))
286 {
287 // For the first indefinite function we find, we write the lower
288 // bounds on each gradient component to min_lower_abs_grad.
289 if (!min_lower_abs_grad)
290 {
291 min_lower_abs_grad.emplace();
292 for (unsigned int d = 0; d < dim; ++d)
293 {
294 (*min_lower_abs_grad)[d] =
295 lower_bound_on_abs(bounds.gradient[d]);
296 }
297 }
298 else
299 {
300 for (unsigned int d = 0; d < dim; ++d)
301 {
302 (*min_lower_abs_grad)[d] =
303 std::min((*min_lower_abs_grad)[d],
304 lower_bound_on_abs(bounds.gradient[d]));
305 }
306 }
307 }
308 }
309
310 if (min_lower_abs_grad)
311 {
312 const auto max_element =
313 std::max_element(min_lower_abs_grad->begin(),
314 min_lower_abs_grad->end());
315
317 data.direction =
318 std::distance(min_lower_abs_grad->begin(), max_element);
319 data.min_abs_dfdx = *max_element;
320
321 return data;
322 }
323
324 return std::optional<HeightDirectionData>();
325 }
326
327
328
334 template <int dim>
335 inline bool
337 const std::vector<FunctionBounds<dim>> &all_function_bounds)
338 {
339 if (all_function_bounds.size() != 2)
340 return false;
341 else
342 {
343 const FunctionBounds<dim> &bounds0 = all_function_bounds.at(0);
344 const FunctionBounds<dim> &bounds1 = all_function_bounds.at(1);
345
346 if (bounds0.value.first > 0 && bounds1.value.second < 0)
347 return true;
348 if (bounds1.value.first > 0 && bounds0.value.second < 0)
349 return true;
350 return false;
351 }
352 }
353
354
355
363 template <int dim>
364 void
366 const BoundingBox<dim> &box,
367 ExtendableQuadrature<dim> &quadrature)
368 {
369 for (unsigned int i = 0; i < unit_quadrature.size(); ++i)
370 {
371 const Point<dim> point = box.unit_to_real(unit_quadrature.point(i));
372 const double weight = unit_quadrature.weight(i) * box.volume();
373
374 quadrature.push_back(point, weight);
375 }
376 }
377
378
379
390 template <int dim>
391 void
393 const std::vector<std::reference_wrapper<const Function<dim>>>
394 &functions,
395 const BoundingBox<dim> &box,
396 const unsigned int direction,
397 std::vector<Functions::CoordinateRestriction<dim - 1>> &restrictions)
398 {
399 AssertIndexRange(direction, dim);
400
401 restrictions.clear();
402 restrictions.reserve(2 * functions.size());
403
404 const double bottom = box.lower_bound(direction);
405 const double top = box.upper_bound(direction);
406
407 for (const auto &function : functions)
408 {
409 restrictions.push_back(Functions::CoordinateRestriction<dim - 1>(
410 function, direction, bottom));
411 restrictions.push_back(Functions::CoordinateRestriction<dim - 1>(
412 function, direction, top));
413 }
414 }
415
416
417
426 template <int dim>
427 void
429 const std::vector<std::reference_wrapper<const Function<dim>>>
430 &functions,
431 const Point<dim - 1> &point,
432 const unsigned int open_direction,
433 std::vector<Functions::PointRestriction<dim - 1>> &restrictions)
434 {
435 AssertIndexRange(open_direction, dim);
436
437 restrictions.clear();
438 restrictions.reserve(functions.size());
439 for (const auto &function : functions)
440 {
441 restrictions.push_back(Functions::PointRestriction<dim - 1>(
442 function, open_direction, point));
443 }
444 }
445
446
447
461 template <int dim>
462 void
464 const Quadrature<1> &quadrature1D,
465 const BoundingBox<1> &interval,
466 const std::vector<double> &roots,
467 const Point<dim - 1> &point,
468 const double weight,
469 const unsigned int height_function_direction,
470 const std::vector<std::reference_wrapper<const Function<1>>>
471 &level_sets,
472 const AdditionalQGeneratorData &additional_data,
473 QPartitioning<dim> &q_partitioning)
474 {
475 // Make this int to avoid a warning signed/unsigned comparison.
476 const int n_roots = roots.size();
477
478 // The number of intervals are roots.size() + 1
479 for (int i = -1; i < n_roots; ++i)
480 {
481 // Start and end point of the subinterval.
482 const double start = i < 0 ? interval.lower_bound(0) : roots[i];
483 const double end =
484 i + 1 < n_roots ? roots[i + 1] : interval.upper_bound(0);
485
486 const double length = end - start;
487 // It might be that the end points of the subinterval are roots.
488 // If this is the case then the subinterval has length zero.
489 // Don't distribute points on the subinterval if it is shorter than
490 // some tolerance.
491 if (length > additional_data.min_distance_between_roots)
492 {
493 // All points on the interval belong to the same region in
494 // the QPartitioning. Determine the quadrature we should add
495 // the points to.
496 const Point<1> center(start + 0.5 * length);
497 const Definiteness definiteness =
498 pointwise_definiteness(level_sets, center);
499 ExtendableQuadrature<dim> &target_quadrature =
500 q_partitioning.quadrature_by_definiteness(definiteness);
501
503 weight,
504 quadrature1D,
505 start,
506 end,
507 height_function_direction,
508 target_quadrature);
509 }
510 }
511 }
512
513
514
516 const double tolerance,
517 const unsigned int max_recursion_depth,
518 const unsigned int max_iterations)
519 : tolerance(tolerance)
520 , max_recursion_depth(max_recursion_depth)
521 , max_iterations(max_iterations)
522 {}
523
524
525
529
530
531
532 void
534 const std::vector<std::reference_wrapper<const Function<1>>> &functions,
535 const BoundingBox<1> &interval,
536 std::vector<double> &roots)
537 {
538 for (const Function<1> &function : functions)
539 {
540 const unsigned int recursion_depth = 0;
541 find_roots(function, interval, recursion_depth, roots);
542 }
543 // Sort and make sure no roots are duplicated
544 std::sort(roots.begin(), roots.end());
545
546 const auto roots_are_equal = [this](const double &a, const double &b) {
547 return std::abs(a - b) < additional_data.tolerance;
548 };
549 roots.erase(unique(roots.begin(), roots.end(), roots_are_equal),
550 roots.end());
551 }
552
553
554
555 void
557 const BoundingBox<1> &interval,
558 const unsigned int recursion_depth,
559 std::vector<double> &roots)
560 {
561 // Compute function values at end points.
562 const double left_value = function.value(interval.vertex(0));
563 const double right_value = function.value(interval.vertex(1));
564
565 // If we have a sign change we solve for the root.
566 if (boost::math::sign(left_value) != boost::math::sign(right_value))
567 {
568 const auto lambda = [&function](const double x) {
569 return function.value(Point<1>(x));
570 };
571
572 const auto stopping_criteria = [this](const double &a,
573 const double &b) {
574 return std::abs(a - b) < additional_data.tolerance;
575 };
576
577 boost::uintmax_t iterations = additional_data.max_iterations;
578
579 const std::pair<double, double> root_bracket =
580 boost::math::tools::toms748_solve(lambda,
581 interval.lower_bound(0),
582 interval.upper_bound(0),
583 left_value,
584 right_value,
585 stopping_criteria,
586 iterations);
587
588 const double root = .5 * (root_bracket.first + root_bracket.second);
589 roots.push_back(root);
590 }
591 else
592 {
593 // Compute bounds on the incoming function to check if there are
594 // roots. If the function is positive or negative on the whole
595 // interval we do nothing.
596 std::pair<double, double> value_bounds;
597 std::array<std::pair<double, double>, 1> gradient_bounds;
598 FunctionTools::taylor_estimate_function_bounds<1>(function,
599 interval,
600 value_bounds,
601 gradient_bounds);
602
603 // Since we already know the function values at the interval ends we
604 // might as well check these for min/max too.
605 const double function_min =
606 std::min(std::min(left_value, right_value), value_bounds.first);
607
608 // If the functions is positive there are no roots.
609 if (function_min > 0)
610 return;
611
612 const double function_max =
613 std::max(std::max(left_value, right_value), value_bounds.second);
614
615 // If the function is negative there are no roots.
616 if (function_max < 0)
617 return;
618
619 // If we can't say that the function is strictly positive/negative
620 // we split the interval in half. We can't split forever, so if we
621 // have reached the max recursion, we stop looking for roots.
622 if (recursion_depth < additional_data.max_recursion_depth)
623 {
624 find_roots(function,
625 interval.child(0),
626 recursion_depth + 1,
627 roots);
628 find_roots(function,
629 interval.child(1),
630 recursion_depth + 1,
631 roots);
632 }
633 }
634 }
635
636
637
638 template <int dim>
640 const Quadrature<dim> &quadrature)
641 : Quadrature<dim>(quadrature)
642 {}
643
644
645
646 template <int dim>
647 inline void
649 {
650 this->quadrature_points.clear();
651 this->weights.clear();
652 }
653
654
655
656 template <int dim>
657 void
659 const double weight)
660 {
661 this->quadrature_points.push_back(point);
662 this->weights.push_back(weight);
663 }
664
665
666
667 template <int dim>
670 const Definiteness definiteness)
671 {
672 switch (definiteness)
673 {
675 return negative;
677 return positive;
678 default:
679 return indefinite;
680 }
681 }
682
683
684
685 template <int dim>
686 void
688 {
689 positive.clear();
690 negative.clear();
691 indefinite.clear();
692 surface.clear();
693 }
694
695
696
704 template <int dim>
707 const unsigned int direction,
708 const BoundingBox<dim> &box,
709 const Function<dim> &level_set)
710 {
711 const Point<dim> bottom_point =
713 direction,
714 box.lower_bound(direction));
715 const double bottom_value = level_set.value(bottom_point);
716
717 const Point<dim> top_point =
719 direction,
720 box.upper_bound(direction));
721 const double top_value = level_set.value(top_point);
722
723 // The end point closest to the zero-contour is the one with smallest
724 // absolute value.
725 return std::abs(bottom_value) < std::abs(top_value) ? bottom_point :
726 top_point;
727 }
728
729
730
731 template <int dim, int spacedim>
733 const hp::QCollection<1> &q_collection1D,
734 const AdditionalQGeneratorData &additional_data)
735 : q_collection1D(&q_collection1D)
736 , additional_data(additional_data)
737 , root_finder(
738 RootFinder::AdditionalData(additional_data.root_finder_tolerance,
739 additional_data.max_root_finder_splits))
740 {
741 q_index = 0;
742 }
743
744
745
746 template <int dim, int spacedim>
747 void
749 const std::vector<std::reference_wrapper<const Function<dim>>>
750 &level_sets,
751 const BoundingBox<dim> &box,
752 const Quadrature<dim - 1> &low_dim_quadrature,
753 const unsigned int height_function_direction,
754 QPartitioning<dim> &q_partitioning)
755 {
756 const Quadrature<1> &quadrature1D = (*q_collection1D)[q_index];
757
758 for (unsigned int q = 0; q < low_dim_quadrature.size(); ++q)
759 {
760 const Point<dim - 1> &point = low_dim_quadrature.point(q);
761 const double weight = low_dim_quadrature.weight(q);
762 restrict_to_point(level_sets,
763 point,
764 height_function_direction,
765 point_restrictions);
766
767 // We need a vector of references to do the recursive call.
768 const std::vector<std::reference_wrapper<const Function<1>>>
769 restrictions(point_restrictions.begin(),
770 point_restrictions.end());
771
772 const BoundingBox<1> bounds_in_direction =
773 box.bounds(height_function_direction);
774
775 roots.clear();
776 root_finder.find_roots(restrictions, bounds_in_direction, roots);
777
779 bounds_in_direction,
780 roots,
781 point,
782 weight,
783 height_function_direction,
784 restrictions,
785 additional_data,
786 q_partitioning);
787
788 if (dim == spacedim)
789 create_surface_point(point,
790 weight,
791 level_sets,
792 box,
793 height_function_direction,
794 q_partitioning.surface);
795 }
796
797 point_restrictions.clear();
798 }
799
800
801
802 template <int dim, int spacedim>
803 void
805 const Point<dim - 1> &point,
806 const double weight,
807 const std::vector<std::reference_wrapper<const Function<dim>>>
808 &level_sets,
809 const BoundingBox<dim> &box,
810 const unsigned int height_function_direction,
811 ImmersedSurfaceQuadrature<dim> &surface_quadrature)
812 {
813 AssertIndexRange(roots.size(), 2);
814 Assert(level_sets.size() == 1, ExcInternalError());
815
816
817 const Function<dim> &level_set = level_sets.at(0);
818
819 Point<dim> surface_point;
820 if (roots.size() == 1)
821 {
823 point, height_function_direction, roots[0]);
824 }
825 else
826 {
827 // If we got here, we have missed roots in the lower dimensional
828 // algorithm. This is a rare event but can happen if the
829 // zero-contour has a high curvature. The problem is that the
830 // incoming point has been incorrectly added to the indefinite
831 // quadrature in QPartitioning<dim-1>. Since we missed a root on
832 // this box, we will likely miss it on the neighboring box too. If
833 // this happens, the point will NOT be in the indefinite quadrature
834 // on the neighbor. The best thing we can do is to compute the
835 // surface point by projecting the lower dimensional point to the
836 // face closest to the zero-contour. We should add a surface point
837 // because the neighbor will not.
839 point, height_function_direction, box, level_set);
840 }
841
842 const Tensor<1, dim> gradient = level_set.gradient(surface_point);
843 Tensor<1, dim> normal = gradient;
844 normal *= 1. / normal.norm();
845
846 // Note that gradient[height_function_direction] is non-zero
847 // because of the implicit function theorem.
848 const double surface_weight =
849 weight * gradient.norm() /
850 std::abs(gradient[height_function_direction]);
851 surface_quadrature.push_back(surface_point, surface_weight, normal);
852 }
853
854
855
856 template <int dim, int spacedim>
857 void
859 unsigned int q_index)
860 {
861 AssertIndexRange(q_index, q_collection1D->size());
862 this->q_index = q_index;
863 }
864
865
866
867 template <int dim, int spacedim>
869 const hp::QCollection<1> &q_collection1D,
870 const AdditionalQGeneratorData &additional_data)
871 : additional_data(additional_data)
872 , q_collection1D(&q_collection1D)
873 {
874 q_index = 0;
875 }
876
877
878
879 template <int dim, int spacedim>
881 const hp::QCollection<1> &q_collection1D,
882 const AdditionalQGeneratorData &additional_data)
883 : QGeneratorBase<dim, spacedim>(q_collection1D, additional_data)
884 , low_dim_algorithm(q_collection1D, additional_data)
885 , up_through_dimension_creator(q_collection1D, additional_data)
886 {
887 for (const auto &quadrature : q_collection1D)
888 tensor_products.push_back(quadrature);
889 }
890
891
892
893 template <int dim, int spacedim>
894 void
896 {
897 q_partitioning.clear();
898 }
899
900
901
902 template <int dim, int spacedim>
903 const QPartitioning<dim> &
905 {
906 return q_partitioning;
907 }
908
909
910
911 template <int dim, int spacedim>
912 void
914 const std::vector<std::reference_wrapper<const Function<dim>>>
915 &level_sets,
916 const BoundingBox<dim> &box,
917 const unsigned int n_box_splits)
918 {
919 std::vector<FunctionBounds<dim>> all_function_bounds;
920 estimate_function_bounds(level_sets, box, all_function_bounds);
921
922 const std::pair<double, double> extreme_values =
923 find_extreme_values(all_function_bounds);
924
925 if (extreme_values.first > this->additional_data.limit_to_be_definite)
926 {
927 map_quadrature_to_box(tensor_products[this->q_index],
928 box,
929 this->q_partitioning.positive);
930 }
931 else if (extreme_values.second <
932 -(this->additional_data.limit_to_be_definite))
933 {
934 map_quadrature_to_box(tensor_products[this->q_index],
935 box,
936 this->q_partitioning.negative);
937 }
938 else if (one_positive_one_negative_definite(all_function_bounds))
939 {
940 map_quadrature_to_box(tensor_products[this->q_index],
941 box,
942 this->q_partitioning.indefinite);
943 }
944 else
945 {
946 const std::optional<HeightDirectionData> data =
947 find_best_height_direction(all_function_bounds);
948
949 // Check larger than a constant to avoid that min_abs_dfdx is only
950 // larger by 0 by floating point precision.
951 if (data && data->min_abs_dfdx >
952 this->additional_data.lower_bound_implicit_function)
953 {
954 create_low_dim_quadratures(data->direction,
955 level_sets,
956 box,
957 n_box_splits);
958 create_high_dim_quadratures(data->direction, level_sets, box);
959 }
960 else if (n_box_splits < this->additional_data.max_box_splits)
961 {
962 split_box_and_recurse(level_sets, box, data, n_box_splits);
963 }
964 else
965 {
966 // We can't split the box recursively forever. Use the midpoint
967 // method as a last resort.
968 use_midpoint_method(level_sets, box);
969 }
970 }
971 }
972
973
974
980 template <int dim>
981 std::optional<unsigned int>
983 {
984 // Get the side lengths for each direction and sort them.
985 std::array<std::pair<double, unsigned int>, dim> side_lengths;
986 for (unsigned int i = 0; i < dim; ++i)
987 {
988 side_lengths[i].first = box.side_length(i);
989 side_lengths[i].second = i;
990 }
991 // Sort is lexicographic, so this sorts based on side length first.
992 std::sort(side_lengths.begin(), side_lengths.end());
993
994 // Check if the two largest side lengths have the same length. This
995 // function isn't called in 1d, so the (dim - 2)-element exists.
996 if (boost::math::epsilon_difference(side_lengths[dim - 1].first,
997 side_lengths[dim - 2].first) < 100)
998 return std::optional<unsigned int>();
999
1000 return side_lengths.back().second;
1001 }
1002
1003
1004
1016 template <int dim>
1017 unsigned int
1019 const BoundingBox<dim> &box,
1020 const std::optional<HeightDirectionData> &height_direction_data)
1021 {
1022 const std::optional<unsigned int> direction =
1024
1025 if (direction)
1026 return *direction;
1027
1028 // This direction is closest to being a height direction, so
1029 // we split in this direction.
1030 if (height_direction_data)
1031 return height_direction_data->direction;
1032
1033 // We have to choose some direction, we might as well take 0.
1034 return 0;
1035 }
1037
1038
1043 template <int dim>
1044 inline BoundingBox<dim>
1045 left_half(const BoundingBox<dim> &box, const unsigned int direction)
1046 {
1047 AssertIndexRange(direction, dim);
1048
1049 // Move the upper corner half a side-length to the left.
1050 std::pair<Point<dim>, Point<dim>> corners = box.get_boundary_points();
1051 corners.second[direction] -= .5 * box.side_length(direction);
1052
1053 return BoundingBox<dim>(corners);
1054 }
1055
1056
1057
1062 template <int dim>
1063 inline BoundingBox<dim>
1064 right_half(const BoundingBox<dim> &box, const unsigned int direction)
1065 {
1066 AssertIndexRange(direction, dim);
1067
1068 // Move the lower corner half a side-length to the right.
1069 std::pair<Point<dim>, Point<dim>> corners = box.get_boundary_points();
1070 corners.first[direction] += .5 * box.side_length(direction);
1071
1072 return BoundingBox<dim>(corners);
1073 }
1074
1075
1076
1077 template <int dim, int spacedim>
1078 void
1080 const std::vector<std::reference_wrapper<const Function<dim>>>
1081 &level_sets,
1082 const BoundingBox<dim> &box,
1083 const std::optional<HeightDirectionData> &direction_data,
1084 const unsigned int n_box_splits)
1085 {
1086 if (this->additional_data.split_in_half)
1087 {
1088 const unsigned int direction =
1089 compute_split_direction(box, direction_data);
1090
1091 const BoundingBox<dim> left_box = left_half(box, direction);
1092 const BoundingBox<dim> right_box = right_half(box, direction);
1093
1094 generate(level_sets, left_box, n_box_splits + 1);
1095 generate(level_sets, right_box, n_box_splits + 1);
1096 }
1097 else
1098 {
1099 for (unsigned int i = 0;
1100 i < GeometryInfo<dim>::max_children_per_cell;
1101 ++i)
1102 {
1103 generate(level_sets, box.child(i), n_box_splits + 1);
1104 }
1105 }
1106 }
1107
1108
1109
1110 template <int dim, int spacedim>
1111 void
1113 const unsigned int height_function_direction,
1114 const std::vector<std::reference_wrapper<const Function<dim>>>
1115 &level_sets,
1116 const BoundingBox<dim> &box,
1117 const unsigned int n_box_splits)
1118 {
1119 std::vector<Functions::CoordinateRestriction<dim - 1>>
1120 face_restrictions;
1121 restrict_to_top_and_bottom(level_sets,
1122 box,
1123 height_function_direction,
1124 face_restrictions);
1125
1126 // We need a vector of references to do the recursive call.
1127 const std::vector<std::reference_wrapper<const Function<dim - 1>>>
1128 restrictions(face_restrictions.begin(), face_restrictions.end());
1129
1130 const BoundingBox<dim - 1> cross_section =
1131 box.cross_section(height_function_direction);
1132
1133 low_dim_algorithm.clear_quadratures();
1134 low_dim_algorithm.generate(restrictions, cross_section, n_box_splits);
1135 }
1136
1137
1138
1139 template <int dim, int spacedim>
1140 void
1142 const unsigned int height_function_direction,
1143 const std::vector<std::reference_wrapper<const Function<dim>>>
1144 &level_sets,
1145 const BoundingBox<dim> &box)
1146 {
1147 const QPartitioning<dim - 1> &low_dim_quadratures =
1148 low_dim_algorithm.get_quadratures();
1149
1150 const Quadrature<1> &quadrature1D =
1151 (*this->q_collection1D)[this->q_index];
1152
1153 add_tensor_product(low_dim_quadratures.negative,
1154 quadrature1D,
1155 box.lower_bound(height_function_direction),
1156 box.upper_bound(height_function_direction),
1157 height_function_direction,
1158 this->q_partitioning.negative);
1160 add_tensor_product(low_dim_quadratures.positive,
1161 quadrature1D,
1162 box.lower_bound(height_function_direction),
1163 box.upper_bound(height_function_direction),
1164 height_function_direction,
1165 this->q_partitioning.positive);
1166
1167 up_through_dimension_creator.generate(level_sets,
1168 box,
1169 low_dim_quadratures.indefinite,
1170 height_function_direction,
1171 this->q_partitioning);
1172 }
1174
1175
1176 template <int dim, int spacedim>
1177 void
1179 const std::vector<std::reference_wrapper<const Function<dim>>>
1180 &level_sets,
1181 const BoundingBox<dim> &box)
1182 {
1183 const Point<dim> center = box.center();
1184 const Definiteness definiteness =
1185 pointwise_definiteness(level_sets, center);
1186
1187 ExtendableQuadrature<dim> &quadrature =
1188 this->q_partitioning.quadrature_by_definiteness(definiteness);
1189
1190 quadrature.push_back(center, box.volume());
1191 }
1192
1194
1195 template <int dim, int spacedim>
1196 void
1198 {
1199 AssertIndexRange(q_index, this->q_collection1D->size());
1200
1201 this->q_index = q_index;
1202 low_dim_algorithm.set_1D_quadrature(q_index);
1203 up_through_dimension_creator.set_1D_quadrature(q_index);
1204 }
1205
1207
1208 template <int spacedim>
1210 const hp::QCollection<1> &q_collection1D,
1211 const AdditionalQGeneratorData &additional_data)
1212 : QGeneratorBase<1, spacedim>(q_collection1D, additional_data)
1213 , root_finder(
1214 RootFinder::AdditionalData(additional_data.root_finder_tolerance,
1215 additional_data.max_root_finder_splits))
1216 {
1217 Assert(q_collection1D.size() > 0,
1218 ExcMessage("Incoming quadrature collection is empty."));
1219 }
1220
1221
1222
1223 template <int spacedim>
1224 void
1226 const std::vector<std::reference_wrapper<const Function<1>>>
1227 &level_sets,
1228 const BoundingBox<1> &box,
1229 const unsigned int n_box_splits)
1230 {
1231 (void)n_box_splits;
1232
1233 roots.clear();
1234 root_finder.find_roots(level_sets, box, roots);
1236 const Quadrature<1> &quadrature1D =
1237 (*this->q_collection1D)[this->q_index];
1238
1240 box,
1241 roots,
1242 zero_dim_point,
1243 unit_weight,
1244 direction,
1245 level_sets,
1246 this->additional_data,
1247 this->q_partitioning);
1248
1249 if (spacedim == 1)
1250 this->create_surface_points(level_sets);
1251 }
1252
1253
1254
1255 template <int spacedim>
1256 void
1258 const std::vector<std::reference_wrapper<const Function<1>>>
1259 &level_sets)
1260 {
1261 Assert(level_sets.size() == 1, ExcInternalError());
1262
1263 for (const double root : roots)
1264 {
1265 // A surface integral in 1d is just a point evaluation,
1266 // so the weight is always 1.
1267 const double weight = 1;
1268 const Point<1> point(root);
1269
1270 Tensor<1, 1> normal = level_sets[0].get().gradient(point);
1271 const double gradient_norm = normal.norm();
1272 Assert(
1273 gradient_norm > 1e-11,
1274 ExcMessage(
1275 "The level set function has a gradient almost equal to 0."));
1276 normal *= 1. / gradient_norm;
1277
1278 this->q_partitioning.surface.push_back(point, weight, normal);
1279 }
1280 }
1281
1282
1283
1284 template <int spacedim>
1285 void
1287 {
1288 AssertIndexRange(q_index, this->q_collection1D->size());
1289 this->q_index = q_index;
1290 }
1291
1292
1293
1294 } // namespace QuadratureGeneratorImplementation
1295
1296
1297
1298 namespace DiscreteQuadratureGeneratorImplementation
1299 {
1302 "The set_active_cell function has to be called before calling this function.");
1303
1306 "The reference cell of the incoming cell must be a hypercube.");
1307
1308
1334 template <int dim, typename VectorType = Vector<double>>
1336 {
1337 public:
1346 const VectorType &dof_values);
1347
1351 void
1353 &cell) override;
1354
1358 void
1359 set_subcell(const std::vector<unsigned int> &mask,
1360 const BoundingBox<dim> &subcell_box) override;
1361
1365 bool
1366 is_fe_q_iso_q1() const override;
1367
1371 unsigned int
1372 n_subdivisions() const override;
1373
1381 double
1382 value(const Point<dim> &point,
1383 const unsigned int component = 0) const override;
1384
1393 gradient(const Point<dim> &point,
1394 const unsigned int component = 0) const override;
1395
1404 hessian(const Point<dim> &point,
1405 const unsigned int component = 0) const override;
1406
1407 private:
1411 bool
1412 cell_is_set() const;
1413
1418
1424
1430
1434 std::vector<types::global_dof_index> local_dof_indices;
1435
1440 std::vector<typename VectorType::value_type> local_dof_values;
1441
1446 std::vector<typename VectorType::value_type> local_dof_values_subcell;
1447
1452
1458
1464 std::vector<Polynomials::Polynomial<double>> poly;
1465
1469 std::vector<unsigned int> renumber;
1470
1475
1479 Lazy<std::unique_ptr<FE_Q<dim>>> fe_q1;
1480 };
1481
1482
1483
1484 template <int dim, typename VectorType>
1486 const DoFHandler<dim> &dof_handler,
1487 const VectorType &dof_values)
1488 : dof_handler(&dof_handler)
1489 , global_dof_values(&dof_values)
1490 , n_subdivisions_per_line(numbers::invalid_unsigned_int)
1491 {
1492 Assert(dof_handler.n_dofs() == dof_values.size(),
1493 ExcDimensionMismatch(dof_handler.n_dofs(), dof_values.size()));
1494 }
1495
1496
1497
1498 template <int dim, typename VectorType>
1499 void
1501 const typename Triangulation<dim>::active_cell_iterator &cell)
1502 {
1503 Assert(
1504 &cell->get_triangulation() == &dof_handler->get_triangulation(),
1505 ExcMessage(
1506 "The incoming cell must belong to the triangulation associated with "
1507 "the DoFHandler passed to the constructor."));
1508
1509 const auto dof_handler_cell =
1510 cell->as_dof_handler_iterator(*dof_handler);
1511
1512 // Save the element and the local dof values, since this is what we need
1513 // to evaluate the function.
1514
1515 // Check if we can use the fast path. In case we have a different
1516 // element from the one used before we need to set up the data
1517 // structures again.
1518 if (element != &dof_handler_cell->get_fe())
1519 {
1520 poly.clear();
1521 element = &dof_handler_cell->get_fe();
1522
1523 const FiniteElement<dim> *fe = &dof_handler_cell->get_fe();
1524
1525 if (const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
1526 dynamic_cast<const FE_Q_iso_Q1<dim> *>(
1527 &dof_handler_cell->get_fe()))
1528 {
1529 this->n_subdivisions_per_line = fe_q_iso_q1->get_degree();
1530
1531 fe = fe_q1
1532 .value_or_initialize(
1533 []() { return std::make_unique<FE_Q<dim>>(1); })
1534 .get();
1535 local_dof_values_subcell.resize(
1536 fe_q1.value()->n_dofs_per_cell());
1537 }
1538 else
1539 this->n_subdivisions_per_line = numbers::invalid_unsigned_int;
1540
1541 if (fe->n_base_elements() == 1 &&
1543 0))
1544 {
1546 shape_info;
1547
1548 shape_info.reinit(QMidpoint<1>(), *element, 0);
1549 renumber = shape_info.lexicographic_numbering;
1550 poly =
1552 fe->base_element(0));
1553
1554 polynomials_are_hat_functions =
1555 (poly.size() == 2 && poly[0].value(0.) == 1. &&
1556 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1557 poly[1].value(1.) == 1.);
1558 }
1559 }
1560 else
1561 element = &dof_handler_cell->get_fe();
1562
1563 local_dof_indices.resize(element->dofs_per_cell);
1564 dof_handler_cell->get_dof_indices(local_dof_indices);
1565
1566 local_dof_values.resize(element->dofs_per_cell);
1567
1568 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1569 local_dof_values[i] =
1571 *global_dof_values, local_dof_indices[i]);
1572 }
1573
1574
1575
1576 template <int dim, typename VectorType>
1577 void
1579 const std::vector<unsigned int> &mask,
1580 const BoundingBox<dim> &subcell_box)
1581 {
1582 for (unsigned int i = 0; i < local_dof_values_subcell.size(); ++i)
1583 local_dof_values_subcell[i] = local_dof_values[renumber[mask[i]]];
1584
1585 this->subcell_box = subcell_box;
1586 }
1587
1588
1589
1590 template <int dim, typename VectorType>
1591 bool
1593 {
1594 return n_subdivisions_per_line != numbers::invalid_unsigned_int;
1595 }
1596
1597
1598
1599 template <int dim, typename VectorType>
1600 unsigned int
1602 {
1603 return n_subdivisions_per_line;
1604 }
1605
1606
1607
1608 template <int dim, typename VectorType>
1609 bool
1611 {
1612 // If set cell hasn't been called the size of local_dof_values will be
1613 // zero.
1614 return local_dof_values.size() > 0;
1615 }
1616
1617
1618
1619 template <int dim, typename VectorType>
1620 double
1622 const Point<dim> &point,
1623 const unsigned int component) const
1624 {
1625 AssertIndexRange(component, this->n_components);
1626 Assert(cell_is_set(), ExcCellNotSet());
1627
1628 if (!poly.empty() && component == 0)
1629 {
1630 // TODO: this could be extended to a component that is not zero
1631 return this->is_fe_q_iso_q1() ?
1633 poly,
1634 make_array_view(local_dof_values_subcell),
1635 subcell_box.real_to_unit(point),
1636 polynomials_are_hat_functions) :
1638 poly,
1639 make_array_view(local_dof_values),
1640 point,
1641 polynomials_are_hat_functions,
1642 renumber);
1643 }
1644 else
1645 {
1646 double value = 0;
1647 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1648 value += local_dof_values[i] *
1649 element->shape_value_component(i, point, component);
1650
1651 return value;
1652 }
1653 }
1654
1655
1656
1657 template <int dim, typename VectorType>
1660 const Point<dim> &point,
1661 const unsigned int component) const
1662 {
1663 AssertIndexRange(component, this->n_components);
1664 Assert(cell_is_set(), ExcCellNotSet());
1665
1666 if (!poly.empty() && component == 0)
1667 {
1668 // TODO: this could be extended to a component that is not zero
1669 return (this->is_fe_q_iso_q1() ?
1670 ::internal::
1671 evaluate_tensor_product_value_and_gradient(
1672 poly,
1673 make_array_view(local_dof_values_subcell),
1674 subcell_box.real_to_unit(point),
1675 polynomials_are_hat_functions) :
1678 poly,
1679 make_array_view(local_dof_values),
1680 point,
1681 polynomials_are_hat_functions,
1682 renumber))
1683 .second;
1684 }
1685 else
1686 {
1687 Tensor<1, dim> gradient;
1688 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1689 gradient += local_dof_values[i] *
1690 element->shape_grad_component(i, point, component);
1691
1692 return gradient;
1693 }
1694 }
1695
1696
1697
1698 template <int dim, typename VectorType>
1701 const Point<dim> &point,
1702 const unsigned int component) const
1703 {
1704 AssertIndexRange(component, this->n_components);
1705 Assert(cell_is_set(), ExcCellNotSet());
1706
1707 if (!poly.empty() && component == 0)
1708 {
1709 // TODO: this could be extended to a component that is not zero
1710 return this->is_fe_q_iso_q1() ?
1712 poly,
1713 make_array_view(local_dof_values_subcell),
1714 subcell_box.real_to_unit(point)) :
1716 poly,
1717 make_array_view(local_dof_values),
1718 point,
1719 renumber);
1720 }
1721 else
1722 {
1723 Tensor<2, dim> hessian;
1724 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1725 hessian +=
1726 local_dof_values[i] *
1727 element->shape_grad_grad_component(i, point, component);
1728
1729 return symmetrize(hessian);
1730 }
1731 }
1732
1733
1734
1735 template <int dim>
1738 const std::array<unsigned int, dim> &subcell_indices,
1739 const unsigned int n_subdivisions)
1740 {
1741 std::pair<Point<dim>, Point<dim>> corners;
1742 for (unsigned int d = 0; d < dim; ++d)
1743 {
1744 AssertIndexRange(subcell_indices[d], n_subdivisions + 1);
1745
1746 const double split_box_side_length =
1747 (1. / n_subdivisions) * unit_box.side_length(d);
1748 corners.first[d] = subcell_indices[d] * split_box_side_length;
1749 corners.second[d] = corners.first[d] + split_box_side_length;
1750 }
1751
1752 return BoundingBox<dim>(corners);
1753 }
1754
1755
1756
1757 template <int dim>
1758 std::vector<unsigned int>
1760 const std::array<unsigned int, dim> &subcell_indices,
1761 unsigned int dofs_per_line)
1762 {
1763 // This code expands the tensor product space, with `mask` increasing
1764 // in size by a factor `2` for each iteration in `d`. Note how the
1765 // indices in the next line are appended at the end of the vector.
1766 std::vector<unsigned int> mask;
1767 mask.reserve(1 << dim);
1768 mask.push_back(0);
1769 unsigned int stride = 1;
1770 for (unsigned int d = 0; d < dim; ++d)
1771 {
1772 const unsigned int previous_mask_size = 1U << d;
1773 AssertDimension(mask.size(), previous_mask_size);
1774 for (unsigned int i = 0; i < previous_mask_size; ++i)
1775 {
1776 mask[i] += subcell_indices[d] * stride;
1777 mask.push_back(stride + mask[i]);
1778 }
1779 stride *= dofs_per_line;
1780 }
1781 return mask;
1782 }
1783 } // namespace DiscreteQuadratureGeneratorImplementation
1784 } // namespace internal
1785
1786
1787
1789 const unsigned int max_box_splits,
1790 const double lower_bound_implicit_function,
1791 const double min_distance_between_roots,
1792 const double limit_to_be_definite,
1793 const double root_finder_tolerance,
1794 const unsigned int max_root_finder_splits,
1795 bool split_in_half)
1796 : max_box_splits(max_box_splits)
1797 , lower_bound_implicit_function(lower_bound_implicit_function)
1798 , min_distance_between_roots(min_distance_between_roots)
1799 , limit_to_be_definite(limit_to_be_definite)
1800 , root_finder_tolerance(root_finder_tolerance)
1801 , max_root_finder_splits(max_root_finder_splits)
1802 , split_in_half(split_in_half)
1803 {}
1804
1805
1806
1807 template <int dim>
1809 const hp::QCollection<1> &q_collection,
1810 const AdditionalData &additional_data)
1811 : q_generator(q_collection, additional_data)
1812 {
1813 Assert(q_collection.size() > 0,
1814 ExcMessage("Incoming hp::QCollection<1> is empty."));
1815 }
1816
1817
1818
1819 template <int dim>
1820 void
1822 {
1823 q_generator.clear_quadratures();
1824 }
1825
1826
1827
1828 template <int dim>
1829 void
1831 const BoundingBox<dim> &box)
1832 {
1833 clear_quadratures();
1834 generate_append(level_set, box);
1835 }
1836
1837
1838
1839 template <int dim>
1840 void
1842 const BoundingBox<dim> &box)
1843 {
1844 Assert(level_set.n_components == 1,
1845 ExcMessage(
1846 "The incoming function should be a scalar level set function,"
1847 " it should have one component."));
1848 Assert(box.volume() > 0, ExcMessage("Incoming box has zero volume."));
1849
1850 std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
1851 level_sets.push_back(level_set);
1852
1853 const unsigned int n_box_splits = 0;
1854 q_generator.generate(level_sets, box, n_box_splits);
1855
1856 // With a single level set function, the "indefinite" quadrature should be
1857 // zero. If you call generate() with a ZeroFunction nothing good can be
1858 // done. You will end up here.
1859 Assert(
1860 q_generator.get_quadratures().indefinite.empty(),
1861 ExcMessage(
1862 "Generation of quadrature rules failed. This can mean that the level "
1863 "set function is degenerate in some way, e.g. oscillating extremely "
1864 "rapidly."));
1865 }
1866
1867
1868
1869 template <int dim>
1870 const Quadrature<dim> &
1872 {
1873 return q_generator.get_quadratures().negative;
1874 }
1875
1876
1877
1878 template <int dim>
1879 const Quadrature<dim> &
1881 {
1882 return q_generator.get_quadratures().positive;
1883 }
1884
1885
1886
1887 template <int dim>
1890 {
1891 return q_generator.get_quadratures().surface;
1892 }
1893
1894
1895 template <int dim>
1896 void
1898 {
1899 q_generator.set_1D_quadrature(q_index);
1900 }
1901
1902
1903
1904 template <int dim>
1906 const hp::QCollection<1> &quadratures1D,
1907 const AdditionalData &additional_data)
1908 : quadrature_generator(quadratures1D, additional_data)
1909 {}
1910
1911
1912
1913 template <int dim>
1914 void
1916 {
1917 quadrature_generator.clear_quadratures();
1918 surface_quadrature.clear();
1919 }
1920
1921
1922
1923 template <int dim>
1924 void
1926 const BoundingBox<dim> &box,
1927 const unsigned int face_index)
1928 {
1929 clear_quadratures();
1930 generate_append(level_set, box, face_index);
1931 }
1932
1933
1934
1935 template <int dim>
1936 void
1938 const BoundingBox<dim> &box,
1939 const unsigned int face_index)
1940 {
1942
1943 // We restrict the level set function to the face, by locking the coordinate
1944 // that is constant over the face. This will be the same as the direction of
1945 // the face normal.
1946 const unsigned int face_normal_direction =
1948
1949 const Point<dim> vertex0 =
1951 const double coordinate_value = vertex0[face_normal_direction];
1952
1953 const Functions::CoordinateRestriction<dim - 1> face_restriction(
1954 level_set, face_normal_direction, coordinate_value);
1955
1956 // Reuse the lower dimensional QuadratureGenerator on the face.
1957 const BoundingBox<dim - 1> cross_section =
1958 box.cross_section(face_normal_direction);
1959
1960 quadrature_generator.generate_append(face_restriction, cross_section);
1961
1962 // We need the dim-dimensional normals of the zero-contour.
1963 // Recompute these.
1964 const ImmersedSurfaceQuadrature<dim - 1, dim - 1>
1965 &surface_quadrature_wrong_normal =
1966 quadrature_generator.get_surface_quadrature();
1967
1968 std::vector<Tensor<1, dim>> normals;
1969 normals.reserve(surface_quadrature_wrong_normal.size());
1970 for (unsigned int i = 0; i < surface_quadrature_wrong_normal.size(); ++i)
1971 {
1973 surface_quadrature_wrong_normal.point(i),
1974 face_normal_direction,
1975 coordinate_value);
1976
1977 Tensor<1, dim> normal = level_set.gradient(point);
1978 normal /= normal.norm();
1979 normals.push_back(normal);
1980 }
1981 surface_quadrature = ImmersedSurfaceQuadrature<dim - 1, dim>(
1982 surface_quadrature_wrong_normal.get_points(),
1983 surface_quadrature_wrong_normal.get_weights(),
1984 normals);
1985 }
1986
1987
1988
1989 template <int dim>
1990 void
1992 {
1993 quadrature_generator.set_1D_quadrature(q_index);
1994 }
1995
1996
1997
1998 template <int dim>
1999 const Quadrature<dim - 1> &
2001 {
2002 return quadrature_generator.get_inside_quadrature();
2003 }
2004
2005
2006 template <int dim>
2007 const Quadrature<dim - 1> &
2009 {
2010 return quadrature_generator.get_outside_quadrature();
2011 }
2012
2013
2014
2015 template <int dim>
2016 const ImmersedSurfaceQuadrature<dim - 1, dim> &
2018 {
2019 return surface_quadrature;
2020 }
2021
2022
2023
2025 const hp::QCollection<1> &quadratures1D,
2026 const AdditionalData &additional_data)
2027 {
2028 (void)quadratures1D;
2029 (void)additional_data;
2030 }
2031
2032
2033
2034 void
2036 {
2037 // do nothing
2038 }
2039
2040
2041
2042 void
2044 const BoundingBox<1> &box,
2045 const unsigned int face_index)
2046 {
2047 clear_quadratures();
2048 generate_append(level_set, box, face_index);
2049 }
2050
2051
2052
2053 void
2055 const BoundingBox<1> &box,
2056 const unsigned int face_index)
2057 {
2059
2060 // The only vertex the 1d-face has.
2061 const Point<1> vertex =
2063
2064 const unsigned int n_points = 1;
2065 const double weight = 1;
2066 const std::vector<Point<0>> points(n_points);
2067 const std::vector<double> weights(n_points, weight);
2068
2069 const double level_set_value = level_set.value(vertex);
2070 if (level_set_value < 0)
2071 {
2072 inside_quadrature = Quadrature<0>(points, weights);
2073 outside_quadrature = Quadrature<0>();
2074 }
2075 else if (level_set_value > 0)
2076 {
2077 inside_quadrature = Quadrature<0>();
2078 outside_quadrature = Quadrature<0>(points, weights);
2079 }
2080 else
2081 {
2082 inside_quadrature = Quadrature<0>();
2083 outside_quadrature = Quadrature<0>();
2084 }
2085 }
2086
2087
2088
2089 void
2091 {
2092 (void)q_index;
2093 }
2094
2095
2096
2097 const Quadrature<0> &
2099 {
2100 return inside_quadrature;
2101 }
2102
2103
2104 const Quadrature<0> &
2106 {
2107 return outside_quadrature;
2108 }
2109
2110
2111
2114 {
2115 return surface_quadrature;
2116 }
2117
2118
2119
2120 template <int dim>
2121 template <typename VectorType>
2123 const hp::QCollection<1> &quadratures1D,
2124 const DoFHandler<dim> &dof_handler,
2125 const VectorType &level_set,
2126 const AdditionalData &additional_data)
2127 : QuadratureGenerator<dim>(quadratures1D, additional_data)
2128 , reference_space_level_set(
2129 std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
2130 RefSpaceFEFieldFunction<dim, VectorType>>(
2131 dof_handler,
2132 level_set))
2133 {}
2134
2135
2136
2137 template <int dim>
2138 void
2140 const BoundingBox<dim> &unit_box)
2141 {
2142 const unsigned int n_subdivisions =
2143 reference_space_level_set->n_subdivisions();
2144
2145 const unsigned int dofs_per_line = n_subdivisions + 1;
2146
2147 this->clear_quadratures();
2148
2149 // The triple nested loop below expands the tensor product of the
2150 // subdivisions, increasing the size of 'all_subcell_indices' by a factor
2151 // 'n_subdivisions' in each iteration. Note how 'n_subdivisions - 1'
2152 // copies of the data are added at the end of the vector in each
2153 // iteration over d. This code avoids a dimension-dependent
2154 // number of nested loops over each direction, yet having the same
2155 // effect.
2156 std::vector<std::array<unsigned int, dim>> all_subcell_indices;
2157 all_subcell_indices.reserve(Utilities::pow(n_subdivisions, dim));
2158 all_subcell_indices.emplace_back();
2159 for (unsigned d = 0; d < dim; ++d)
2160 {
2161 const unsigned int old_size = all_subcell_indices.size();
2162 for (unsigned int i = 1; i < n_subdivisions; ++i)
2163 {
2164 for (unsigned int j = 0; j < old_size; ++j)
2165 {
2166 std::array<unsigned int, dim> next_entry =
2167 all_subcell_indices[j];
2168 next_entry[d] = i;
2169 all_subcell_indices.push_back(next_entry);
2170 }
2171 }
2172 }
2173
2174
2175 for (const auto &subcell_indices : all_subcell_indices)
2176 {
2177 const std::vector<unsigned int> lexicographic_mask =
2178 internal::DiscreteQuadratureGeneratorImplementation::
2179 setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
2180
2181 const auto subcell_box =
2182 internal::DiscreteQuadratureGeneratorImplementation::
2183 create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
2184
2185 reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
2186
2187 QuadratureGenerator<dim>::generate_append(*reference_space_level_set,
2188 subcell_box);
2189 }
2190 }
2191
2192
2193
2194 template <int dim>
2195 void
2197 const typename Triangulation<dim>::active_cell_iterator &cell)
2198 {
2199 Assert(cell->reference_cell().is_hyper_cube(),
2200 internal::DiscreteQuadratureGeneratorImplementation::
2201 ExcReferenceCellNotHypercube());
2202
2203 reference_space_level_set->set_active_cell(cell);
2204 const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
2205 if (reference_space_level_set->is_fe_q_iso_q1())
2206 generate_fe_q_iso_q1(unit_box);
2207 else
2208 QuadratureGenerator<dim>::generate(*reference_space_level_set, unit_box);
2209 }
2210
2211
2212
2213 template <int dim>
2214 template <typename VectorType>
2216 const hp::QCollection<1> &quadratures1D,
2217 const DoFHandler<dim> &dof_handler,
2218 const VectorType &level_set,
2219 const AdditionalData &additional_data)
2220 : FaceQuadratureGenerator<dim>(quadratures1D, additional_data)
2221 , reference_space_level_set(
2222 std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
2223 RefSpaceFEFieldFunction<dim, VectorType>>(
2224 dof_handler,
2225 level_set))
2226 {}
2227
2228
2229
2230 template <int dim>
2231 void
2233 const BoundingBox<dim> &unit_box,
2234 unsigned int face_index)
2235 {
2236 const unsigned int n_subdivisions =
2237 reference_space_level_set->n_subdivisions();
2238
2239 const unsigned int dofs_per_line = n_subdivisions + 1;
2240
2241 this->clear_quadratures();
2242
2243 const unsigned int coordinate_direction_face = face_index / 2;
2244 const unsigned int coordinate_orientation_face = face_index % 2;
2245
2246 // The triple nested loop below expands the tensor product of the
2247 // subdivisions, increasing the size of 'all_subcell_indices' by a factor
2248 // 'n_subdivisions' in each iteration. Note how 'n_subdivisions - 1'
2249 // copies of the data are added at the end of the vector in each
2250 // iteration over d. This code avoids a dimension-dependent
2251 // number of nested loops over each direction, yet having the same
2252 // effect. To account for the face all subcell indices in the coordinate of
2253 // the face direction are filled with the appropriate index depending on the
2254 // coordinate orientation. The other two indices (on the face) are
2255 // set accordingly in the nested loops (with a warp around).
2256 std::vector<std::array<unsigned int, dim>> all_subcell_indices;
2257 all_subcell_indices.reserve(Utilities::pow(n_subdivisions, dim));
2258 all_subcell_indices.emplace_back();
2259 all_subcell_indices[0][coordinate_direction_face] =
2260 coordinate_orientation_face == 0 ? 0 : n_subdivisions - 1;
2261 for (unsigned d = 0; d < dim - 1; ++d)
2262 {
2263 const unsigned int old_size = all_subcell_indices.size();
2264 for (unsigned int i = 1; i < n_subdivisions; ++i)
2265 {
2266 for (unsigned int j = 0; j < old_size; ++j)
2267 {
2268 std::array<unsigned int, dim> next_entry =
2269 all_subcell_indices[j];
2270 next_entry[(coordinate_direction_face + d + 1) % dim] = i;
2271 all_subcell_indices.push_back(next_entry);
2272 }
2273 }
2274 }
2275
2276 for (const auto &subcell_indices : all_subcell_indices)
2277 {
2278 const std::vector<unsigned int> lexicographic_mask =
2279 internal::DiscreteQuadratureGeneratorImplementation::
2280 setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
2281
2282 const auto subcell_box =
2283 internal::DiscreteQuadratureGeneratorImplementation::
2284 create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
2285
2286 reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
2287
2289 *reference_space_level_set, subcell_box, face_index);
2290 }
2291 }
2292
2293
2294
2295 template <int dim>
2296 void
2298 const typename Triangulation<dim>::active_cell_iterator &cell,
2299 const unsigned int face_index)
2300 {
2301 Assert(cell->reference_cell().is_hyper_cube(),
2302 internal::DiscreteQuadratureGeneratorImplementation::
2303 ExcReferenceCellNotHypercube());
2304
2305 reference_space_level_set->set_active_cell(cell);
2306 const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
2307 if (reference_space_level_set->is_fe_q_iso_q1())
2308 generate_fe_q_iso_q1(unit_box, face_index);
2309 else
2310 FaceQuadratureGenerator<dim>::generate(*reference_space_level_set,
2311 unit_box,
2312 face_index);
2313 }
2314} // namespace NonMatching
2315#include "quadrature_generator.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
BoundingBox< 1, Number > bounds(const unsigned int direction) const
Point< spacedim, Number > center() const
Number lower_bound(const unsigned int direction) const
std::pair< Point< spacedim, Number >, Point< spacedim, Number > > & get_boundary_points()
double volume() const
Number side_length(const unsigned int direction) const
BoundingBox< spacedim, Number > child(const unsigned int index) const
Point< spacedim, Number > vertex(const unsigned int index) const
BoundingBox< spacedim - 1, Number > cross_section(const unsigned int direction) const
Number upper_bound(const unsigned int direction) const
Point< spacedim, Number > unit_to_real(const Point< spacedim, Number > &point) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
unsigned int n_base_elements() const
virtual Tensor< 1, dim, RangeNumberType > gradient(const Point< dim > &p, const unsigned int component=0) const
const unsigned int n_components
Definition function.h:163
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
DiscreteFaceQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const DoFHandler< dim > &dof_handler, const VectorType &level_set, const AdditionalData &additional_data=AdditionalData())
void generate(const typename Triangulation< dim >::active_cell_iterator &cell, const unsigned int face_index)
void generate_fe_q_iso_q1(const BoundingBox< dim > &unit_box, unsigned int face_index)
void generate(const typename Triangulation< dim >::active_cell_iterator &cell)
void generate_fe_q_iso_q1(const BoundingBox< dim > &unit_box)
DiscreteQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const DoFHandler< dim > &dof_handler, const VectorType &level_set, const AdditionalData &additional_data=AdditionalData())
const Quadrature< dim - 1 > & get_inside_quadrature() const
void generate(const Function< dim > &level_set, const BoundingBox< dim > &box, const unsigned int face_index)
void set_1D_quadrature(const unsigned int q_index)
void generate_append(const Function< dim > &level_set, const BoundingBox< dim > &box, const unsigned int face_index)
FaceQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const AdditionalData &additional_data=AdditionalData())
const ImmersedSurfaceQuadrature< dim - 1, dim > & get_surface_quadrature() const
const Quadrature< dim - 1 > & get_outside_quadrature() const
void push_back(const Point< dim > &point, const double weight, const Tensor< 1, spacedim > &normal)
void generate_append(const Function< dim > &level_set, const BoundingBox< dim > &box)
void set_1D_quadrature(const unsigned int q_index)
QuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const AdditionalData &additional_data=AdditionalData())
const Quadrature< dim > & get_inside_quadrature() const
const Quadrature< dim > & get_outside_quadrature() const
const ImmersedSurfaceQuadrature< dim > & get_surface_quadrature() const
void generate(const Function< dim > &level_set, const BoundingBox< dim > &box)
void set_active_cell(const typename Triangulation< dim >::active_cell_iterator &cell) override
Tensor< 1, dim > gradient(const Point< dim > &point, const unsigned int component=0) const override
double value(const Point< dim > &point, const unsigned int component=0) const override
void set_subcell(const std::vector< unsigned int > &mask, const BoundingBox< dim > &subcell_box) override
RefSpaceFEFieldFunction(const DoFHandler< dim > &dof_handler, const VectorType &dof_values)
SymmetricTensor< 2, dim > hessian(const Point< dim > &point, const unsigned int component=0) const override
QGeneratorBase(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
QGenerator(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
void create_high_dim_quadratures(const unsigned int height_function_direction, const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box)
void generate(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const unsigned int n_box_splits)
void split_box_and_recurse(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const std::optional< HeightDirectionData > &direction_data, const unsigned int n_box_splits)
void use_midpoint_method(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box)
void create_low_dim_quadratures(const unsigned int height_function_direction, const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const unsigned int n_box_splits)
ExtendableQuadrature< dim > & quadrature_by_definiteness(const Definiteness definiteness)
void find_roots(const std::vector< std::reference_wrapper< const Function< 1 > > > &functions, const BoundingBox< 1 > &interval, std::vector< double > &roots)
void generate(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const Quadrature< dim - 1 > &low_dim_quadrature, const unsigned int height_function_direction, QPartitioning< dim > &q_partitioning)
void create_surface_point(const Point< dim - 1 > &point, const double weight, const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const unsigned int height_function_direction, ImmersedSurfaceQuadrature< dim > &surface_quadrature)
UpThroughDimensionCreator(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
Definition point.h:111
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
unsigned int n_vertices() const
numbers::NumberTraits< Number >::real_type norm() const
unsigned int size() const
Definition collection.h:264
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > center
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< unsigned int > setup_lexicographic_mask(const std::array< unsigned int, dim > &subcell_indices, unsigned int dofs_per_line)
BoundingBox< dim > create_subcell_box(const BoundingBox< dim > &unit_box, const std::array< unsigned int, dim > &subcell_indices, const unsigned int n_subdivisions)
void distribute_points_between_roots(const Quadrature< 1 > &quadrature1D, const BoundingBox< 1 > &interval, const std::vector< double > &roots, const Point< dim - 1 > &point, const double weight, const unsigned int height_function_direction, const std::vector< std::reference_wrapper< const Function< 1 > > > &level_sets, const AdditionalQGeneratorData &additional_data, QPartitioning< dim > &q_partitioning)
double lower_bound_on_abs(const std::pair< double, double > &function_bounds)
void estimate_function_bounds(const std::vector< std::reference_wrapper< const Function< dim > > > &functions, const BoundingBox< dim > &box, std::vector< FunctionBounds< dim > > &all_function_bounds)
BoundingBox< dim > left_half(const BoundingBox< dim > &box, const unsigned int direction)
std::optional< unsigned int > direction_of_largest_extent(const BoundingBox< dim > &box)
void add_tensor_product(const Quadrature< dim - 1 > &lower, const Quadrature< 1 > &quadrature1D, const double start, const double end, const unsigned int component_in_dim, ExtendableQuadrature< dim > &quadrature)
bool one_positive_one_negative_definite(const std::vector< FunctionBounds< dim > > &all_function_bounds)
void map_quadrature_to_box(const Quadrature< dim > &unit_quadrature, const BoundingBox< dim > &box, ExtendableQuadrature< dim > &quadrature)
std::pair< double, double > find_extreme_values(const std::vector< FunctionBounds< dim > > &all_function_bounds)
std::optional< HeightDirectionData > find_best_height_direction(const std::vector< FunctionBounds< dim > > &all_function_bounds)
void tensor_point_with_1D_quadrature(const Point< dim - 1 > &point, const double weight, const Quadrature< 1 > &quadrature1D, const double start, const double end, const unsigned int component_in_dim, ExtendableQuadrature< dim > &quadrature)
void restrict_to_point(const std::vector< std::reference_wrapper< const Function< dim > > > &functions, const Point< dim - 1 > &point, const unsigned int open_direction, std::vector< Functions::PointRestriction< dim - 1 > > &restrictions)
BoundingBox< dim > right_half(const BoundingBox< dim > &box, const unsigned int direction)
bool is_indefinite(const std::pair< double, double > &function_bounds)
void take_min_max_at_vertices(const Function< dim > &function, const BoundingBox< dim > &box, std::pair< double, double > &value_bounds)
unsigned int compute_split_direction(const BoundingBox< dim > &box, const std::optional< HeightDirectionData > &height_direction_data)
void restrict_to_top_and_bottom(const std::vector< std::reference_wrapper< const Function< dim > > > &functions, const BoundingBox< dim > &box, const unsigned int direction, std::vector< Functions::CoordinateRestriction< dim - 1 > > &restrictions)
Point< dim > face_projection_closest_zero_contour(const Point< dim - 1 > &point, const unsigned int direction, const BoundingBox< dim > &box, const Function< dim > &level_set)
Definiteness pointwise_definiteness(const std::vector< std::reference_wrapper< const Function< dim > > > &functions, const Point< dim > &point)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
std::vector< Polynomials::Polynomial< double > > get_polynomial_space(const FiniteElement< dim, spacedim > &fe)
bool is_fast_path_supported(const FiniteElement< dim, spacedim > &fe, const unsigned int base_element_number)
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
Point< dim+1 > create_higher_dim_point(const Point< dim > &point, const unsigned int component_in_dim_plus_1, const double coordinate_value)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
static const unsigned int invalid_unsigned_int
Definition types.h:220
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
AdditionalQGeneratorData(const unsigned int max_box_splits=4, const double lower_bound_implicit_function=1e-11, const double min_distance_between_roots=1e-12, const double limit_to_be_definite=1e-11, const double root_finder_tolerance=1e-12, const unsigned int max_root_finder_splits=2, bool split_in_half=true)
AdditionalData(const double tolerance=1e-12, const unsigned int max_recursion_depth=2, const unsigned int max_iterations=500)
void reinit(const Quadrature< dim_q > &quad, const FiniteElement< dim, spacedim > &fe_dim, const unsigned int base_element=0)
std::vector< unsigned int > lexicographic_numbering
Definition shape_info.h:479
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)