Reference documentation for deal.II version 9.4.1
Searching...
No Matches
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2021 - 2022 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
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 {
50 {
51 template <int dim>
52 void
54 const double weight,
56 const double start,
57 const double end,
58 const unsigned int component_in_dim,
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
84 const double start,
85 const double end,
86 const unsigned int component_in_dim,
88 {
89 for (unsigned int j = 0; j < lower.size(); ++j)
90 {
92 lower.weight(j),
94 start,
95 end,
96 component_in_dim,
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>(
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.
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
289 {
291 for (unsigned int d = 0; d < dim; ++d)
292 {
295 }
296 }
297 else
298 {
299 for (unsigned int d = 0; d < dim; ++d)
300 {
304 }
305 }
306 }
307 }
308
310 {
311 const auto max_element =
314
316 data.direction =
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
365 const BoundingBox<dim> & box,
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
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
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,
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.
491 {
492 // All points on the interval belong to the same region in
494 // the points to.
495 const Point<1> center(start + 0.5 * length);
496 const Definiteness definiteness =
497 pointwise_definiteness(level_sets, center);
500
502 weight,
504 start,
505 end,
506 height_function_direction,
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
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
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;
597 FunctionTools::taylor_estimate_function_bounds<1>(function,
598 interval,
599 value_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.
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>
641 {}
642
643
644
645 template <int dim>
646 void
648 const double weight)
649 {
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,
712 : q_collection1D(&q_collection1D)
714 , root_finder(
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,
730 const unsigned int height_function_direction,
731 QPartitioning<dim> & q_partitioning)
732 {
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,
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,
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
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 =
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,
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,
863 {
864 for (unsigned int i = 0; i < q_collection1D.size(); ++i)
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
903 {
905 box,
906 this->q_partitioning.positive);
907 }
908 else if (extreme_values.second <
910 {
912 box,
913 this->q_partitioning.negative);
914 }
915 else if (one_positive_one_negative_definite(all_function_bounds))
916 {
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 >
930 {
932 level_sets,
933 box,
934 n_box_splits);
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 {
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
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 =
1126
1128 (*this->q_collection1D)[this->q_index];
1129
1132 box.lower_bound(height_function_direction),
1133 box.upper_bound(height_function_direction),
1134 height_function_direction,
1135 this->q_partitioning.negative);
1136
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,
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
1166
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;
1181 }
1182
1183
1184
1185 template <int spacedim>
1187 const hp::QCollection<1> & q_collection1D,
1190 , root_finder(
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
1214 (*this->q_collection1D)[this->q_index];
1215
1217 box,
1218 roots,
1219 zero_dim_point,
1220 unit_weight,
1221 direction,
1222 level_sets,
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(
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
1272
1273
1274
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
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
1520 poly,
1521 local_dof_values,
1522 point,
1523 polynomials_are_hat_functions,
1524 renumber)
1525 .first;
1526 }
1527 else
1528 {
1529 double value = 0;
1530 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1531 value += local_dof_values[i] *
1532 element->shape_value_component(i, point, component);
1533
1534 return value;
1535 }
1536 }
1537
1538
1539
1540 template <int dim, class VectorType>
1543 const Point<dim> & point,
1544 const unsigned int component) const
1545 {
1546 AssertIndexRange(component, this->n_components);
1547 Assert(cell_is_set(), ExcCellNotSet());
1548
1549 if (!poly.empty() && component == 0)
1550 {
1551 // TODO: this could be extended to a component that is not zero
1553 poly,
1554 local_dof_values,
1555 point,
1556 polynomials_are_hat_functions,
1557 renumber)
1558 .second;
1559 }
1560 else
1561 {
1563 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1566
1568 }
1569 }
1570
1571
1572
1573 template <int dim, class VectorType>
1576 const Point<dim> & point,
1577 const unsigned int component) const
1578 {
1579 AssertIndexRange(component, this->n_components);
1580 Assert(cell_is_set(), ExcCellNotSet());
1581
1582 if (!poly.empty() && component == 0)
1583 {
1584 // TODO: this could be extended to a component that is not zero
1585 return ::internal::evaluate_tensor_product_hessian(
1586 poly, local_dof_values, point, renumber);
1587 }
1588 else
1589 {
1590 Tensor<2, dim> hessian;
1591 for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1592 hessian +=
1593 local_dof_values[i] *
1595
1596 return symmetrize(hessian);
1597 }
1598 }
1600 } // namespace internal
1601
1602
1603
1605 const unsigned int max_box_splits,
1606 const double lower_bound_implicit_function,
1607 const double min_distance_between_roots,
1608 const double limit_to_be_definite,
1609 const double root_finder_tolerance,
1610 const unsigned int max_root_finder_splits,
1611 bool split_in_half)
1612 : max_box_splits(max_box_splits)
1613 , lower_bound_implicit_function(lower_bound_implicit_function)
1614 , min_distance_between_roots(min_distance_between_roots)
1615 , limit_to_be_definite(limit_to_be_definite)
1616 , root_finder_tolerance(root_finder_tolerance)
1617 , max_root_finder_splits(max_root_finder_splits)
1618 , split_in_half(split_in_half)
1619 {}
1620
1621
1622
1623 template <int dim>
1625 const hp::QCollection<1> &q_collection,
1628 {
1629 Assert(q_collection.size() > 0,
1630 ExcMessage("Incoming hp::QCollection<1> is empty."));
1631 }
1632
1633
1634
1635 template <int dim>
1636 void
1638 const BoundingBox<dim> &box)
1639 {
1640 Assert(level_set.n_components == 1,
1641 ExcMessage(
1642 "The incoming function should be a scalar level set function,"
1643 " it should have one component."));
1644 Assert(box.volume() > 0, ExcMessage("Incoming box has zero volume."));
1645
1647
1648 std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
1649 level_sets.push_back(level_set);
1650
1651 const unsigned int n_box_splits = 0;
1652 q_generator.generate(level_sets, box, n_box_splits);
1653
1654 // With a single level set function, the "indefinite" quadrature should be
1655 // zero. If you call generate() with a ZeroFunction nothing good can be
1656 // done. You will end up here.
1657 Assert(
1659 ExcMessage(
1660 "Generation of quadrature rules failed. This can mean that the level "
1661 "set function is degenerate in some way, e.g. oscillating extremely "
1662 "rapidly."));
1663 }
1664
1665
1666
1667 template <int dim>
1670 {
1672 }
1673
1674
1675
1676 template <int dim>
1679 {
1681 }
1682
1683
1684
1685 template <int dim>
1688 {
1690 }
1691
1692
1693 template <int dim>
1694 void
1696 {
1698 }
1699
1700
1701
1702 template <int dim>
1707 {}
1708
1709
1710
1711 template <int dim>
1712 void
1714 const BoundingBox<dim> &box,
1715 const unsigned int face_index)
1716 {
1718
1719 // We restrict the level set function to the face, by locking the coordinate
1720 // that is constant over the face. This will be the same as the direction of
1721 // the face normal.
1722 const unsigned int face_normal_direction =
1724
1725 const Point<dim> vertex0 =
1727 const double coordinate_value = vertex0(face_normal_direction);
1728
1729 const Functions::CoordinateRestriction<dim - 1> face_restriction(
1730 level_set, face_normal_direction, coordinate_value);
1731
1732 // Reuse the lower dimensional QuadratureGenerator on the face.
1733 const BoundingBox<dim - 1> cross_section =
1734 box.cross_section(face_normal_direction);
1736
1737 // We need the dim-dimensional normals of the zero-contour.
1738 // Recompute these.
1739 const ImmersedSurfaceQuadrature<dim - 1, dim - 1>
1742
1743 std::vector<Tensor<1, dim>> normals;
1745 for (unsigned int i = 0; i < surface_quadrature_wrong_normal.size(); ++i)
1746 {
1749 face_normal_direction,
1750 coordinate_value);
1751
1752 Tensor<1, dim> normal = level_set.gradient(point);
1753 normal /= normal.norm();
1754 normals.push_back(normal);
1755 }
1759 normals);
1760 }
1761
1762
1763
1764 template <int dim>
1765 void
1767 {
1769 }
1770
1771
1772
1773 template <int dim>
1774 const Quadrature<dim - 1> &
1776 {
1778 }
1779
1780
1781 template <int dim>
1782 const Quadrature<dim - 1> &
1784 {
1786 }
1787
1788
1789
1790 template <int dim>
1791 const ImmersedSurfaceQuadrature<dim - 1, dim> &
1793 {
1795 }
1796
1797
1798
1802 {
1805 }
1806
1807
1808
1809 void
1811 const BoundingBox<1> &box,
1812 const unsigned int face_index)
1813 {
1815
1816 // The only vertex the 1D-face has.
1817 const Point<1> vertex =
1819
1820 const unsigned int n_points = 1;
1821 const double weight = 1;
1822 const std::vector<Point<0>> points(n_points);
1823 const std::vector<double> weights(n_points, weight);
1824
1825 const double level_set_value = level_set.value(vertex);
1826 if (level_set_value < 0)
1827 {
1830 }
1831 else if (level_set_value > 0)
1832 {
1835 }
1836 else
1837 {
1840 }
1841 }
1842
1843
1844
1845 void
1847 {
1848 (void)q_index;
1849 }
1850
1851
1852
1855 {
1857 }
1858
1859
1862 {
1864 }
1865
1866
1867
1870 {
1872 }
1873
1874
1875
1876 template <int dim>
1877 template <class VectorType>
1880 const DoFHandler<dim> & dof_handler,
1881 const VectorType & level_set,
1884 , reference_space_level_set(
1886 RefSpaceFEFieldFunction<dim, VectorType>>(
1887 dof_handler,
1888 level_set))
1889 {}
1890
1891
1892
1893 template <int dim>
1894 void
1896 const typename Triangulation<dim>::active_cell_iterator &cell)
1897 {
1898 Assert(cell->reference_cell().is_hyper_cube(),
1900 ExcReferenceCellNotHypercube());
1901
1902 reference_space_level_set->set_active_cell(cell);
1903 const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
1905 }
1906
1907
1908
1909 template <int dim>
1910 template <class VectorType>
1913 const DoFHandler<dim> & dof_handler,
1914 const VectorType & level_set,
1917 , reference_space_level_set(
1919 RefSpaceFEFieldFunction<dim, VectorType>>(
1920 dof_handler,
1921 level_set))
1922 {}
1923
1924
1925
1926 template <int dim>
1927 void
1929 const typename Triangulation<dim>::active_cell_iterator &cell,
1930 const unsigned int face_index)
1931 {
1932 Assert(cell->reference_cell().is_hyper_cube(),
1934 ExcReferenceCellNotHypercube());
1935
1936 reference_space_level_set->set_active_cell(cell);
1937 const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
1939 unit_box,
1940 face_index);
1941 }
1942} // namespace NonMatching
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 unsigned int 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
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)
void generate(const Function< dim > &level_set, const BoundingBox< dim > &box, const unsigned int face_index)
void push_back(const Point< dim > &point, const double weight, const Tensor< 1, spacedim > &normal)
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
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)
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)
Definition: point.h:111
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
unsigned int n_vertices() const
Definition: tensor.h:503
numbers::NumberTraits< Number >::real_type norm() const
unsigned int size() const
Definition: collection.h:264
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > center
Point< 2 > first
Definition: grid_out.cc:4603
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
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
Definition: dof_handler.h:438
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)
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)
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:201
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:417
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)