Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
quadrature_generator.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 
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  {
51  namespace QuadratureGeneratorImplementation
52  {
53  template <int dim>
54  void
56  const double weight,
57  const Quadrature<1> &quadrature1D,
58  const double start,
59  const double end,
60  const unsigned int component_in_dim,
61  ExtendableQuadrature<dim> &quadrature)
62  {
63  Assert(start < end,
64  ExcMessage("Interval start must be less than interval end."));
65 
66  const double length = end - start;
67  for (unsigned int j = 0; j < quadrature1D.size(); ++j)
68  {
69  const double x = start + length * quadrature1D.point(j)[0];
71  point, component_in_dim, x),
72  length * weight * quadrature1D.weight(j));
73  }
74  }
75 
76 
77 
82  template <int dim>
83  void
85  const Quadrature<1> &quadrature1D,
86  const double start,
87  const double end,
88  const unsigned int component_in_dim,
89  ExtendableQuadrature<dim> &quadrature)
90  {
91  for (unsigned int j = 0; j < lower.size(); ++j)
92  {
94  lower.weight(j),
95  quadrature1D,
96  start,
97  end,
98  component_in_dim,
99  quadrature);
100  }
101  }
102 
103 
104 
105  template <int dim>
108  const std::vector<std::reference_wrapper<const Function<dim>>>
109  &functions,
110  const Point<dim> &point)
111  {
112  Assert(functions.size() > 0,
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>(
192  function, box, bounds.value, bounds.gradient);
193  take_min_max_at_vertices(function, box, bounds.value);
194 
195  all_function_bounds.push_back(bounds);
196  }
197  }
198 
199 
200 
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."));
260 
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.
282  std::optional<std::array<double, dim>> min_lower_abs_grad;
283 
284  for (const FunctionBounds<dim> &bounds : all_function_bounds)
285  {
286  if (is_indefinite(bounds.value))
287  {
288  // For the first indefinite function we find, we write the lower
289  // bounds on each gradient component to min_lower_abs_grad.
290  if (!min_lower_abs_grad)
291  {
292  min_lower_abs_grad.emplace();
293  for (unsigned int d = 0; d < dim; ++d)
294  {
295  (*min_lower_abs_grad)[d] =
296  lower_bound_on_abs(bounds.gradient[d]);
297  }
298  }
299  else
300  {
301  for (unsigned int d = 0; d < dim; ++d)
302  {
303  (*min_lower_abs_grad)[d] =
304  std::min((*min_lower_abs_grad)[d],
305  lower_bound_on_abs(bounds.gradient[d]));
306  }
307  }
308  }
309  }
310 
311  if (min_lower_abs_grad)
312  {
313  const auto max_element =
314  std::max_element(min_lower_abs_grad->begin(),
315  min_lower_abs_grad->end());
316 
317  HeightDirectionData data;
318  data.direction =
319  std::distance(min_lower_abs_grad->begin(), max_element);
320  data.min_abs_dfdx = *max_element;
321 
322  return data;
323  }
324 
325  return std::optional<HeightDirectionData>();
326  }
327 
328 
329 
335  template <int dim>
336  inline bool
338  const std::vector<FunctionBounds<dim>> &all_function_bounds)
339  {
340  if (all_function_bounds.size() != 2)
341  return false;
342  else
343  {
344  const FunctionBounds<dim> &bounds0 = all_function_bounds.at(0);
345  const FunctionBounds<dim> &bounds1 = all_function_bounds.at(1);
346 
347  if (bounds0.value.first > 0 && bounds1.value.second < 0)
348  return true;
349  if (bounds1.value.first > 0 && bounds0.value.second < 0)
350  return true;
351  return false;
352  }
353  }
354 
355 
356 
364  template <int dim>
365  void
366  map_quadrature_to_box(const Quadrature<dim> &unit_quadrature,
367  const BoundingBox<dim> &box,
368  ExtendableQuadrature<dim> &quadrature)
369  {
370  for (unsigned int i = 0; i < unit_quadrature.size(); ++i)
371  {
372  const Point<dim> point = box.unit_to_real(unit_quadrature.point(i));
373  const double weight = unit_quadrature.weight(i) * box.volume();
374 
375  quadrature.push_back(point, weight);
376  }
377  }
378 
379 
380 
391  template <int dim>
392  void
394  const std::vector<std::reference_wrapper<const Function<dim>>>
395  &functions,
396  const BoundingBox<dim> &box,
397  const unsigned int direction,
398  std::vector<Functions::CoordinateRestriction<dim - 1>> &restrictions)
399  {
400  AssertIndexRange(direction, dim);
401 
402  restrictions.clear();
403  restrictions.reserve(2 * functions.size());
404 
405  const double bottom = box.lower_bound(direction);
406  const double top = box.upper_bound(direction);
407 
408  for (const auto &function : functions)
409  {
410  restrictions.push_back(Functions::CoordinateRestriction<dim - 1>(
411  function, direction, bottom));
412  restrictions.push_back(Functions::CoordinateRestriction<dim - 1>(
413  function, direction, top));
414  }
415  }
416 
417 
418 
427  template <int dim>
428  void
430  const std::vector<std::reference_wrapper<const Function<dim>>>
431  &functions,
432  const Point<dim - 1> &point,
433  const unsigned int open_direction,
434  std::vector<Functions::PointRestriction<dim - 1>> &restrictions)
435  {
436  AssertIndexRange(open_direction, dim);
437 
438  restrictions.clear();
439  restrictions.reserve(functions.size());
440  for (const auto &function : functions)
441  {
442  restrictions.push_back(Functions::PointRestriction<dim - 1>(
443  function, open_direction, point));
444  }
445  }
446 
447 
448 
462  template <int dim>
463  void
465  const Quadrature<1> &quadrature1D,
466  const BoundingBox<1> &interval,
467  const std::vector<double> &roots,
468  const Point<dim - 1> &point,
469  const double weight,
470  const unsigned int height_function_direction,
471  const std::vector<std::reference_wrapper<const Function<1>>>
472  &level_sets,
473  const AdditionalQGeneratorData &additional_data,
474  QPartitioning<dim> &q_partitioning)
475  {
476  // Make this int to avoid a warning signed/unsigned comparison.
477  const int n_roots = roots.size();
478 
479  // The number of intervals are roots.size() + 1
480  for (int i = -1; i < n_roots; ++i)
481  {
482  // Start and end point of the subinterval.
483  const double start = i < 0 ? interval.lower_bound(0) : roots[i];
484  const double end =
485  i + 1 < n_roots ? roots[i + 1] : interval.upper_bound(0);
486 
487  const double length = end - start;
488  // It might be that the end points of the subinterval are roots.
489  // If this is the case then the subinterval has length zero.
490  // Don't distribute points on the subinterval if it is shorter than
491  // some tolerance.
492  if (length > additional_data.min_distance_between_roots)
493  {
494  // All points on the interval belong to the same region in
495  // the QPartitioning. Determine the quadrature we should add
496  // the points to.
497  const Point<1> center(start + 0.5 * length);
498  const Definiteness definiteness =
499  pointwise_definiteness(level_sets, center);
500  ExtendableQuadrature<dim> &target_quadrature =
501  q_partitioning.quadrature_by_definiteness(definiteness);
502 
504  weight,
505  quadrature1D,
506  start,
507  end,
508  height_function_direction,
509  target_quadrature);
510  }
511  }
512  }
513 
514 
515 
517  const double tolerance,
518  const unsigned int max_recursion_depth,
519  const unsigned int max_iterations)
520  : tolerance(tolerance)
521  , max_recursion_depth(max_recursion_depth)
522  , max_iterations(max_iterations)
523  {}
524 
525 
526 
528  : additional_data(data)
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 
578  boost::uintmax_t iterations = additional_data.max_iterations;
579 
580  const std::pair<double, double> root_bracket =
581  boost::math::tools::toms748_solve(lambda,
582  interval.lower_bound(0),
583  interval.upper_bound(0),
584  left_value,
585  right_value,
586  stopping_criteria,
587  iterations);
588 
589  const double root = .5 * (root_bracket.first + root_bracket.second);
590  roots.push_back(root);
591  }
592  else
593  {
594  // Compute bounds on the incoming function to check if there are
595  // roots. If the function is positive or negative on the whole
596  // interval we do nothing.
597  std::pair<double, double> value_bounds;
598  std::array<std::pair<double, double>, 1> gradient_bounds;
599  FunctionTools::taylor_estimate_function_bounds<1>(function,
600  interval,
601  value_bounds,
602  gradient_bounds);
603 
604  // Since we already know the function values at the interval ends we
605  // might as well check these for min/max too.
606  const double function_min =
607  std::min(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.
623  if (recursion_depth < additional_data.max_recursion_depth)
624  {
625  find_roots(function,
626  interval.child(0),
627  recursion_depth + 1,
628  roots);
629  find_roots(function,
630  interval.child(1),
631  recursion_depth + 1,
632  roots);
633  }
634  }
635  }
636 
637 
638 
639  template <int dim>
641  const Quadrature<dim> &quadrature)
642  : Quadrature<dim>(quadrature)
643  {}
644 
645 
646 
647  template <int dim>
648  inline void
650  {
651  this->quadrature_points.clear();
652  this->weights.clear();
653  }
654 
655 
656 
657  template <int dim>
658  void
660  const double weight)
661  {
662  this->quadrature_points.push_back(point);
663  this->weights.push_back(weight);
664  }
665 
666 
667 
668  template <int dim>
671  const Definiteness definiteness)
672  {
673  switch (definiteness)
674  {
676  return negative;
678  return positive;
679  default:
680  return indefinite;
681  }
682  }
683 
684 
685 
686  template <int dim>
687  void
689  {
690  positive.clear();
691  negative.clear();
692  indefinite.clear();
693  surface.clear();
694  }
695 
696 
697 
705  template <int dim>
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,
735  const AdditionalQGeneratorData &additional_data)
736  : q_collection1D(&q_collection1D)
737  , additional_data(additional_data)
738  , root_finder(
739  RootFinder::AdditionalData(additional_data.root_finder_tolerance,
740  additional_data.max_root_finder_splits))
741  {
742  q_index = 0;
743  }
744 
745 
746 
747  template <int dim, int spacedim>
748  void
750  const std::vector<std::reference_wrapper<const Function<dim>>>
751  &level_sets,
752  const BoundingBox<dim> &box,
753  const Quadrature<dim - 1> &low_dim_quadrature,
754  const unsigned int height_function_direction,
755  QPartitioning<dim> &q_partitioning)
756  {
757  const Quadrature<1> &quadrature1D = (*q_collection1D)[q_index];
758 
759  for (unsigned int q = 0; q < low_dim_quadrature.size(); ++q)
760  {
761  const Point<dim - 1> &point = low_dim_quadrature.point(q);
762  const double weight = low_dim_quadrature.weight(q);
763  restrict_to_point(level_sets,
764  point,
765  height_function_direction,
766  point_restrictions);
767 
768  // We need a vector of references to do the recursive call.
769  const std::vector<std::reference_wrapper<const Function<1>>>
770  restrictions(point_restrictions.begin(),
771  point_restrictions.end());
772 
773  const BoundingBox<1> bounds_in_direction =
774  box.bounds(height_function_direction);
775 
776  roots.clear();
777  root_finder.find_roots(restrictions, bounds_in_direction, roots);
778 
779  distribute_points_between_roots(quadrature1D,
780  bounds_in_direction,
781  roots,
782  point,
783  weight,
784  height_function_direction,
785  restrictions,
786  additional_data,
787  q_partitioning);
788 
789  if (dim == spacedim)
790  create_surface_point(point,
791  weight,
792  level_sets,
793  box,
794  height_function_direction,
795  q_partitioning.surface);
796  }
797 
798  point_restrictions.clear();
799  }
800 
801 
802 
803  template <int dim, int spacedim>
804  void
806  const Point<dim - 1> &point,
807  const double weight,
808  const std::vector<std::reference_wrapper<const Function<dim>>>
809  &level_sets,
810  const BoundingBox<dim> &box,
811  const unsigned int height_function_direction,
812  ImmersedSurfaceQuadrature<dim> &surface_quadrature)
813  {
814  AssertIndexRange(roots.size(), 2);
815  Assert(level_sets.size() == 1, ExcInternalError());
816 
817 
818  const Function<dim> &level_set = level_sets.at(0);
819 
820  Point<dim> surface_point;
821  if (roots.size() == 1)
822  {
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 
843  const Tensor<1, dim> gradient = level_set.gradient(surface_point);
844  Tensor<1, dim> normal = gradient;
845  normal *= 1. / normal.norm();
846 
847  // Note that gradient[height_function_direction] is non-zero
848  // because of the implicit function theorem.
849  const double surface_weight =
850  weight * gradient.norm() /
851  std::abs(gradient[height_function_direction]);
852  surface_quadrature.push_back(surface_point, surface_weight, normal);
853  }
854 
855 
856 
857  template <int dim, int spacedim>
858  void
860  unsigned int q_index)
861  {
862  AssertIndexRange(q_index, q_collection1D->size());
863  this->q_index = q_index;
864  }
865 
866 
867 
868  template <int dim, int spacedim>
870  const hp::QCollection<1> &q_collection1D,
871  const AdditionalQGeneratorData &additional_data)
872  : additional_data(additional_data)
873  , q_collection1D(&q_collection1D)
874  {
875  q_index = 0;
876  }
877 
878 
879 
880  template <int dim, int spacedim>
882  const hp::QCollection<1> &q_collection1D,
883  const AdditionalQGeneratorData &additional_data)
884  : QGeneratorBase<dim, spacedim>(q_collection1D, additional_data)
885  , low_dim_algorithm(q_collection1D, additional_data)
886  , up_through_dimension_creator(q_collection1D, additional_data)
887  {
888  for (const auto &quadrature : q_collection1D)
889  tensor_products.push_back(quadrature);
890  }
891 
892 
893 
894  template <int dim, int spacedim>
895  void
897  {
898  q_partitioning.clear();
899  }
900 
901 
902 
903  template <int dim, int spacedim>
904  const QPartitioning<dim> &
906  {
907  return q_partitioning;
908  }
909 
910 
911 
912  template <int dim, int spacedim>
913  void
915  const std::vector<std::reference_wrapper<const Function<dim>>>
916  &level_sets,
917  const BoundingBox<dim> &box,
918  const unsigned int n_box_splits)
919  {
920  std::vector<FunctionBounds<dim>> all_function_bounds;
921  estimate_function_bounds(level_sets, box, all_function_bounds);
922 
923  const std::pair<double, double> extreme_values =
924  find_extreme_values(all_function_bounds);
925 
926  if (extreme_values.first > this->additional_data.limit_to_be_definite)
927  {
928  map_quadrature_to_box(tensor_products[this->q_index],
929  box,
930  this->q_partitioning.positive);
931  }
932  else if (extreme_values.second <
933  -(this->additional_data.limit_to_be_definite))
934  {
935  map_quadrature_to_box(tensor_products[this->q_index],
936  box,
937  this->q_partitioning.negative);
938  }
939  else if (one_positive_one_negative_definite(all_function_bounds))
940  {
941  map_quadrature_to_box(tensor_products[this->q_index],
942  box,
943  this->q_partitioning.indefinite);
944  }
945  else
946  {
947  const std::optional<HeightDirectionData> data =
948  find_best_height_direction(all_function_bounds);
949 
950  // Check larger than a constant to avoid that min_abs_dfdx is only
951  // larger by 0 by floating point precision.
952  if (data && data->min_abs_dfdx >
953  this->additional_data.lower_bound_implicit_function)
954  {
955  create_low_dim_quadratures(data->direction,
956  level_sets,
957  box,
958  n_box_splits);
959  create_high_dim_quadratures(data->direction, level_sets, box);
960  }
961  else if (n_box_splits < this->additional_data.max_box_splits)
962  {
963  split_box_and_recurse(level_sets, box, data, n_box_splits);
964  }
965  else
966  {
967  // We can't split the box recursively forever. Use the midpoint
968  // method as a last resort.
969  use_midpoint_method(level_sets, box);
970  }
971  }
972  }
973 
974 
975 
981  template <int dim>
982  std::optional<unsigned int>
984  {
985  // Get the side lengths for each direction and sort them.
986  std::array<std::pair<double, unsigned int>, dim> side_lengths;
987  for (unsigned int i = 0; i < dim; ++i)
988  {
989  side_lengths[i].first = box.side_length(i);
990  side_lengths[i].second = i;
991  }
992  // Sort is lexicographic, so this sorts based on side length first.
993  std::sort(side_lengths.begin(), side_lengths.end());
994 
995  // Check if the two largest side lengths have the same length. This
996  // function isn't called in 1d, so the (dim - 2)-element exists.
997  if (boost::math::epsilon_difference(side_lengths[dim - 1].first,
998  side_lengths[dim - 2].first) < 100)
999  return std::optional<unsigned int>();
1000 
1001  return side_lengths.back().second;
1002  }
1003 
1004 
1005 
1017  template <int dim>
1018  unsigned int
1020  const BoundingBox<dim> &box,
1021  const std::optional<HeightDirectionData> &height_direction_data)
1022  {
1023  const std::optional<unsigned int> direction =
1025 
1026  if (direction)
1027  return *direction;
1028 
1029  // This direction is closest to being a height direction, so
1030  // we split in this direction.
1031  if (height_direction_data)
1032  return height_direction_data->direction;
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  {
1087  if (this->additional_data.split_in_half)
1088  {
1089  const unsigned int direction =
1090  compute_split_direction(box, direction_data);
1091 
1092  const BoundingBox<dim> left_box = left_half(box, direction);
1093  const BoundingBox<dim> right_box = right_half(box, direction);
1094 
1095  generate(level_sets, left_box, n_box_splits + 1);
1096  generate(level_sets, right_box, n_box_splits + 1);
1097  }
1098  else
1099  {
1100  for (unsigned int i = 0;
1101  i < GeometryInfo<dim>::max_children_per_cell;
1102  ++i)
1103  {
1104  generate(level_sets, box.child(i), n_box_splits + 1);
1105  }
1106  }
1107  }
1108 
1109 
1110 
1111  template <int dim, int spacedim>
1112  void
1114  const unsigned int height_function_direction,
1115  const std::vector<std::reference_wrapper<const Function<dim>>>
1116  &level_sets,
1117  const BoundingBox<dim> &box,
1118  const unsigned int n_box_splits)
1119  {
1120  std::vector<Functions::CoordinateRestriction<dim - 1>>
1121  face_restrictions;
1122  restrict_to_top_and_bottom(level_sets,
1123  box,
1124  height_function_direction,
1125  face_restrictions);
1126 
1127  // We need a vector of references to do the recursive call.
1128  const std::vector<std::reference_wrapper<const Function<dim - 1>>>
1129  restrictions(face_restrictions.begin(), face_restrictions.end());
1130 
1131  const BoundingBox<dim - 1> cross_section =
1132  box.cross_section(height_function_direction);
1133 
1134  low_dim_algorithm.clear_quadratures();
1135  low_dim_algorithm.generate(restrictions, cross_section, n_box_splits);
1136  }
1137 
1138 
1139 
1140  template <int dim, int spacedim>
1141  void
1143  const unsigned int height_function_direction,
1144  const std::vector<std::reference_wrapper<const Function<dim>>>
1145  &level_sets,
1146  const BoundingBox<dim> &box)
1147  {
1148  const QPartitioning<dim - 1> &low_dim_quadratures =
1149  low_dim_algorithm.get_quadratures();
1150 
1151  const Quadrature<1> &quadrature1D =
1152  (*this->q_collection1D)[this->q_index];
1153 
1154  add_tensor_product(low_dim_quadratures.negative,
1155  quadrature1D,
1156  box.lower_bound(height_function_direction),
1157  box.upper_bound(height_function_direction),
1158  height_function_direction,
1159  this->q_partitioning.negative);
1161  add_tensor_product(low_dim_quadratures.positive,
1162  quadrature1D,
1163  box.lower_bound(height_function_direction),
1164  box.upper_bound(height_function_direction),
1165  height_function_direction,
1166  this->q_partitioning.positive);
1167 
1168  up_through_dimension_creator.generate(level_sets,
1169  box,
1170  low_dim_quadratures.indefinite,
1171  height_function_direction,
1172  this->q_partitioning);
1173  }
1175 
1176 
1177  template <int dim, int spacedim>
1178  void
1180  const std::vector<std::reference_wrapper<const Function<dim>>>
1181  &level_sets,
1182  const BoundingBox<dim> &box)
1183  {
1184  const Point<dim> center = box.center();
1185  const Definiteness definiteness =
1186  pointwise_definiteness(level_sets, center);
1187 
1188  ExtendableQuadrature<dim> &quadrature =
1189  this->q_partitioning.quadrature_by_definiteness(definiteness);
1190 
1191  quadrature.push_back(center, box.volume());
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;
1203  low_dim_algorithm.set_1D_quadrature(q_index);
1204  up_through_dimension_creator.set_1D_quadrature(q_index);
1205  }
1206 
1208 
1209  template <int spacedim>
1211  const hp::QCollection<1> &q_collection1D,
1212  const AdditionalQGeneratorData &additional_data)
1213  : QGeneratorBase<1, spacedim>(q_collection1D, additional_data)
1214  , root_finder(
1215  RootFinder::AdditionalData(additional_data.root_finder_tolerance,
1216  additional_data.max_root_finder_splits))
1217  {
1218  Assert(q_collection1D.size() > 0,
1219  ExcMessage("Incoming quadrature collection is empty."));
1220  }
1221 
1222 
1223 
1224  template <int spacedim>
1225  void
1227  const std::vector<std::reference_wrapper<const Function<1>>>
1228  &level_sets,
1229  const BoundingBox<1> &box,
1230  const unsigned int n_box_splits)
1231  {
1232  (void)n_box_splits;
1233 
1234  roots.clear();
1235  root_finder.find_roots(level_sets, box, roots);
1237  const Quadrature<1> &quadrature1D =
1238  (*this->q_collection1D)[this->q_index];
1239 
1240  distribute_points_between_roots(quadrature1D,
1241  box,
1242  roots,
1243  zero_dim_point,
1244  unit_weight,
1245  direction,
1246  level_sets,
1247  this->additional_data,
1248  this->q_partitioning);
1249 
1250  if (spacedim == 1)
1251  this->create_surface_points(level_sets);
1252  }
1253 
1254 
1255 
1256  template <int spacedim>
1257  void
1259  const std::vector<std::reference_wrapper<const Function<1>>>
1260  &level_sets)
1261  {
1262  Assert(level_sets.size() == 1, ExcInternalError());
1263 
1264  for (const double root : roots)
1265  {
1266  // A surface integral in 1d is just a point evaluation,
1267  // so the weight is always 1.
1268  const double weight = 1;
1269  const Point<1> point(root);
1270 
1271  Tensor<1, 1> normal = level_sets[0].get().gradient(point);
1272  const double gradient_norm = normal.norm();
1273  Assert(
1274  gradient_norm > 1e-11,
1275  ExcMessage(
1276  "The level set function has a gradient almost equal to 0."));
1277  normal *= 1. / gradient_norm;
1278 
1279  this->q_partitioning.surface.push_back(point, weight, normal);
1280  }
1281  }
1282 
1283 
1284 
1285  template <int spacedim>
1286  void
1288  {
1290  this->q_index = q_index;
1291  }
1292 
1293 
1294 
1295  } // namespace QuadratureGeneratorImplementation
1296 
1297 
1298 
1299  namespace DiscreteQuadratureGeneratorImplementation
1300  {
1302  ExcCellNotSet,
1303  "The set_active_cell function has to be called before calling this function.");
1304 
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
1360  set_subcell(const std::vector<unsigned int> &mask,
1361  const BoundingBox<dim> &subcell_box) override;
1362 
1366  bool
1367  is_fe_q_iso_q1() const override;
1368 
1372  unsigned int
1373  n_subdivisions() const override;
1374 
1382  double
1383  value(const Point<dim> &point,
1384  const unsigned int component = 0) const override;
1385 
1394  gradient(const Point<dim> &point,
1395  const unsigned int component = 0) const override;
1396 
1405  hessian(const Point<dim> &point,
1406  const unsigned int component = 0) const override;
1407 
1408  private:
1412  bool
1413  cell_is_set() const;
1414 
1419 
1425 
1431 
1435  std::vector<types::global_dof_index> local_dof_indices;
1436 
1441  std::vector<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
1572  const std::vector<unsigned int> &mask,
1573  const BoundingBox<dim> &subcell_box)
1574  {
1575  for (unsigned int i = 0; i < local_dof_values_subcell.size(); ++i)
1576  local_dof_values_subcell[i] = local_dof_values[renumber[mask[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  {
1670  Tensor<1, dim> gradient;
1671  for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
1672  gradient += local_dof_values[i] *
1673  element->shape_grad_component(i, point, component);
1674 
1675  return gradient;
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] *
1706  element->shape_grad_grad_component(i, point, component);
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.
1745  std::vector<unsigned int> mask;
1746  mask.reserve(1 << dim);
1747  mask.push_back(0);
1748  unsigned int stride = 1;
1749  for (unsigned int d = 0; d < dim; ++d)
1750  {
1751  const unsigned int previous_mask_size = 1U << d;
1752  AssertDimension(mask.size(), previous_mask_size);
1753  for (unsigned int i = 0; i < previous_mask_size; ++i)
1754  {
1755  mask[i] += subcell_indices[d] * stride;
1756  mask.push_back(stride + mask[i]);
1757  }
1758  stride *= dofs_per_line;
1759  }
1760  return mask;
1761  }
1762  } // namespace DiscreteQuadratureGeneratorImplementation
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,
1789  const AdditionalData &additional_data)
1790  : q_generator(q_collection, additional_data)
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  {
1802  q_generator.clear_quadratures();
1803  }
1804 
1805 
1806 
1807  template <int dim>
1808  void
1810  const BoundingBox<dim> &box)
1811  {
1812  clear_quadratures();
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(
1839  q_generator.get_quadratures().indefinite.empty(),
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>
1849  const Quadrature<dim> &
1851  {
1852  return q_generator.get_quadratures().negative;
1853  }
1854 
1855 
1856 
1857  template <int dim>
1858  const Quadrature<dim> &
1860  {
1861  return q_generator.get_quadratures().positive;
1862  }
1863 
1864 
1865 
1866  template <int dim>
1869  {
1870  return q_generator.get_quadratures().surface;
1871  }
1872 
1873 
1874  template <int dim>
1875  void
1877  {
1878  q_generator.set_1D_quadrature(q_index);
1879  }
1880 
1881 
1882 
1883  template <int dim>
1885  const hp::QCollection<1> &quadratures1D,
1886  const AdditionalData &additional_data)
1887  : quadrature_generator(quadratures1D, additional_data)
1888  {}
1889 
1890 
1891 
1892  template <int dim>
1893  void
1895  {
1896  quadrature_generator.clear_quadratures();
1897  surface_quadrature.clear();
1898  }
1899 
1900 
1901 
1902  template <int dim>
1903  void
1905  const BoundingBox<dim> &box,
1906  const unsigned int face_index)
1907  {
1908  clear_quadratures();
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 
1939  quadrature_generator.generate_append(face_restriction, cross_section);
1940 
1941  // We need the dim-dimensional normals of the zero-contour.
1942  // Recompute these.
1943  const ImmersedSurfaceQuadrature<dim - 1, dim - 1>
1944  &surface_quadrature_wrong_normal =
1945  quadrature_generator.get_surface_quadrature();
1946 
1947  std::vector<Tensor<1, dim>> normals;
1948  normals.reserve(surface_quadrature_wrong_normal.size());
1949  for (unsigned int i = 0; i < surface_quadrature_wrong_normal.size(); ++i)
1950  {
1952  surface_quadrature_wrong_normal.point(i),
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  }
1960  surface_quadrature = ImmersedSurfaceQuadrature<dim - 1, dim>(
1961  surface_quadrature_wrong_normal.get_points(),
1962  surface_quadrature_wrong_normal.get_weights(),
1963  normals);
1964  }
1965 
1966 
1967 
1968  template <int dim>
1969  void
1971  {
1972  quadrature_generator.set_1D_quadrature(q_index);
1973  }
1974 
1975 
1976 
1977  template <int dim>
1978  const Quadrature<dim - 1> &
1980  {
1981  return quadrature_generator.get_inside_quadrature();
1982  }
1983 
1984 
1985  template <int dim>
1986  const Quadrature<dim - 1> &
1988  {
1989  return quadrature_generator.get_outside_quadrature();
1990  }
1991 
1992 
1993 
1994  template <int dim>
1995  const ImmersedSurfaceQuadrature<dim - 1, dim> &
1997  {
1998  return surface_quadrature;
1999  }
2000 
2001 
2002 
2004  const hp::QCollection<1> &quadratures1D,
2005  const AdditionalData &additional_data)
2006  {
2007  (void)quadratures1D;
2008  (void)additional_data;
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  {
2051  inside_quadrature = Quadrature<0>(points, weights);
2052  outside_quadrature = Quadrature<0>();
2053  }
2054  else if (level_set_value > 0)
2055  {
2056  inside_quadrature = Quadrature<0>();
2057  outside_quadrature = Quadrature<0>(points, weights);
2058  }
2059  else
2060  {
2061  inside_quadrature = Quadrature<0>();
2062  outside_quadrature = Quadrature<0>();
2063  }
2064  }
2065 
2066 
2067 
2068  void
2070  {
2071  (void)q_index;
2072  }
2073 
2074 
2075 
2076  const Quadrature<0> &
2078  {
2079  return inside_quadrature;
2080  }
2081 
2082 
2083  const Quadrature<0> &
2085  {
2086  return outside_quadrature;
2087  }
2088 
2089 
2090 
2093  {
2094  return surface_quadrature;
2095  }
2096 
2097 
2098 
2099  template <int dim>
2100  template <typename VectorType>
2102  const hp::QCollection<1> &quadratures1D,
2103  const DoFHandler<dim> &dof_handler,
2104  const VectorType &level_set,
2105  const AdditionalData &additional_data)
2106  : QuadratureGenerator<dim>(quadratures1D, additional_data)
2107  , reference_space_level_set(
2108  std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
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 
2126  this->clear_quadratures();
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 =
2157  internal::DiscreteQuadratureGeneratorImplementation::
2158  setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
2159 
2160  const auto subcell_box =
2161  internal::DiscreteQuadratureGeneratorImplementation::
2162  create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
2163 
2164  reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
2165 
2166  QuadratureGenerator<dim>::generate_append(*reference_space_level_set,
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(),
2179  internal::DiscreteQuadratureGeneratorImplementation::
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
2187  QuadratureGenerator<dim>::generate(*reference_space_level_set, unit_box);
2188  }
2189 
2190 
2191 
2192  template <int dim>
2193  template <typename VectorType>
2195  const hp::QCollection<1> &quadratures1D,
2196  const DoFHandler<dim> &dof_handler,
2197  const VectorType &level_set,
2198  const AdditionalData &additional_data)
2199  : FaceQuadratureGenerator<dim>(quadratures1D, additional_data)
2200  , reference_space_level_set(
2201  std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
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 
2220  this->clear_quadratures();
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 =
2258  internal::DiscreteQuadratureGeneratorImplementation::
2259  setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
2260 
2261  const auto subcell_box =
2262  internal::DiscreteQuadratureGeneratorImplementation::
2263  create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
2264 
2265  reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
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(),
2281  internal::DiscreteQuadratureGeneratorImplementation::
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
2289  FaceQuadratureGenerator<dim>::generate(*reference_space_level_set,
2290  unit_box,
2291  face_index);
2292  }
2293 } // namespace NonMatching
2294 #include "quadrature_generator.inst"
BoundingBox< 1, Number > bounds(const unsigned int direction) const
Point< spacedim, Number > center() const
Number lower_bound(const unsigned int direction) const
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
DiscreteFaceQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const DoFHandler< dim > &dof_handler, const VectorType &level_set, const AdditionalData &additional_data=AdditionalData())
void generate(const typename Triangulation< dim >::active_cell_iterator &cell, const unsigned int face_index)
void generate_fe_q_iso_q1(const BoundingBox< dim > &unit_box, unsigned int face_index)
void generate(const typename Triangulation< dim >::active_cell_iterator &cell)
void generate_fe_q_iso_q1(const BoundingBox< dim > &unit_box)
DiscreteQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const DoFHandler< dim > &dof_handler, const VectorType &level_set, const AdditionalData &additional_data=AdditionalData())
const Quadrature< dim - 1 > & get_inside_quadrature() const
void generate(const Function< dim > &level_set, const BoundingBox< dim > &box, const unsigned int face_index)
ImmersedSurfaceQuadrature< dim - 1, dim > surface_quadrature
void set_1D_quadrature(const unsigned int q_index)
void generate_append(const Function< dim > &level_set, const BoundingBox< dim > &box, const unsigned int face_index)
FaceQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const AdditionalData &additional_data=AdditionalData())
const ImmersedSurfaceQuadrature< dim - 1, dim > & get_surface_quadrature() const
const Quadrature< dim - 1 > & get_outside_quadrature() const
void push_back(const Point< dim > &point, const double weight, const Tensor< 1, spacedim > &normal)
void generate_append(const Function< dim > &level_set, const BoundingBox< dim > &box)
void set_1D_quadrature(const unsigned int q_index)
QuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const AdditionalData &additional_data=AdditionalData())
const Quadrature< dim > & get_inside_quadrature() const
const Quadrature< dim > & get_outside_quadrature() const
const ImmersedSurfaceQuadrature< dim > & get_surface_quadrature() const
void generate(const Function< dim > &level_set, const BoundingBox< dim > &box)
void set_active_cell(const typename Triangulation< dim >::active_cell_iterator &cell) override
Tensor< 1, dim > gradient(const Point< dim > &point, const unsigned int component=0) const override
double value(const Point< dim > &point, const unsigned int component=0) const override
void set_subcell(const std::vector< unsigned int > &mask, const BoundingBox< dim > &subcell_box) override
RefSpaceFEFieldFunction(const DoFHandler< dim > &dof_handler, const VectorType &dof_values)
SymmetricTensor< 2, dim > hessian(const Point< dim > &point, const unsigned int component=0) const override
QGeneratorBase(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
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)
QGenerator(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
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)
ExtendableQuadrature< dim > & quadrature_by_definiteness(const Definiteness definiteness)
void find_roots(const std::vector< std::reference_wrapper< const Function< 1 >>> &functions, const BoundingBox< 1 > &interval, std::vector< double > &roots)
void 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)
UpThroughDimensionCreator(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
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)
void add_tensor_product(const Quadrature< dim - 1 > &lower, const Quadrature< 1 > &quadrature1D, const double start, const double end, const unsigned int component_in_dim, ExtendableQuadrature< dim > &quadrature)
BoundingBox< dim > right_half(const BoundingBox< dim > &box, const unsigned int direction)
void map_quadrature_to_box(const Quadrature< dim > &unit_quadrature, const BoundingBox< dim > &box, ExtendableQuadrature< dim > &quadrature)
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)