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