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