Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2021 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
19
30#include <deal.II/lac/vector.h>
31
33
35
36#include <boost/math/special_functions/relative_difference.hpp>
37#include <boost/math/special_functions/sign.hpp>
38#include <boost/math/tools/roots.hpp>
39
40#include <algorithm>
41#include <vector>
42
44
45namespace NonMatching
46{
47 namespace internal
48 {
49 namespace QuadratureGeneratorImplementation
50 {
51 template <int dim>
52 void
54 const double weight,
55 const Quadrature<1> & quadrature1D,
56 const double start,
57 const double end,
58 const unsigned int component_in_dim,
59 ExtendableQuadrature<dim> &quadrature)
60 {
61 Assert(start < end,
62 ExcMessage("Interval start must be less than interval end."));
63
64 const double length = end - start;
65 for (unsigned int j = 0; j < quadrature1D.size(); ++j)
66 {
67 const double x = start + length * quadrature1D.point(j)[0];
69 point, component_in_dim, x),
70 length * weight * quadrature1D.weight(j));
71 }
72 }
73
74
75
80 template <int dim>
81 void
83 const Quadrature<1> & quadrature1D,
84 const double start,
85 const double end,
86 const unsigned int component_in_dim,
87 ExtendableQuadrature<dim> &quadrature)
88 {
89 for (unsigned int j = 0; j < lower.size(); ++j)
90 {
92 lower.weight(j),
93 quadrature1D,
94 start,
95 end,
96 component_in_dim,
97 quadrature);
98 }
99 }
100
101
102
103 template <int dim>
106 const std::vector<std::reference_wrapper<const Function<dim>>>
107 & functions,
108 const Point<dim> &point)
109 {
110 Assert(functions.size() > 0,
112 "The incoming vector must contain at least one function."));
113
114 const int sign_of_first =
115 boost::math::sign(functions[0].get().value(point));
116
117 if (sign_of_first == 0)
119
120 for (unsigned int j = 1; j < functions.size(); ++j)
121 {
122 const int sign = boost::math::sign(functions[j].get().value(point));
123
124 if (sign != sign_of_first)
126 }
127 // If we got here all functions have the same sign.
128 if (sign_of_first < 0)
130 else
132 }
133
134
135
148 template <int dim>
149 void
151 const BoundingBox<dim> & box,
152 std::pair<double, double> &value_bounds)
153 {
154 const ReferenceCell &cube = ReferenceCells::get_hypercube<dim>();
155 for (unsigned int i = 0; i < cube.n_vertices(); ++i)
156 {
157 const double vertex_value = function.value(box.vertex(i));
158
159 value_bounds.first = std::min(value_bounds.first, vertex_value);
160 value_bounds.second = std::max(value_bounds.second, vertex_value);
161 }
162 }
163
164
165
176 template <int dim>
177 void
179 const std::vector<std::reference_wrapper<const Function<dim>>>
180 & functions,
181 const BoundingBox<dim> & box,
182 std::vector<FunctionBounds<dim>> &all_function_bounds)
183 {
184 all_function_bounds.clear();
185 all_function_bounds.reserve(functions.size());
186 for (const Function<dim> &function : functions)
187 {
188 FunctionBounds<dim> bounds;
189 FunctionTools::taylor_estimate_function_bounds<dim>(
190 function, box, bounds.value, bounds.gradient);
191 take_min_max_at_vertices(function, box, bounds.value);
192
193 all_function_bounds.push_back(bounds);
194 }
195 }
196
197
198
199 template <int dim>
200 std::pair<double, double>
201 find_extreme_values(const std::vector<FunctionBounds<dim>> &bounds)
202 {
203 Assert(bounds.size() > 0, ExcMessage("The incoming vector is empty."));
204
205 std::pair<double, double> extremes = bounds[0].value;
206 for (unsigned int i = 1; i < bounds.size(); ++i)
207 {
208 extremes.first = std::min(extremes.first, bounds[i].value.first);
209 extremes.second = std::max(extremes.second, bounds[i].value.second);
210 }
211
212 return extremes;
213 }
214
215
216
221 inline bool
222 is_indefinite(const std::pair<double, double> &function_bounds)
223 {
224 if (function_bounds.first > 0)
225 return false;
226 if (function_bounds.second < 0)
227 return false;
228 return true;
229 }
230
231
232
253 inline double
254 lower_bound_on_abs(const std::pair<double, double> &function_bounds)
255 {
256 Assert(function_bounds.first <= function_bounds.second,
257 ExcMessage("Function bounds reversed, max < min."));
258
259 return 0.5 * (std::abs(function_bounds.second + function_bounds.first) -
260 (function_bounds.second - function_bounds.first));
261 }
262
263
264
266 {
268 min_abs_dfdx = 0;
269 }
270
271
272
273 template <int dim>
274 std_cxx17::optional<HeightDirectionData>
276 const std::vector<FunctionBounds<dim>> &all_function_bounds)
277 {
278 // Minimum (taken over the indefinite functions) on the lower bound on
279 // each component of the gradient.
280 std_cxx17::optional<std::array<double, dim>> min_lower_abs_grad;
281
282 for (const FunctionBounds<dim> &bounds : all_function_bounds)
283 {
284 if (is_indefinite(bounds.value))
285 {
286 // For the first indefinite function we find, we write the lower
287 // bounds on each gradient component to min_lower_abs_grad.
288 if (!min_lower_abs_grad)
289 {
290 min_lower_abs_grad.emplace();
291 for (unsigned int d = 0; d < dim; ++d)
292 {
293 (*min_lower_abs_grad)[d] =
294 lower_bound_on_abs(bounds.gradient[d]);
295 }
296 }
297 else
298 {
299 for (unsigned int d = 0; d < dim; ++d)
300 {
301 (*min_lower_abs_grad)[d] =
302 std::min((*min_lower_abs_grad)[d],
303 lower_bound_on_abs(bounds.gradient[d]));
304 }
305 }
306 }
307 }
308
309 if (min_lower_abs_grad)
310 {
311 const auto max_element =
312 std::max_element(min_lower_abs_grad->begin(),
313 min_lower_abs_grad->end());
314
316 data.direction =
317 std::distance(min_lower_abs_grad->begin(), max_element);
318 data.min_abs_dfdx = *max_element;
319
320 return data;
321 }
322
323 return std_cxx17::optional<HeightDirectionData>();
324 }
325
326
327
333 template <int dim>
334 inline bool
336 const std::vector<FunctionBounds<dim>> &all_function_bounds)
337 {
338 if (all_function_bounds.size() != 2)
339 return false;
340 else
341 {
342 const FunctionBounds<dim> &bounds0 = all_function_bounds.at(0);
343 const FunctionBounds<dim> &bounds1 = all_function_bounds.at(1);
344
345 if (bounds0.value.first > 0 && bounds1.value.second < 0)
346 return true;
347 if (bounds1.value.first > 0 && bounds0.value.second < 0)
348 return true;
349 return false;
350 }
351 }
352
353
354
362 template <int dim>
363 void
364 map_quadrature_to_box(const Quadrature<dim> & unit_quadrature,
365 const BoundingBox<dim> & box,
366 ExtendableQuadrature<dim> &quadrature)
367 {
368 for (unsigned int i = 0; i < unit_quadrature.size(); ++i)
369 {
370 const Point<dim> point = box.unit_to_real(unit_quadrature.point(i));
371 const double weight = unit_quadrature.weight(i) * box.volume();
372
373 quadrature.push_back(point, weight);
374 }
375 }
376
377
378
389 template <int dim>
390 void
392 const std::vector<std::reference_wrapper<const Function<dim>>>
393 & functions,
394 const BoundingBox<dim> & box,
395 const unsigned int direction,
396 std::vector<Functions::CoordinateRestriction<dim - 1>> &restrictions)
397 {
398 AssertIndexRange(direction, dim);
399
400 restrictions.clear();
401 restrictions.reserve(2 * functions.size());
402
403 const double bottom = box.lower_bound(direction);
404 const double top = box.upper_bound(direction);
405
406 for (const auto &function : functions)
407 {
408 restrictions.push_back(Functions::CoordinateRestriction<dim - 1>(
409 function, direction, bottom));
410 restrictions.push_back(Functions::CoordinateRestriction<dim - 1>(
411 function, direction, top));
412 }
413 }
414
415
416
425 template <int dim>
426 void
428 const std::vector<std::reference_wrapper<const Function<dim>>>
429 & functions,
430 const Point<dim - 1> & point,
431 const unsigned int open_direction,
432 std::vector<Functions::PointRestriction<dim - 1>> &restrictions)
433 {
434 AssertIndexRange(open_direction, dim);
435
436 restrictions.clear();
437 restrictions.reserve(functions.size());
438 for (const auto &function : functions)
439 {
440 restrictions.push_back(Functions::PointRestriction<dim - 1>(
441 function, open_direction, point));
442 }
443 }
444
445
446
460 template <int dim>
461 void
463 const Quadrature<1> & quadrature1D,
464 const BoundingBox<1> & interval,
465 const std::vector<double> &roots,
466 const Point<dim - 1> & point,
467 const double weight,
468 const unsigned int height_function_direction,
469 const std::vector<std::reference_wrapper<const Function<1>>>
470 & level_sets,
471 const AdditionalQGeneratorData &additional_data,
472 QPartitioning<dim> & q_partitioning)
473 {
474 // Make this int to avoid a warning signed/unsigned comparison.
475 const int n_roots = roots.size();
476
477 // The number of intervals are roots.size() + 1
478 for (int i = -1; i < n_roots; ++i)
479 {
480 // Start and end point of the subinterval.
481 const double start = i < 0 ? interval.lower_bound(0) : roots[i];
482 const double end =
483 i + 1 < n_roots ? roots[i + 1] : interval.upper_bound(0);
484
485 const double length = end - start;
486 // It might be that the end points of the subinterval are roots.
487 // If this is the case then the subinterval has length zero.
488 // Don't distribute points on the subinterval if it is shorter than
489 // some tolerance.
490 if (length > additional_data.min_distance_between_roots)
491 {
492 // All points on the interval belong to the same region in
493 // the QPartitioning. Determine the quadrature we should add
494 // the points to.
495 const Point<1> center(start + 0.5 * length);
496 const Definiteness definiteness =
497 pointwise_definiteness(level_sets, center);
498 ExtendableQuadrature<dim> &target_quadrature =
499 q_partitioning.quadrature_by_definiteness(definiteness);
500
502 weight,
503 quadrature1D,
504 start,
505 end,
506 height_function_direction,
507 target_quadrature);
508 }
509 }
510 }
511
512
513
515 const double tolerance,
516 const unsigned int max_recursion_depth,
517 const unsigned int max_iterations)
518 : tolerance(tolerance)
519 , max_recursion_depth(max_recursion_depth)
520 , max_iterations(max_iterations)
521 {}
522
523
524
526 : additional_data(data)
527 {}
528
529
530
531 void
533 const std::vector<std::reference_wrapper<const Function<1>>> &functions,
534 const BoundingBox<1> & interval,
535 std::vector<double> & roots)
536 {
537 for (const Function<1> &function : functions)
538 {
539 const unsigned int recursion_depth = 0;
540 find_roots(function, interval, recursion_depth, roots);
541 }
542 // Sort and make sure no roots are duplicated
543 std::sort(roots.begin(), roots.end());
544
545 const auto roots_are_equal = [this](const double &a, const double &b) {
546 return std::abs(a - b) < additional_data.tolerance;
547 };
548 roots.erase(unique(roots.begin(), roots.end(), roots_are_equal),
549 roots.end());
550 }
551
552
553
554 void
556 const BoundingBox<1> &interval,
557 const unsigned int recursion_depth,
558 std::vector<double> & roots)
559 {
560 // Compute function values at end points.
561 const double left_value = function.value(interval.vertex(0));
562 const double right_value = function.value(interval.vertex(1));
563
564 // If we have a sign change we solve for the root.
565 if (boost::math::sign(left_value) != boost::math::sign(right_value))
566 {
567 const auto lambda = [&function](const double x) {
568 return function.value(Point<1>(x));
569 };
570
571 const auto stopping_criteria = [this](const double &a,
572 const double &b) {
573 return std::abs(a - b) < additional_data.tolerance;
574 };
575
576 boost::uintmax_t iterations = additional_data.max_iterations;
577
578 const std::pair<double, double> root_bracket =
579 boost::math::tools::toms748_solve(lambda,
580 interval.lower_bound(0),
581 interval.upper_bound(0),
582 left_value,
583 right_value,
584 stopping_criteria,
585 iterations);
586
587 const double root = .5 * (root_bracket.first + root_bracket.second);
588 roots.push_back(root);
589 }
590 else
591 {
592 // Compute bounds on the incoming function to check if there are
593 // roots. If the function is positive or negative on the whole
594 // interval we do nothing.
595 std::pair<double, double> value_bounds;
596 std::array<std::pair<double, double>, 1> gradient_bounds;
597 FunctionTools::taylor_estimate_function_bounds<1>(function,
598 interval,
599 value_bounds,
600 gradient_bounds);
601
602 // Since we already know the function values at the interval ends we
603 // might as well check these for min/max too.
604 const double function_min =
605 std::min(std::min(left_value, right_value), value_bounds.first);
606
607 // If the functions is positive there are no roots.
608 if (function_min > 0)
609 return;
610
611 const double function_max =
612 std::max(std::max(left_value, right_value), value_bounds.second);
613
614 // If the function is negative there are no roots.
615 if (function_max < 0)
616 return;
617
618 // If we can't say that the function is strictly positive/negative
619 // we split the interval in half. We can't split forever, so if we
620 // have reached the max recursion, we stop looking for roots.
621 if (recursion_depth < additional_data.max_recursion_depth)
622 {
623 find_roots(function,
624 interval.child(0),
625 recursion_depth + 1,
626 roots);
627 find_roots(function,
628 interval.child(1),
629 recursion_depth + 1,
630 roots);
631 }
632 }
633 }
634
635
636
637 template <int dim>
639 const Quadrature<dim> &quadrature)
640 : Quadrature<dim>(quadrature)
641 {}
642
643
644
645 template <int dim>
646 void
648 const double weight)
649 {
650 this->quadrature_points.push_back(point);
651 this->weights.push_back(weight);
652 }
653
654
655
656 template <int dim>
659 const Definiteness definiteness)
660 {
661 switch (definiteness)
662 {
664 return negative;
666 return positive;
667 default:
668 return indefinite;
669 }
670 }
671
672
673
681 template <int dim>
684 const unsigned int direction,
685 const BoundingBox<dim> &box,
686 const Function<dim> & level_set)
687 {
688 const Point<dim> bottom_point =
690 direction,
691 box.lower_bound(direction));
692 const double bottom_value = level_set.value(bottom_point);
693
694 const Point<dim> top_point =
696 direction,
697 box.upper_bound(direction));
698 const double top_value = level_set.value(top_point);
699
700 // The end point closest to the zero-contour is the one with smallest
701 // absolute value.
702 return std::abs(bottom_value) < std::abs(top_value) ? bottom_point :
703 top_point;
704 }
705
706
707
708 template <int dim, int spacedim>
710 const hp::QCollection<1> & q_collection1D,
711 const AdditionalQGeneratorData &additional_data)
712 : q_collection1D(&q_collection1D)
713 , additional_data(additional_data)
714 , root_finder(
715 RootFinder::AdditionalData(additional_data.root_finder_tolerance,
716 additional_data.max_root_finder_splits))
717 {
718 q_index = 0;
719 }
720
721
722
723 template <int dim, int spacedim>
724 void
726 const std::vector<std::reference_wrapper<const Function<dim>>>
727 & level_sets,
728 const BoundingBox<dim> & box,
729 const Quadrature<dim - 1> &low_dim_quadrature,
730 const unsigned int height_function_direction,
731 QPartitioning<dim> & q_partitioning)
732 {
733 const Quadrature<1> &quadrature1D = (*q_collection1D)[q_index];
734
735 for (unsigned int q = 0; q < low_dim_quadrature.size(); ++q)
736 {
737 const Point<dim - 1> &point = low_dim_quadrature.point(q);
738 const double weight = low_dim_quadrature.weight(q);
739 restrict_to_point(level_sets,
740 point,
741 height_function_direction,
742 point_restrictions);
743
744 // We need a vector of references to do the recursive call.
745 const std::vector<std::reference_wrapper<const Function<1>>>
746 restrictions(point_restrictions.begin(),
747 point_restrictions.end());
748
749 const BoundingBox<1> bounds_in_direction =
750 box.bounds(height_function_direction);
751
752 roots.clear();
753 root_finder.find_roots(restrictions, bounds_in_direction, roots);
754
756 bounds_in_direction,
757 roots,
758 point,
759 weight,
760 height_function_direction,
761 restrictions,
762 additional_data,
763 q_partitioning);
764
765 if (dim == spacedim)
766 create_surface_point(point,
767 weight,
768 level_sets,
769 box,
770 height_function_direction,
771 q_partitioning.surface);
772 }
773
774 point_restrictions.clear();
775 }
776
777
778
779 template <int dim, int spacedim>
780 void
782 const Point<dim - 1> &point,
783 const double weight,
784 const std::vector<std::reference_wrapper<const Function<dim>>>
785 & level_sets,
786 const BoundingBox<dim> & box,
787 const unsigned int height_function_direction,
788 ImmersedSurfaceQuadrature<dim> &surface_quadrature)
789 {
790 AssertIndexRange(roots.size(), 2);
791 Assert(level_sets.size() == 1, ExcInternalError());
792
793
794 const Function<dim> &level_set = level_sets.at(0);
795
796 Point<dim> surface_point;
797 if (roots.size() == 1)
798 {
800 point, height_function_direction, roots[0]);
801 }
802 else
803 {
804 // If we got here, we have missed roots in the lower dimensional
805 // algorithm. This is a rare event but can happen if the
806 // zero-contour has a high curvature. The problem is that the
807 // incoming point has been incorrectly added to the indefinite
808 // quadrature in QPartitioning<dim-1>. Since we missed a root on
809 // this box, we will likely miss it on the neighboring box too. If
810 // this happens, the point will NOT be in the indefinite quadrature
811 // on the neighbor. The best thing we can do is to compute the
812 // surface point by projecting the lower dimensional point to the
813 // face closest to the zero-contour. We should add a surface point
814 // because the neighbor will not.
816 point, height_function_direction, box, level_set);
817 }
818
819 const Tensor<1, dim> gradient = level_set.gradient(surface_point);
820 Tensor<1, dim> normal = gradient;
821 normal *= 1. / normal.norm();
822
823 // Note that gradient[height_function_direction] is non-zero
824 // because of the implicit function theorem.
825 const double surface_weight =
826 weight * gradient.norm() /
827 std::abs(gradient[height_function_direction]);
828 surface_quadrature.push_back(surface_point, surface_weight, normal);
829 }
830
831
832
833 template <int dim, int spacedim>
834 void
836 unsigned int q_index)
837 {
838 AssertIndexRange(q_index, q_collection1D->size());
839 this->q_index = q_index;
840 }
841
842
843
844 template <int dim, int spacedim>
846 const hp::QCollection<1> & q_collection1D,
847 const AdditionalQGeneratorData &additional_data)
848 : additional_data(additional_data)
849 , q_collection1D(&q_collection1D)
850 {
851 q_index = 0;
852 }
853
854
855
856 template <int dim, int spacedim>
858 const hp::QCollection<1> & q_collection1D,
859 const AdditionalQGeneratorData &additional_data)
860 : QGeneratorBase<dim, spacedim>(q_collection1D, additional_data)
861 , low_dim_algorithm(q_collection1D, additional_data)
862 , up_through_dimension_creator(q_collection1D, additional_data)
863 {
864 for (const auto &quadrature : q_collection1D)
865 tensor_products.push_back(quadrature);
866 }
867
868
869
870 template <int dim, int spacedim>
871 void
873 {
874 q_partitioning = QPartitioning<dim>();
875 }
876
877
878
879 template <int dim, int spacedim>
880 const QPartitioning<dim> &
882 {
883 return q_partitioning;
884 }
885
886
887
888 template <int dim, int spacedim>
889 void
891 const std::vector<std::reference_wrapper<const Function<dim>>>
892 & level_sets,
893 const BoundingBox<dim> &box,
894 const unsigned int n_box_splits)
895 {
896 std::vector<FunctionBounds<dim>> all_function_bounds;
897 estimate_function_bounds(level_sets, box, all_function_bounds);
898
899 const std::pair<double, double> extreme_values =
900 find_extreme_values(all_function_bounds);
901
902 if (extreme_values.first > this->additional_data.limit_to_be_definite)
903 {
904 map_quadrature_to_box(tensor_products[this->q_index],
905 box,
906 this->q_partitioning.positive);
907 }
908 else if (extreme_values.second <
909 -(this->additional_data.limit_to_be_definite))
910 {
911 map_quadrature_to_box(tensor_products[this->q_index],
912 box,
913 this->q_partitioning.negative);
914 }
915 else if (one_positive_one_negative_definite(all_function_bounds))
916 {
917 map_quadrature_to_box(tensor_products[this->q_index],
918 box,
919 this->q_partitioning.indefinite);
920 }
921 else
922 {
923 const std_cxx17::optional<HeightDirectionData> data =
924 find_best_height_direction(all_function_bounds);
925
926 // Check larger than a constant to avoid that min_abs_dfdx is only
927 // larger by 0 by floating point precision.
928 if (data && data->min_abs_dfdx >
929 this->additional_data.lower_bound_implicit_function)
930 {
931 create_low_dim_quadratures(data->direction,
932 level_sets,
933 box,
934 n_box_splits);
935 create_high_dim_quadratures(data->direction, level_sets, box);
936 }
937 else if (n_box_splits < this->additional_data.max_box_splits)
938 {
939 split_box_and_recurse(level_sets, box, data, n_box_splits);
940 }
941 else
942 {
943 // We can't split the box recursively forever. Use the midpoint
944 // method as a last resort.
945 use_midpoint_method(level_sets, box);
946 }
947 }
948 }
949
950
951
957 template <int dim>
958 std_cxx17::optional<unsigned int>
960 {
961 // Get the side lengths for each direction and sort them.
962 std::array<std::pair<double, unsigned int>, dim> side_lengths;
963 for (unsigned int i = 0; i < dim; ++i)
964 {
965 side_lengths[i].first = box.side_length(i);
966 side_lengths[i].second = i;
967 }
968 // Sort is lexicographic, so this sorts based on side length first.
969 std::sort(side_lengths.begin(), side_lengths.end());
970
971 // Check if the two largest side lengths have the same length. This
972 // function isn't called in 1d, so the (dim - 2)-element exists.
973 if (boost::math::epsilon_difference(side_lengths[dim - 1].first,
974 side_lengths[dim - 2].first) < 100)
975 return std_cxx17::optional<unsigned int>();
976
977 return side_lengths.back().second;
978 }
979
980
981
993 template <int dim>
994 unsigned int
996 const BoundingBox<dim> & box,
997 const std_cxx17::optional<HeightDirectionData> &height_direction_data)
998 {
999 const std_cxx17::optional<unsigned int> direction =
1001
1002 if (direction)
1003 return *direction;
1004
1005 // This direction is closest to being a height direction, so
1006 // we split in this direction.
1007 if (height_direction_data)
1008 return height_direction_data->direction;
1009
1010 // We have to choose some direction, we might as well take 0.
1011 return 0;
1012 }
1013
1014
1015
1020 template <int dim>
1021 inline BoundingBox<dim>
1022 left_half(const BoundingBox<dim> &box, const unsigned int direction)
1023 {
1024 AssertIndexRange(direction, dim);
1025
1026 // Move the upper corner half a side-length to the left.
1027 std::pair<Point<dim>, Point<dim>> corners = box.get_boundary_points();
1028 corners.second[direction] -= .5 * box.side_length(direction);
1029
1030 return BoundingBox<dim>(corners);
1031 }
1032
1033
1034
1039 template <int dim>
1040 inline BoundingBox<dim>
1041 right_half(const BoundingBox<dim> &box, const unsigned int direction)
1042 {
1043 AssertIndexRange(direction, dim);
1044
1045 // Move the lower corner half a side-length to the right.
1046 std::pair<Point<dim>, Point<dim>> corners = box.get_boundary_points();
1047 corners.first[direction] += .5 * box.side_length(direction);
1048
1049 return BoundingBox<dim>(corners);
1050 }
1051
1052
1053
1054 template <int dim, int spacedim>
1055 void
1057 const std::vector<std::reference_wrapper<const Function<dim>>>
1058 & level_sets,
1059 const BoundingBox<dim> & box,
1060 const std_cxx17::optional<HeightDirectionData> &direction_data,
1061 const unsigned int n_box_splits)
1062 {
1063 if (this->additional_data.split_in_half)
1064 {
1065 const unsigned int direction =
1066 compute_split_direction(box, direction_data);
1067
1068 const BoundingBox<dim> left_box = left_half(box, direction);
1069 const BoundingBox<dim> right_box = right_half(box, direction);
1070
1071 generate(level_sets, left_box, n_box_splits + 1);
1072 generate(level_sets, right_box, n_box_splits + 1);
1073 }
1074 else
1075 {
1076 for (unsigned int i = 0;
1077 i < GeometryInfo<dim>::max_children_per_cell;
1078 ++i)
1079 {
1080 generate(level_sets, box.child(i), n_box_splits + 1);
1081 }
1082 }
1083 }
1084
1085
1086
1087 template <int dim, int spacedim>
1088 void
1090 const unsigned int height_function_direction,
1091 const std::vector<std::reference_wrapper<const Function<dim>>>
1092 & level_sets,
1093 const BoundingBox<dim> &box,
1094 const unsigned int n_box_splits)
1095 {
1096 std::vector<Functions::CoordinateRestriction<dim - 1>>
1097 face_restrictions;
1098 restrict_to_top_and_bottom(level_sets,
1099 box,
1100 height_function_direction,
1101 face_restrictions);
1102
1103 // We need a vector of references to do the recursive call.
1104 const std::vector<std::reference_wrapper<const Function<dim - 1>>>
1105 restrictions(face_restrictions.begin(), face_restrictions.end());
1106
1107 const BoundingBox<dim - 1> cross_section =
1108 box.cross_section(height_function_direction);
1109
1110 low_dim_algorithm.clear_quadratures();
1111 low_dim_algorithm.generate(restrictions, cross_section, n_box_splits);
1112 }
1113
1114
1116 template <int dim, int spacedim>
1117 void
1119 const unsigned int height_function_direction,
1120 const std::vector<std::reference_wrapper<const Function<dim>>>
1121 & level_sets,
1122 const BoundingBox<dim> &box)
1123 {
1124 const QPartitioning<dim - 1> &low_dim_quadratures =
1125 low_dim_algorithm.get_quadratures();
1126
1127 const Quadrature<1> &quadrature1D =
1128 (*this->q_collection1D)[this->q_index];
1129
1130 add_tensor_product(low_dim_quadratures.negative,
1131 quadrature1D,
1132 box.lower_bound(height_function_direction),
1133 box.upper_bound(height_function_direction),
1134 height_function_direction,
1135 this->q_partitioning.negative);
1136
1137 add_tensor_product(low_dim_quadratures.positive,
1138 quadrature1D,
1139 box.lower_bound(height_function_direction),
1140 box.upper_bound(height_function_direction),
1141 height_function_direction,
1142 this->q_partitioning.positive);
1143
1144 up_through_dimension_creator.generate(level_sets,
1145 box,
1146 low_dim_quadratures.indefinite,
1147 height_function_direction,
1148 this->q_partitioning);
1149 }
1151
1152
1153 template <int dim, int spacedim>
1154 void
1156 const std::vector<std::reference_wrapper<const Function<dim>>>
1157 & level_sets,
1158 const BoundingBox<dim> &box)
1159 {
1160 const Point<dim> center = box.center();
1161 const Definiteness definiteness =
1162 pointwise_definiteness(level_sets, center);
1163
1164 ExtendableQuadrature<dim> &quadrature =
1165 this->q_partitioning.quadrature_by_definiteness(definiteness);
1166
1167 quadrature.push_back(center, box.volume());
1168 }
1169
1170
1171
1172 template <int dim, int spacedim>
1173 void
1175 {
1176 AssertIndexRange(q_index, this->q_collection1D->size());
1177
1178 this->q_index = q_index;
1179 low_dim_algorithm.set_1D_quadrature(q_index);
1180 up_through_dimension_creator.set_1D_quadrature(q_index);
1181 }
1182
1183
1184
1185 template <int spacedim>
1187 const hp::QCollection<1> & q_collection1D,
1188 const AdditionalQGeneratorData &additional_data)
1189 : QGeneratorBase<1, spacedim>(q_collection1D, additional_data)
1190 , root_finder(
1191 RootFinder::AdditionalData(additional_data.root_finder_tolerance,
1192 additional_data.max_root_finder_splits))
1193 {
1194 Assert(q_collection1D.size() > 0,
1195 ExcMessage("Incoming quadrature collection is empty."));
1196 }
1197
1198
1199
1200 template <int spacedim>
1201 void
1203 const std::vector<std::reference_wrapper<const Function<1>>>
1204 & level_sets,
1205 const BoundingBox<1> &box,
1206 const unsigned int n_box_splits)
1207 {
1208 (void)n_box_splits;
1209
1210 roots.clear();
1211 root_finder.find_roots(level_sets, box, roots);
1212
1213 const Quadrature<1> &quadrature1D =
1214 (*this->q_collection1D)[this->q_index];
1215
1217 box,
1218 roots,
1219 zero_dim_point,
1220 unit_weight,
1221 direction,
1222 level_sets,
1223 this->additional_data,
1224 this->q_partitioning);
1225
1226 if (spacedim == 1)
1227 this->create_surface_points(level_sets);
1228 }
1229
1230
1231
1232 template <int spacedim>
1233 void
1235 const std::vector<std::reference_wrapper<const Function<1>>>
1236 &level_sets)
1237 {
1238 Assert(level_sets.size() == 1, ExcInternalError());
1239
1240 for (const double root : roots)
1241 {
1242 // A surface integral in 1d is just a point evaluation,
1243 // so the weight is always 1.
1244 const double weight = 1;
1245 const Point<1> point(root);
1246
1247 Tensor<1, 1> normal = level_sets[0].get().gradient(point);
1248 const double gradient_norm = normal.norm();
1249 Assert(
1250 gradient_norm > 1e-11,
1251 ExcMessage(
1252 "The level set function has a gradient almost equal to 0."));
1253 normal *= 1. / gradient_norm;
1254
1255 this->q_partitioning.surface.push_back(point, weight, normal);
1256 }
1257 }
1258
1259
1260
1261 template <int spacedim>
1262 void
1264 {
1265 AssertIndexRange(q_index, this->q_collection1D->size());
1266 this->q_index = q_index;
1267 }
1268
1269
1270
1271 } // namespace QuadratureGeneratorImplementation
1272
1273
1274
1275 namespace DiscreteQuadratureGeneratorImplementation
1276 {
1279 "The set_active_cell function has to be called before calling this function.");
1280
1283 "The reference cell of the incoming cell must be a hypercube.");
1284
1285
1311 template <int dim, class VectorType = Vector<double>>
1313 {
1314 public:
1323 const VectorType & dof_values);
1324
1328 void
1330 &cell) override;
1331
1339 double
1340 value(const Point<dim> & point,
1341 const unsigned int component = 0) const override;
1342
1351 gradient(const Point<dim> & point,
1352 const unsigned int component = 0) const override;
1353
1362 hessian(const Point<dim> & point,
1363 const unsigned int component = 0) const override;
1364
1365 private:
1369 bool
1370 cell_is_set() const;
1371
1376
1382
1388
1392 std::vector<types::global_dof_index> local_dof_indices;
1393
1398 std::vector<typename VectorType::value_type> local_dof_values;
1399
1405 std::vector<Polynomials::Polynomial<double>> poly;
1406
1410 std::vector<unsigned int> renumber;
1411
1416 };
1417
1418
1419
1420 template <int dim, class VectorType>
1422 const DoFHandler<dim> &dof_handler,
1423 const VectorType & dof_values)
1424 : dof_handler(&dof_handler)
1425 , global_dof_values(&dof_values)
1426 {
1427 Assert(dof_handler.n_dofs() == dof_values.size(),
1428 ExcDimensionMismatch(dof_handler.n_dofs(), dof_values.size()));
1429 }
1430
1431
1432
1433 template <int dim, class VectorType>
1434 void
1436 const typename Triangulation<dim>::active_cell_iterator &cell)
1437 {
1438 Assert(
1439 &cell->get_triangulation() == &dof_handler->get_triangulation(),
1440 ExcMessage(
1441 "The incoming cell must belong to the triangulation associated with "
1442 "the DoFHandler passed to the constructor."));
1443
1444 const typename DoFHandler<dim>::active_cell_iterator dof_handler_cell(
1445 &dof_handler->get_triangulation(),
1446 cell->level(),
1447 cell->index(),
1448 dof_handler);
1449
1450 // Save the element and the local dof values, since this is what we need
1451 // to evaluate the function.
1452
1453 // Check if we can use the fast path. In case we have a different
1454 // element from the one used before we need to set up the data
1455 // structures again.
1456 if (element != &dof_handler_cell->get_fe())
1457 {
1458 poly.clear();
1459 element = &dof_handler_cell->get_fe();
1460
1461 if (element->n_base_elements() == 1 &&
1463 *element, 0))
1464 {
1466 shape_info;
1467
1468 shape_info.reinit(QMidpoint<1>(), *element, 0);
1469 renumber = shape_info.lexicographic_numbering;
1470 poly =
1472 element->base_element(0));
1473
1474 polynomials_are_hat_functions =
1475 (poly.size() == 2 && poly[0].value(0.) == 1. &&
1476 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1477 poly[1].value(1.) == 1.);
1478 }
1479 }
1480 else
1481 element = &dof_handler_cell->get_fe();
1482
1483 local_dof_indices.resize(element->dofs_per_cell);
1484 dof_handler_cell->get_dof_indices(local_dof_indices);
1485
1486 local_dof_values.resize(element->dofs_per_cell);
1487
1488 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1489 local_dof_values[i] =
1491 *global_dof_values, local_dof_indices[i]);
1492 }
1493
1494
1495
1496 template <int dim, class VectorType>
1497 bool
1499 {
1500 // If set cell hasn't been called the size of local_dof_values will be
1501 // zero.
1502 return local_dof_values.size() > 0;
1503 }
1504
1505
1506
1507 template <int dim, class VectorType>
1508 double
1510 const Point<dim> & point,
1511 const unsigned int component) const
1512 {
1513 AssertIndexRange(component, this->n_components);
1514 Assert(cell_is_set(), ExcCellNotSet());
1515
1516 if (!poly.empty() && component == 0)
1517 {
1518 // TODO: this could be extended to a component that is not zero
1519 return ::internal::evaluate_tensor_product_value(
1520 poly,
1521 local_dof_values,
1522 point,
1523 polynomials_are_hat_functions,
1524 renumber);
1525 }
1526 else
1527 {
1528 double value = 0;
1529 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1530 value += local_dof_values[i] *
1531 element->shape_value_component(i, point, component);
1532
1533 return value;
1534 }
1535 }
1536
1537
1538
1539 template <int dim, class VectorType>
1542 const Point<dim> & point,
1543 const unsigned int component) const
1544 {
1545 AssertIndexRange(component, this->n_components);
1546 Assert(cell_is_set(), ExcCellNotSet());
1547
1548 if (!poly.empty() && component == 0)
1549 {
1550 // TODO: this could be extended to a component that is not zero
1551 return ::internal::evaluate_tensor_product_value_and_gradient(
1552 poly,
1553 local_dof_values,
1554 point,
1555 polynomials_are_hat_functions,
1556 renumber)
1557 .second;
1558 }
1559 else
1560 {
1561 Tensor<1, dim> gradient;
1562 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1563 gradient += local_dof_values[i] *
1564 element->shape_grad_component(i, point, component);
1565
1566 return gradient;
1567 }
1568 }
1569
1570
1571
1572 template <int dim, class VectorType>
1575 const Point<dim> & point,
1576 const unsigned int component) const
1577 {
1578 AssertIndexRange(component, this->n_components);
1579 Assert(cell_is_set(), ExcCellNotSet());
1580
1581 if (!poly.empty() && component == 0)
1582 {
1583 // TODO: this could be extended to a component that is not zero
1584 return ::internal::evaluate_tensor_product_hessian(
1585 poly, local_dof_values, point, renumber);
1586 }
1587 else
1588 {
1589 Tensor<2, dim> hessian;
1590 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1591 hessian +=
1592 local_dof_values[i] *
1593 element->shape_grad_grad_component(i, point, component);
1594
1595 return symmetrize(hessian);
1596 }
1597 }
1598 } // namespace DiscreteQuadratureGeneratorImplementation
1599 } // namespace internal
1600
1601
1602
1604 const unsigned int max_box_splits,
1605 const double lower_bound_implicit_function,
1606 const double min_distance_between_roots,
1607 const double limit_to_be_definite,
1608 const double root_finder_tolerance,
1609 const unsigned int max_root_finder_splits,
1610 bool split_in_half)
1611 : max_box_splits(max_box_splits)
1612 , lower_bound_implicit_function(lower_bound_implicit_function)
1613 , min_distance_between_roots(min_distance_between_roots)
1614 , limit_to_be_definite(limit_to_be_definite)
1615 , root_finder_tolerance(root_finder_tolerance)
1616 , max_root_finder_splits(max_root_finder_splits)
1617 , split_in_half(split_in_half)
1618 {}
1619
1620
1621
1622 template <int dim>
1624 const hp::QCollection<1> &q_collection,
1625 const AdditionalData & additional_data)
1626 : q_generator(q_collection, additional_data)
1627 {
1628 Assert(q_collection.size() > 0,
1629 ExcMessage("Incoming hp::QCollection<1> is empty."));
1630 }
1631
1632
1633
1634 template <int dim>
1635 void
1637 const BoundingBox<dim> &box)
1638 {
1639 Assert(level_set.n_components == 1,
1640 ExcMessage(
1641 "The incoming function should be a scalar level set function,"
1642 " it should have one component."));
1643 Assert(box.volume() > 0, ExcMessage("Incoming box has zero volume."));
1644
1645 q_generator.clear_quadratures();
1646
1647 std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
1648 level_sets.push_back(level_set);
1649
1650 const unsigned int n_box_splits = 0;
1651 q_generator.generate(level_sets, box, n_box_splits);
1652
1653 // With a single level set function, the "indefinite" quadrature should be
1654 // zero. If you call generate() with a ZeroFunction nothing good can be
1655 // done. You will end up here.
1656 Assert(
1657 q_generator.get_quadratures().indefinite.size() == 0,
1658 ExcMessage(
1659 "Generation of quadrature rules failed. This can mean that the level "
1660 "set function is degenerate in some way, e.g. oscillating extremely "
1661 "rapidly."));
1662 }
1663
1664
1665
1666 template <int dim>
1667 const Quadrature<dim> &
1669 {
1670 return q_generator.get_quadratures().negative;
1671 }
1672
1673
1674
1675 template <int dim>
1676 const Quadrature<dim> &
1678 {
1679 return q_generator.get_quadratures().positive;
1680 }
1681
1682
1683
1684 template <int dim>
1687 {
1688 return q_generator.get_quadratures().surface;
1689 }
1690
1691
1692 template <int dim>
1693 void
1695 {
1696 q_generator.set_1D_quadrature(q_index);
1697 }
1698
1699
1700
1701 template <int dim>
1703 const hp::QCollection<1> &quadratures1D,
1704 const AdditionalData & additional_data)
1705 : quadrature_generator(quadratures1D, additional_data)
1706 {}
1707
1708
1709
1710 template <int dim>
1711 void
1713 const BoundingBox<dim> &box,
1714 const unsigned int face_index)
1715 {
1717
1718 // We restrict the level set function to the face, by locking the coordinate
1719 // that is constant over the face. This will be the same as the direction of
1720 // the face normal.
1721 const unsigned int face_normal_direction =
1723
1724 const Point<dim> vertex0 =
1726 const double coordinate_value = vertex0(face_normal_direction);
1727
1728 const Functions::CoordinateRestriction<dim - 1> face_restriction(
1729 level_set, face_normal_direction, coordinate_value);
1730
1731 // Reuse the lower dimensional QuadratureGenerator on the face.
1732 const BoundingBox<dim - 1> cross_section =
1733 box.cross_section(face_normal_direction);
1734 quadrature_generator.generate(face_restriction, cross_section);
1735
1736 // We need the dim-dimensional normals of the zero-contour.
1737 // Recompute these.
1738 const ImmersedSurfaceQuadrature<dim - 1, dim - 1>
1739 &surface_quadrature_wrong_normal =
1740 quadrature_generator.get_surface_quadrature();
1741
1742 std::vector<Tensor<1, dim>> normals;
1743 normals.reserve(surface_quadrature_wrong_normal.size());
1744 for (unsigned int i = 0; i < surface_quadrature_wrong_normal.size(); ++i)
1745 {
1747 surface_quadrature_wrong_normal.point(i),
1748 face_normal_direction,
1749 coordinate_value);
1750
1751 Tensor<1, dim> normal = level_set.gradient(point);
1752 normal /= normal.norm();
1753 normals.push_back(normal);
1754 }
1755 surface_quadrature = ImmersedSurfaceQuadrature<dim - 1, dim>(
1756 surface_quadrature_wrong_normal.get_points(),
1757 surface_quadrature_wrong_normal.get_weights(),
1758 normals);
1759 }
1760
1761
1762
1763 template <int dim>
1764 void
1766 {
1767 quadrature_generator.set_1D_quadrature(q_index);
1768 }
1769
1770
1771
1772 template <int dim>
1773 const Quadrature<dim - 1> &
1775 {
1776 return quadrature_generator.get_inside_quadrature();
1777 }
1778
1779
1780 template <int dim>
1781 const Quadrature<dim - 1> &
1783 {
1784 return quadrature_generator.get_outside_quadrature();
1785 }
1786
1787
1788
1789 template <int dim>
1790 const ImmersedSurfaceQuadrature<dim - 1, dim> &
1792 {
1793 return surface_quadrature;
1794 }
1795
1796
1797
1799 const hp::QCollection<1> &quadratures1D,
1800 const AdditionalData & additional_data)
1801 {
1802 (void)quadratures1D;
1803 (void)additional_data;
1804 }
1805
1806
1807
1808 void
1810 const BoundingBox<1> &box,
1811 const unsigned int face_index)
1812 {
1814
1815 // The only vertex the 1d-face has.
1816 const Point<1> vertex =
1818
1819 const unsigned int n_points = 1;
1820 const double weight = 1;
1821 const std::vector<Point<0>> points(n_points);
1822 const std::vector<double> weights(n_points, weight);
1823
1824 const double level_set_value = level_set.value(vertex);
1825 if (level_set_value < 0)
1826 {
1827 inside_quadrature = Quadrature<0>(points, weights);
1828 outside_quadrature = Quadrature<0>();
1829 }
1830 else if (level_set_value > 0)
1831 {
1832 inside_quadrature = Quadrature<0>();
1833 outside_quadrature = Quadrature<0>(points, weights);
1834 }
1835 else
1836 {
1837 inside_quadrature = Quadrature<0>();
1838 outside_quadrature = Quadrature<0>();
1839 }
1840 }
1841
1842
1843
1844 void
1846 {
1847 (void)q_index;
1848 }
1849
1850
1851
1852 const Quadrature<0> &
1854 {
1855 return inside_quadrature;
1856 }
1857
1858
1859 const Quadrature<0> &
1861 {
1862 return outside_quadrature;
1863 }
1864
1865
1866
1869 {
1870 return surface_quadrature;
1871 }
1872
1873
1874
1875 template <int dim>
1876 template <class VectorType>
1878 const hp::QCollection<1> &quadratures1D,
1879 const DoFHandler<dim> & dof_handler,
1880 const VectorType & level_set,
1881 const AdditionalData & additional_data)
1882 : QuadratureGenerator<dim>(quadratures1D, additional_data)
1883 , reference_space_level_set(
1884 std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
1885 RefSpaceFEFieldFunction<dim, VectorType>>(
1886 dof_handler,
1887 level_set))
1888 {}
1889
1890
1891
1892 template <int dim>
1893 void
1895 const typename Triangulation<dim>::active_cell_iterator &cell)
1896 {
1897 Assert(cell->reference_cell().is_hyper_cube(),
1898 internal::DiscreteQuadratureGeneratorImplementation::
1899 ExcReferenceCellNotHypercube());
1900
1901 reference_space_level_set->set_active_cell(cell);
1902 const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
1903 QuadratureGenerator<dim>::generate(*reference_space_level_set, unit_box);
1904 }
1905
1906
1907
1908 template <int dim>
1909 template <class VectorType>
1911 const hp::QCollection<1> &quadratures1D,
1912 const DoFHandler<dim> & dof_handler,
1913 const VectorType & level_set,
1914 const AdditionalData & additional_data)
1915 : FaceQuadratureGenerator<dim>(quadratures1D, additional_data)
1916 , reference_space_level_set(
1917 std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
1918 RefSpaceFEFieldFunction<dim, VectorType>>(
1919 dof_handler,
1920 level_set))
1921 {}
1922
1923
1924
1925 template <int dim>
1926 void
1928 const typename Triangulation<dim>::active_cell_iterator &cell,
1929 const unsigned int face_index)
1930 {
1931 Assert(cell->reference_cell().is_hyper_cube(),
1932 internal::DiscreteQuadratureGeneratorImplementation::
1933 ExcReferenceCellNotHypercube());
1934
1935 reference_space_level_set->set_active_cell(cell);
1936 const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
1937 FaceQuadratureGenerator<dim>::generate(*reference_space_level_set,
1938 unit_box,
1939 face_index);
1940 }
1941} // namespace NonMatching
1942#include "quadrature_generator.inst"
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
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) 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
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(const typename Triangulation< dim >::active_cell_iterator &cell)
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)
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 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
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 use_midpoint_method(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box)
void split_box_and_recurse(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const std_cxx17::optional< HeightDirectionData > &direction_data, const unsigned int n_box_splits)
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:112
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:265
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 3 > center
Point< 2 > first
Definition grid_out.cc:4615
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
unsigned int compute_split_direction(const BoundingBox< dim > &box, const std_cxx17::optional< HeightDirectionData > &height_direction_data)
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)
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)
std_cxx17::optional< HeightDirectionData > find_best_height_direction(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)
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)
std_cxx17::optional< unsigned int > direction_of_largest_extent(const BoundingBox< dim > &box)
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)
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)
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)
Point< dim+1 > create_higher_dim_point(const Point< dim > &point, const unsigned int component_in_dim_plus_1, const double coordinate_value)
static const unsigned int invalid_unsigned_int
Definition types.h:213
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:435
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)