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