From: Simon Sticko Date: Mon, 3 May 2021 09:46:14 +0000 (+0200) Subject: Add class generating high-order immersed quadrature rules X-Git-Tag: v9.4.0-rc1~1234^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=169092199e25b3e285bef6bfb425092e8d2570fa;p=dealii.git Add class generating high-order immersed quadrature rules Add a class NonMatching::QuadratureGenerator that generates high-order immersed quadrature rules over a BoundingBox, B, when the domain is described by a level set function, \psi. Thus creating quadrature rules for the following 3 regions of the box: {x \in B : \psi(x) < 0}, {x \in B : \psi(x) > 0}, {x \in B : \psi(x) = 0}. --- diff --git a/doc/doxygen/images/immersed_quadratures.svg b/doc/doxygen/images/immersed_quadratures.svg new file mode 100644 index 0000000000..063bb8b212 --- /dev/null +++ b/doc/doxygen/images/immersed_quadratures.svg @@ -0,0 +1,462 @@ + + + + + + + + + 2021-05-19T11:42:51.997656 + image/svg+xml + + + Matplotlib v3.3.2, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/include/deal.II/non_matching/quadrature_generator.h b/include/deal.II/non_matching/quadrature_generator.h new file mode 100644 index 0000000000..6d105f2ec7 --- /dev/null +++ b/include/deal.II/non_matching/quadrature_generator.h @@ -0,0 +1,1046 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_non_matching_quadrature_generator_h +#define dealii_non_matching_quadrature_generator_h + +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN +namespace NonMatching +{ + namespace internal + { + namespace QuadratureGeneratorImplementation + { + template + class QGenerator; + } + } // namespace internal + + + /** + * Struct storing settings for the QuadratureGenerator class. + */ + struct AdditionalQGeneratorData + { + /** + * Constructor. + */ + 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); + + /** + * The number of times we are allowed to split the incoming box + * and recurse on each child. + */ + unsigned int max_box_splits; + + /** + * For a level set function, $\psi$, the implicit function theorem states + * that it is possible to write one of the coordinates $x_i$ as a function + * of the others if + * + * $|\frac{\partial \psi}{\partial x_i}| > 0$. + * + * In practice, it might happend the bound we have for the expression in + * the left-hand side is only floating-point close to zero. + * + * This constant is a safety margin, $C$, that states that the implicit + * function theorem can be used when + * + * $|\frac{\partial \psi}{\partial x_i}| > C$ + * + * Thus this constant must be non-negative. + */ + double lower_bound_implicit_function; + + /** + * If two roots are closer to each other than this distance they are + * merged to one. + */ + double min_distance_between_roots; + + /** + * A constant, $C$, controlling when a level set function, $\psi$, is + * considered positive or negative definite: + * + * $\psi(x) > C \Rightarrow \text{Positive definite}$, + * $\psi(x) < -C \Rightarrow \text{Negative definite}$. + */ + double limit_to_be_definite; + + /** + * Tolerance for convergence of the underlying root finder. + */ + double root_finder_tolerance; + + /** + * The number of times the underlying rootfinder is allowed to split + * an interval, while trying to find multiple roots. + */ + unsigned int max_root_finder_splits; + + /** + * This determines how a box is split when this is necessary. If true, the + * box is split in two, if set to false the box is split into its $2^{dim}$ + * children. + */ + bool split_in_half; + }; + + + + /** + * This class creates immersed quadrature rules over a BoundingBox, + * $B \subset \mathbb{R}^{dim}$, when the domain is described by a level set + * function, $\psi$. + * + * This class creates quadrature rules for the intersections between the box + * and the three different regions defined by the level set function. That is, + * it creates quadrature rules to integrate over the following regions + * @f[ + * N = \{x \in B : \psi(x) < 0 \}, \\ + * P = \{x \in B : \psi(x) > 0 \}, \\ + * S = \{x \in B : \psi(x) = 0 \}. + * @f] + * @image html immersed_quadratures.svg + * + * When working with level set functions, the most common is to describe a + * domain, $\Omega$, as + * @f[ + * \Omega = \{ x \in \mathbb{R}^{dim} : \psi(x) < 0 \}. + * @f] + * Given this, we shall use the name convention that $N$ is the "inside" + * region (i.e. inside $\Omega$), $P$ is the "outside" region and $S$ is + * the "surface" region. The "inside" and "outside" quadratures will also be + * referred to as the "bulk"-quadratures. + * + * The underlying algorithm use a 1-dimensional quadrature rule as base for + * creating the immersed quadrature rules. Gauss-Legendre quadrature + * (QGauss) is recommended. The constructor takes an hp::QCollection<1>. + * One can select which 1D-quadrature in the collection should be used + * through the set_1D_quadrature() function. The number of quadrature points + * in the constructed quadratures will vary depending on the level set + * function. More quadrature points will be created if the intersection is + * "bad", for example, if the zero-contour has a high curvature compared to + * the size of the box. However, if the number of points in the 1D quadrature + * is $n$ the number of points will be proportional to $n^{dim}$ in the bulk + * quadratures and to $n^{dim-1}$ in the surface quadrature. For example, + * in the 2D-example in the above figure, there are 2 points in the + * 1D-quadrature. If the 1D-quadrature is a Gauss-Legendre quadrature and the + * grid has size $h$, the immersed quadratures typically give global errors + * proportional to $h^{2n}$, both for the bulk and surface integrals. If the + * 1D-quadrature has positive weights, the weights of the immersed quadratures + * will also be positive. + * + * A detailed description of the underlying algorithm can be found in + * "High-Order %Quadrature Methods for Implicitly Defined Surfaces and + * Volumes in Hyperrectangles, R. I. Saye, SIAM J. Sci. Comput., 37(2), + * doi:10.1137/140966290". This implementation has some modifications + * compared to the algorithm description in the paper. In particular, it + * builds the three different types of quadratures (inside, outside and + * surface) simultaneously. Further, the so-called "pruning" step is not yet + * implemented. + */ + template + class QuadratureGenerator + { + public: + using AdditionalData = AdditionalQGeneratorData; + + /** + * Constructor. Each Quadrature<1> in @p quadratures1D can be chosen as base + * for generating the immersed quadrature rules. + * + * @note It is important that each 1D-quadrature rule in the + * hp::QCollection does not contain the points 0 and 1. + */ + QuadratureGenerator( + const hp::QCollection<1> &quadratures1D, + const AdditionalData & additional_data = AdditionalData()); + + /** + * Construct immersed quadratures rules for the incoming level set + * function over the BoundingBox. + * + * To get the constructed quadratures, use the functions + * get_inside_quadrature(), + * get_outside_quadrature(), + * get_surface_quadrature(). + * + * @note Both value, gradient and hessian need to be implemented on the + * incoming function. + */ + void + generate(const Function &level_set, const BoundingBox &box); + + /** + * Return the quadrature rule for the region + * $\{x \in B : \psi(x) < 0 \}$ + * created in the previous call to generate(). + * Here, $B$ is BoundingBox passed to generate(). + */ + const Quadrature & + get_inside_quadrature() const; + + /** + * Return the quadrature rule for the region + * $\{x \in B : \psi(x) > 0 \}$ + * created in the previous call to generate(). + * Here, $B$ is BoundingBox passed to generate(). + */ + const Quadrature & + get_outside_quadrature() const; + + /** + * Return the quadrature rule for the region + * $\{x \in B : \psi(x) = 0 \}$ + * created in the previous call to generate(). + * Here, $B$ is BoundingBox passed to generate(). + * + * @note The normal at the quadrature points will be parallel to $\nabla \psi$. + */ + const ImmersedSurfaceQuadrature & + get_surface_quadrature() const; + + /** + * Set which 1D-quadrature in the collection passed to the constructor + * should be used to create the immersed quadratures. + */ + void + set_1D_quadrature(const unsigned int q_index); + + private: + /** + * QuadratureGenerator is mainly used to start up the recursive + * algorithm. This is the object that actually generates the quadratures. + */ + internal::QuadratureGeneratorImplementation::QGenerator + q_generator; + }; + + namespace internal + { + namespace QuadratureGeneratorImplementation + { + /** + * A class that attempts to find multiple distinct roots of a function, + * $f(x)$, over an interval, $[l, r]$. This is done as follows. If there + * is a sign change in function value between the interval end points, + * we solve for the root. If there is no sign change, we attempt to + * bound the function value away from zero on $[a, b]$, to conclude that + * no roots exist. If we can't exclude that there are roots, we split + * the interval in two: $[l, (r + l) / 2]$, $[(r + l) / 2, r]$, and use + * the same algorithm recursively on each interval. This means that we + * can typically find 2 distinct roots, but not 3. + * + * The bounds on the functions values are estimated using the function + * taylor_estimate_function_bounds, which approximates the function as a + * second order Taylor-polynomial around the interval midpoint. + * When we have a sign change on an interval, this class uses + * boost::math::tools::toms748_solve for finding roots . + */ + class RootFinder + { + public: + /** + * Struct storing settings for the RootFinder class. + */ + struct AdditionalData + { + /** + * Constructor. + */ + AdditionalData(const double tolerance = 1e-12, + const unsigned int max_recursion_depth = 2, + const unsigned int max_iterations = 500); + + /** + * The tolerance in the stopping criteria for the underlying root + * finding algorithm boost::math::tools::toms748_solve. + */ + double tolerance; + + /** + * The number of times we are allowed to split the interval where we + * seek roots. + */ + unsigned int max_recursion_depth; + + /** + * The maximum number of iterations in + * boost::math::tools::toms748_solve. + */ + unsigned int max_iterations; + }; + + + /** + * Constructor. + */ + RootFinder(const AdditionalData &data = AdditionalData()); + + /** + * For each of the incoming @p functions, attempt to find the roots over + * the interval defined by @p interval and add these to @p roots. + * The returned roots will be sorted in ascending order: + * $x_0 < x_1 <...$ and duplicate roots (with respect to the tolerance + * in AdditionalData) will be removed. + */ + void + find_roots(const std::vector>> + & functions, + const BoundingBox<1> &interval, + std::vector & roots); + + private: + /** + * Attempt to find the roots of the @p function over the interval defined by + * @p interval and add these to @p roots. @p recursion_depth holds the number + * of times this function has been called recursively. + */ + void + find_roots(const Function<1> & function, + const BoundingBox<1> &interval, + const unsigned int recursion_depth, + std::vector & roots); + + const AdditionalData additional_data; + }; + + + /** + * This is a special Quadrature class with a push_back() method for + * conveniently adding a point with an associated weight. + * + * Since we build the quadrature rules in step-wise fashion, + * it's easier to use this class than to pass around two vectors: + * std::vector>, + * std::vector. + * Further, two std::vectors could accidentally end up with different + * sizes. Using push_back we make sure that the number of points and + * weights are the same. + */ + template + class ExtendableQuadrature : public Quadrature + { + public: + /** + * Constructor, creates an empty quadrature rule with no + * points. + */ + ExtendableQuadrature() = default; + + /** + * Constructor, copies the incoming Quadrature. + */ + ExtendableQuadrature(const Quadrature &quadrature); + + /** + * Add a point with an associated weight to the quadrature. + */ + void + push_back(const Point &point, const double weight); + }; + + + /** + * Type that describes the definiteness of a function over a region. + */ + enum class Definiteness + { + negative, + positive, + indefinite + }; + + + /** + * Class that stores quadrature rules to integrate over 4 different + * regions of a single BoundingBox, $B$. Given multiple level set + * functions, + * + * $\psi_i : \mathbb{R}^{dim} \rightarrow \mathbb{R}$, $i = 0, 1, ...$, + * + * the box, $B \subset \mathbb{R}^{dim}$, is partitioned into a + * "negative", "positive", and "indefinite" region, $B = N \cup P \cup I$, + * according to the signs of $\psi_i$ over each region: + * + * @f[ + * N = \{x \in B : \psi_i(x) < 0, \forall i \}, \\ + * P = \{x \in B : \psi_i(x) > 0, \forall i \}, \\ + * I = B \setminus (\overline{N} \cup \overline{P}). + * @f] + * + * Thus, all $\psi_i$ are positive over $P$ and negative over $N$. Over + * $I$ the level set functions differ in sign. This class holds quadrature + * rules for each of these regions. In addition, when there is a single + * level set function, $\psi$, it holds a surface quadrature for the zero + * contour of $\psi$: + * + * $S = \{x \in B : \psi(x) = 0 \}$. + * + * Note that when there is a single level set function, $I$ is empty + * and $N$ and $P$ are the regions that one typically integrates over in + * an immersed finite element method. + */ + template + class QPartitioning + { + public: + /** + * Return a reference to the "bulk" quadrature with the same name as the + * member in Definiteness. + */ + ExtendableQuadrature & + quadrature_by_definiteness(const Definiteness definiteness); + + /** + * Quadrature for the region $\{x \in B : \psi_i(x) < 0 \forall i \}$ of + * the box, $B$. + */ + ExtendableQuadrature negative; + + /** + * Quadrature for the region $\{x \in B : \psi_i(x) > 0 \forall i \}$ of + * the box, $B$. + */ + ExtendableQuadrature positive; + + /** + * Quadrature for a region where the level set functions have different + * sign. + */ + ExtendableQuadrature indefinite; + + /** + * Quadrature for the region $\{x \in B : \psi(x) = 0 \}$ of the + * box, $B$. + */ + ImmersedSurfaceQuadrature surface; + }; + + + /** + * This class is responsible for creating quadrature points for + * the $dim$-dimensional quadrature partitioning from an + * $(dim - 1)$-dimensional "indefinite" quadrature (see + * QPartitioning documentation). + * + * To be precise, let $[L, R]$ be the extents of the box in the height + * function direction and let $I \subset \mathbb{R}^{dim-1}$ be the lower + * dimensional indefinite region. This class will create quadrature points + * over $I \times [L, R] \subset \mathbb{R}^{dim}$ and in the case + * $dim=spacedim$, points for the surface quadrature. + * + * For each lower dimensional quadrature point, $(x_I, w_I)$ in the + * indefinite quadrature, we create several 1D-level set functions by + * restricting $\psi_j$ to $x_I$. We then partition the interval $[L, R]$ + * into $[y_0, y_1, ..., y_n]$, where $y_0 = L$, $y_n = R$, and the + * remaining $y_i$ are the roots of the 1D-level set functions in + * $[L, R]$. Since the level set functions change sign between the + * roots, each interval belong to different regions in the quadrature + * partitioning. + * + * In each interval, $[y_i, y_{i+1}]$, we distribute points + * according to the 1D-base quadrature, $(x_q, w_q)$ and take the + * cartesian product with $(x_I, w_I)$ to create the $dim$-dimensional + * quadrature points, $(X_q, W_q)$: + * $X_q = x_I \times (y_i + (y_{i+1} - y_i) x_q)$, + * $W_q = w_I (y_{i+1} - y_i) w_q$. + * + * When $dim=spacedim$, we have a single level set function, $\psi$. Since + * we have fulfilled the implicit function theorem, there is a single root + * $y_1 \in [L, R]$. The point, $x_s = x_I \times y_1$, will be added as a + * point in the surface quadrature. One can show that the correct weight + * of this point is + * + * $w_s = \frac{\|\nabla \psi(x_s)\|}{|\partial_i \psi(x_s)|} w_I$, + * + * where $i$ is the height function direction. + */ + template + class UpThroughDimensionCreator + { + public: + /** + * Constructor. Takes the same parameters as QuadratureGenerator. + */ + UpThroughDimensionCreator( + const hp::QCollection<1> & q_collection1D, + const AdditionalQGeneratorData &additional_data); + + /** + * Create $dim$-dimensional immersed quadratures from the incoming + * $(dim-1)$-dimensional quadratures and add these to + * @p q_partitioning. + */ + void + generate(const std::vector>> + & level_sets, + const BoundingBox & box, + const Quadrature &low_dim_quadrature, + const unsigned int height_function_direction, + QPartitioning & q_partitioning); + + /** + * Set which 1D-quadrature in the collection passed to the constructor + * should be used to create the immersed quadratures. + */ + void + set_1D_quadrature(const unsigned int q_index); + + private: + /** + * Create a surface quadrature point from the lower-dimensional point + * and add it to surface_quadrature. + * + * This function is only called when $dim=spacedim$ and there is a + * single level set function. At this point there should only be a + * single root in the interval $[L, R]$ + */ + void + create_surface_point( + const Point &point, + const double weight, + const std::vector>> + & level_sets, + const BoundingBox & box, + const unsigned int height_function_direction, + ImmersedSurfaceQuadrature &surface_quadrature); + + /** + * One dimensional quadrature rules used to create the immersed + * quadratures. + */ + const SmartPointer> q_collection1D; + + /** + * Stores options/settings for the algorithm. + */ + const AdditionalQGeneratorData additional_data; + + /** + * Which quadrature rule in the above collection that is used to + * create the immersed quadrature rules. + */ + unsigned int q_index; + + /** + * 1D-functions, that are restrictions of each dim-dimensional level set + * function passed to generate() to some $(dim-1)$-dimensional point. + */ + std::vector> point_restrictions; + + /** + * Class used to find the roots of the above 1D-restictions. + */ + RootFinder root_finder; + + /** + * The roots of the functions in point_restrictions. + * This will be the values of the height functions, $\{H_i(x_I)\}$ at + * some lower dimensional quadrature point, + * $x_I \in \mathbb{R}^{dim-1}}$. + */ + std::vector roots; + }; + + + /** + * Data representing the best choice of height-function direction, + * which is returned by the function find_best_height_direction. + * + * This data consists of a coordinate direction + * + * $i \in \{0, ..., dim - 1 \}$, + * + * and lower bound on the absolute value of the derivative of some + * associated function, f, taken in the above coordinate direction. That + * is, a bound $C$ such that + * + * $|\frac{\partial f}{\partial x_i}| > C$, + * + * holding over some subset of $\mathbb{R}^{dim}$. + */ + struct HeightDirectionData + { + /** + * Constructor. Initializes the direction to invalid_unsigned_int and + * the bound to 0. + */ + HeightDirectionData(); + + + /** + * The height-function direction, described above. + */ + unsigned int direction; + + /** + * The lower bound on $|\frac{\partial f}{\partial x_i}|$, described + * above. + */ + double min_abs_dfdx; + }; + + + /** + * Base class for the class QGenerator and the + * one-dimensional specialization QGenerator<1, spacedim>. + */ + template + class QGeneratorBase + { + public: + QGeneratorBase(const hp::QCollection<1> & q_collection1D, + const AdditionalQGeneratorData &additional_data); + + /** + * Clear the quadratures created by the previous call to generate(). + */ + void + clear_quadratures(); + + /** + * Return the created quadratures. + */ + const QPartitioning & + get_quadratures() const; + + protected: + /** + * Stores options/settings for the algorithm. + */ + const AdditionalQGeneratorData additional_data; + + /** + * Which 1D-quadrature in the collection we should use to generate + * the immersed quadrature. + */ + unsigned int q_index; + + /** + * Index of the quadrature in q_collection1D that should use to + * generate the immersed quadrature rules. + */ + const SmartPointer> q_collection1D; + + /** + * Quadratures that the derived classes create. + */ + QPartitioning q_partitioning; + }; + + + /** + * This class implements the Saye-algorithm cited in the documentation of + * the QuadratureGenerator class. + * + * The generate function takes a number of $dim$-dimensional level set + * functions, $\psi_i$, and a BoundingBox, and builds a partitioning + * of quadratures, as defined in documentation of the QPartitioning class. + * That is, this class builds an object of type QPartitioning. + * + * If all $\psi_i$ passed to generate can be determined to be positive or + * negative definite, the QPartitioning will consist of a single + * quadrature forming a tensor product. + * + * If this is not the case, the algorithm uses recursion over the spatial + * dimension. The spacedim template parameter denotes the dimension we + * started with and dim denotes on what level we are in the recursion. + * That is, we first construct a QPartitioning and then + * build the higher dimensional quadratures from these. What we in the end + * actually want is a spacedim-dimensional partitioning of quadratures, + * for a single level set function, $\psi$. + * + * The algorithm is based on the implicit function theorem. Starting with + * a single level set function, $\psi$, we try to find a direction $i$, + * such that + * + * $|\frac{\partial \psi}{\partial x_i}| > 0$. + * + * throughout the whole box. This means that the zero-contour of the + * level set function can be parameterized by an implicit function + * + * $H = H(x_0, ..., x_{i-1}, x_{i+1}, ..., x_{dim-1})$, + * + * so that + * + * $\psi(..., x_{i-1}, H(..., x_{i-1}, x_{i+1}, ...), x_{i+1}, ...) = 0$, + * + * over a subset, $I \subset C \subset \mathbb{R}^{dim-1}$, of the cross + * section, $C$, of the box (see BoundingBox::cross_section). Here, $I$ is + * the "indefinite"-region defined in the QPartitioning class. To follow + * convention in the original paper, we will -refer to $H$ as the + * "height-function" and to $i$ as the "height-function direction". + * + * If a height function direction can be found, we go down in dimension by + * creating two new level set functions, $\{\psi_0, \psi_1\}$, which are + * the restriction of $\psi$ to the top and bottom faces of the box (in + * the height function direction). We then delegate to + * QGenerator to create a QPartitioning over + * the cross section. + * + * When we reach the base case, $dim = 1$, the creation of + * QPartitioning<1> is simple. See the documentation in specialized + * class: QGenerator<1, spacedim>. + * + * As we go up through the dimensions and create the higher dimensional + * quadratures, we need to know the function value of the height + * functions at the lower dimensional quadrature points. Since the + * functions are implicit, we need to do root-finding on the level set + * functions to find the function values. For this we use the class + * UpThroughDimensionCreator, see documentation there. + * + * When we have $n$ level set functions (i.e. after having gone + * down in dimension), we try to find a height function direction, + * which works for all those $\psi_i$ which are intersected by the zero + * contour (i.e. those not positive or negative definite). + * If such a direction exist, we will have a maximum of $n$ associated + * implicit height functions, $H_j$. Each $H_j$ parametrize the + * $x_i$-coordinate of the zero-contour over a region, $I_j$. The + * indefinite region in the lower dimensional partitioning is the union of + * these $I = \cup_j I_j$. + * + * As we try to find a height function direction, we estimate bounds on + * the gradient components by approximating each component as a 1st-order + * Taylor-polynomial. If a direction can not be found, the box is split + * and we recurse on each smaller box. This makes an implicit function + * more likely to exist since we seek it over a smaller portion of the + * zero contour. It also makes the estimated bounds tighter since we + * extrapolate the Taylor-polynomial a shorter distance. + * + * Since we can not split a box forever, there is an maximum number of + * allowed splits on the additional data struct passed to the constructor. + * If this is reached, the algorithm uses the midpoint method as a last + * resort. + */ + template + class QGenerator : public QGeneratorBase + { + public: + /** + * Constructor. Takes the same parameters QuadratureGenerator. + */ + QGenerator(const hp::QCollection<1> & q_collection1D, + const AdditionalQGeneratorData &additional_data); + + /** + * Create immersed quadrature rules over the incoming @p box and add + * these to the internal QPartitioning object in the base class. + * These quadratures can then be obtained using the + * get_quadratures-function. + * + * This function calls itself if the incoming box need to be split. + * @p n_box_splits counts the number of times this function has called + * itself. + */ + void + generate(const std::vector>> + & level_sets, + const BoundingBox &box, + const unsigned int n_box_splits); + + /** + * Set which 1D-quadrature in the collection passed to the constructor + * should be used to create the immersed quadratures. + */ + void + set_1D_quadrature(const unsigned int q_index); + + private: + /** + * Restricts the incoming level set functions to the top and bottom of + * the incoming box (w.r.t @p height_function_direction). Then call the + * lower dimensional QGenerator with the cross section of the box + * to generate the lower dimensional immersed quadrature rules. + */ + void + create_low_dim_quadratures( + const unsigned int height_function_direction, + const std::vector>> + & level_sets, + const BoundingBox &box, + const unsigned int n_box_splits); + + /** + * Gets the $(dim - 1)$-dimensional quadratures from the lower + * dimensional algorithm and creates the $dim$-dimensional quadrature + * rules over the box from the lower dimensional ones. + */ + void + create_high_dim_quadratures( + const unsigned int height_function_direction, + const std::vector>> + & level_sets, + const BoundingBox &box); + + /** + * Split the incoming box and call generate() recursively with each box. + * The box is split in 2 or 4 parts depending on the value of + * AdditionalQGeneratorData::split_in_half. + */ + void + split_box_and_recurse( + const std::vector>> + & level_sets, + const BoundingBox & box, + const std_cxx17::optional &direction_data, + const unsigned int n_box_splits); + + /** + * Uses the midpoint-method to create a quadrature over the box. + * That is, add a single quadrature point at the center of the box + * with weight corresponding to the volume of the box. + * + * The point is added to the region defined in QPartitioning + * according to the signs of the level set functions at the center of + * the box. + */ + void + use_midpoint_method( + const std::vector>> + & level_sets, + const BoundingBox &box); + + /** + * The same algorithm as this, but creating immersed quadratures + * in one dimension lower. + */ + QGenerator low_dim_algorithm; + + /** + * Object responsible for creating the $dim$-dimensional quadratures + * from + */ + UpThroughDimensionCreator up_through_dimension_creator; + + /** + * Stores tensor products of each of the Quadrature<1>'s in + * q_collection1D. + */ + hp::QCollection tensor_products; + }; + + + /** + * The 1D-base case of the recursive algorithm QGenerator. + * + * Let $L$ and $R$ be the left and right bounds of the one-dimensional + * BoundingBox. This interval is partitioned into $[x_0, x_1, ..., x_n]$ + * where $x_0 = L$, $x_n = R$, and the remaining $x_i$ are the roots + * of the level set functions in the interval $[L, R]$. In each interval, + * $[x_i, x_{i+1}]$, quadrature points are distributed according to a + * 1D-quadrature rule. These points are added to one of the regions of + * QPartitioning determined from the signs of the level set + * functions on the interval (see documentation of QPartitioning). + * + * If spacedim = 1 the points $[x_1, x_{n-1}]$ are also added as surface + * quadrature points to QPartitioning::surface. + */ + template + class QGenerator<1, spacedim> : public QGeneratorBase<1, spacedim> + { + public: + /** + * Constructor. Takes the same parameters QuadratureGenerator. + */ + QGenerator(const hp::QCollection<1> & quadratures1D, + const AdditionalQGeneratorData &additional_data); + + /** + * Creates quadrature points over the interval defined by the incoming + * box and adds these quadrature points to the internally stored + * QPartitioning. These quadratures can then be obtained using + * the get_quadratures-function. + */ + void + generate(const std::vector>> + & level_sets, + const BoundingBox<1> &box, + const unsigned int n_box_splits); + + /** + * Set which 1D-quadrature in the collection passed to the constructor + * should be used to create the immersed quadratures. + */ + void + set_1D_quadrature(const unsigned int q_index); + + private: + /** + * Adds the point defined by coordinate to the surface quadrature of + * ImmersedQuadrature with unit weight. + */ + void + create_surface_points( + const std::vector>> + &level_sets); + + /** + * Class used to find the roots of the functions passed to generate(). + */ + RootFinder root_finder; + + /** + * Roots of the functions passed to generate(). + */ + std::vector roots; + + /** + * This would be the height-function direction in higher dimensions, + * but in 1D there is only one coordinate direction. + */ + const unsigned int direction = 0; + + /** + * To reuse the distribute_points_between_roots()-function + * we need a zero-dimensional quadrature point with unit weight. + */ + const Point<0> zero_dim_point; + const double unit_weight = 1; + }; + + + /** + * Take the tensor product between (point, weight) and @p quadrature1D + * scaled over [start, end] and add the resulting dim-dimensional + * quadrature points to @p quadrature. + * + * @p component_in_dim specifies which dim-dimensional coordinate + * quadrature1D should be written to. + */ + template + void + tensor_point_with_1D_quadrature(const Point &point, + const double weight, + const Quadrature<1> & quadrature1D, + const double start, + const double end, + const unsigned int component_in_dim, + ExtendableQuadrature &quadrature); + + + /** + * Checks the sign of the incoming Functions at the incoming point and + * returns Definiteness::positive/negative if all the functions are + * positive/negative at the point, otherwise returns + * Definiteness::indefinite. + */ + template + Definiteness + pointwise_definiteness( + const std::vector>> + & functions, + const Point &point); + + + /** + * A struct storing the bounds on the function value and bounds + * on each component of the gradient. + */ + template + struct FunctionBounds + { + public: + /** + * Lower and upper bounds on the functions value. + */ + std::pair value; + + /** + * Lower and upper bounds on each component of the gradient. + */ + std::array, dim> gradient; + }; + + + /** + * Returns the max/min bounds on the value, taken over all the entries + * in the incoming vector of FunctionBounds. That is, given the incoming + * function bounds, $[L_j, U_j]$, this function returns + * $[L, U]$, + * where $L = \min_{j} L_j$ and $U = \max_{j} U_j$. + */ + template + std::pair + find_extreme_values( + const std::vector> &all_function_bounds); + + + /** + * Finds the best choice of height function direction, given the + * FunctionBounds for a number of functions $\{\psi_j\}_{j=0}^{n-1}$. + * Here, "best" is meant in the sense of the implicit function theorem. + * + * Let $J_I$ be the index set of the indefinite functions: + * + * $J_I = \{0,..., n - 1\} \setminus \{ j : |\psi_j| > 0 \}$. + * + * This function converts the incoming bounds to a lower bound, $L_{ij}$, + * on the absolute value of each component of the gradient: + * + * $|\partial_k \psi_j| > L_{jk}$. + * + * and then returns a coordindate direction, $i$, and a lower bound $L$, + * such that + * + * @f[ + * i = \arg \max_{k} \min_{j \in J_I} L_{jk}, \\ + * L = \max_{k} \min_{j \in J_I} L_{jk}. + * @f] + * + * This means $i$ is a coordinate direction such that all functions + * intersected by the zero contour (i.e. those belonging to $J_I$) fulfill + * + * $|\partial_i \psi_j| > L$. + * + * Note that the estimated lower bound, $L$, can be zero or negative. This + * means that no suitable height function direction exists. If all of the + * incoming functions are positive or negative definite the returned + * std::optional is non-set. + */ + template + std_cxx17::optional + find_best_height_direction( + const std::vector> &all_function_bounds); + + } // namespace QuadratureGeneratorImplementation + } // namespace internal + +} // namespace NonMatching +DEAL_II_NAMESPACE_CLOSE + +#endif /* dealii_non_matching_quadrature_generator_h */ diff --git a/source/non_matching/CMakeLists.txt b/source/non_matching/CMakeLists.txt index 8ea4a0eb48..2d1b9a8a76 100644 --- a/source/non_matching/CMakeLists.txt +++ b/source/non_matching/CMakeLists.txt @@ -16,11 +16,13 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_src + quadrature_generator.cc coupling.cc immersed_surface_quadrature.cc ) SET(_inst + quadrature_generator.inst.in coupling.inst.in ) diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc new file mode 100644 index 0000000000..ea0245e069 --- /dev/null +++ b/source/non_matching/quadrature_generator.cc @@ -0,0 +1,1358 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN +namespace NonMatching +{ + namespace internal + { + namespace QuadratureGeneratorImplementation + { + template + void + tensor_point_with_1D_quadrature(const Point &point, + const double weight, + const Quadrature<1> & quadrature1D, + const double start, + const double end, + const unsigned int component_in_dim, + ExtendableQuadrature &quadrature) + { + Assert(start < end, + ExcMessage("Interval start must be less than interval end.")); + + const double length = end - start; + for (unsigned int j = 0; j < quadrature1D.size(); ++j) + { + const double x = start + length * quadrature1D.point(j)[0]; + quadrature.push_back(dealii::internal::create_higher_dim_point( + point, component_in_dim, x), + length * weight * quadrature1D.weight(j)); + } + } + + + + /** + * For each (point, weight) in lower create a dim-dimensional quadrature + * using tensor_point_with_1D_quadrature and add the results to @p quadrature. + */ + template + void + add_tensor_product(const Quadrature &lower, + const Quadrature<1> & quadrature1D, + const double start, + const double end, + const unsigned int component_in_dim, + ExtendableQuadrature &quadrature) + { + for (unsigned int j = 0; j < lower.size(); ++j) + { + tensor_point_with_1D_quadrature(lower.point(j), + lower.weight(j), + quadrature1D, + start, + end, + component_in_dim, + quadrature); + } + } + + + + template + Definiteness + pointwise_definiteness( + const std::vector>> + & functions, + const Point &point) + { + Assert(functions.size() > 0, + ExcMessage( + "The incoming vector must contain at least one function.")); + + const int sign_of_first = + boost::math::sign(functions[0].get().value(point)); + + if (sign_of_first == 0) + return Definiteness::indefinite; + + for (unsigned int j = 1; j < functions.size(); ++j) + { + const int sign = boost::math::sign(functions[j].get().value(point)); + + if (sign != sign_of_first) + return Definiteness::indefinite; + } + // If we got here all functions have the same sign. + if (sign_of_first < 0) + return Definiteness::negative; + else + return Definiteness::positive; + } + + + + /** + * Given the incoming lower and upper bounds on the value of a function + * $[L, U]$, return the minimum/maximum of $[L, U]$ and the function + * values at the vertices. That is, this function returns + * + * $[\min(L, L_f), \max(U, U_f)]$, + * + * where $L_f = \min_{v} f(x_v)$, $U_f = \max_{v} f(x_v)|$, + * and $x_v$ is a vertex. + * + * It is assumed that the incoming function is scalar valued. + */ + template + void + take_min_max_at_vertices(const Function & function, + const BoundingBox & box, + std::pair &value_bounds) + { + const ReferenceCell &cube = ReferenceCells::get_hypercube(); + for (unsigned int i = 0; i < cube.n_vertices(); ++i) + { + const double vertex_value = function.value(box.vertex(i)); + + value_bounds.first = std::min(value_bounds.first, vertex_value); + value_bounds.second = std::max(value_bounds.second, vertex_value); + } + } + + + + /** + * Estimate bounds on each of the functions in the incoming vector over + * the incoming box. + * + * Bounds on the functions value and the gradient components are first + * computed using FunctionTools::taylor_estimate_function_bounds. + * In addition, the function value is checked for min/max at the at + * the vertices of the box. The gradient is not checked at the box + * vertices. + */ + template + void + estimate_function_bounds( + const std::vector>> + & functions, + const BoundingBox & box, + std::vector> &all_function_bounds) + { + all_function_bounds.clear(); + all_function_bounds.reserve(functions.size()); + for (const Function &function : functions) + { + FunctionBounds bounds; + FunctionTools::taylor_estimate_function_bounds( + function, box, bounds.value, bounds.gradient); + take_min_max_at_vertices(function, box, bounds.value); + + all_function_bounds.push_back(bounds); + } + } + + + + template + std::pair + find_extreme_values(const std::vector> &bounds) + { + Assert(bounds.size() > 0, ExcMessage("The incoming vector is empty.")); + + std::pair extremes = bounds[0].value; + for (unsigned int i = 1; i < bounds.size(); ++i) + { + extremes.first = std::min(extremes.first, bounds[i].value.first); + extremes.second = std::max(extremes.second, bounds[i].value.second); + } + + return extremes; + } + + + + /** + * Return true if the incoming function bounds correspond to a function + * which is indefinite, i.e., that is not negative or positive definite. + */ + inline bool + is_indefinite(const std::pair &function_bounds) + { + if (function_bounds.first > 0) + return false; + if (function_bounds.second < 0) + return false; + return true; + } + + + + /** + * Return a lower bound, $L_a$, on the absolute value of a function, + * $f(x)$: + * + * $L_a \leq |f(x)|$, + * + * by estimating it from the incoming lower and upper bounds: + * $L \leq f(x) \leq U$. + * + * By rewriting the lower and upper bounds as + * $F - C \leq f(x) \leq F + C$, + * where $L = F - C$, $U = F + C$ (or $F = (U + L)/2$, $C = (U - L)/2$), + * we get $|f(x) - F| \leq C$. + * Using the inverse triangle inequality gives + * $|F| - |f(x)| \leq |f(x) - F| \leq C$. + * Thus, $L_a = |F| - C$. + * + * Note that the returned value can be negative. This is used to indicate + * "how far away" a function is from being definite. + */ + inline double + lower_bound_on_abs(const std::pair &function_bounds) + { + Assert(function_bounds.first <= function_bounds.second, + ExcMessage("Function bounds reversed, max < min.")); + + return 0.5 * (std::abs(function_bounds.second + function_bounds.first) - + (function_bounds.second - function_bounds.first)); + } + + + + HeightDirectionData::HeightDirectionData() + { + direction = numbers::invalid_unsigned_int; + min_abs_dfdx = 0; + } + + + + template + std_cxx17::optional + find_best_height_direction( + const std::vector> &all_function_bounds) + { + // Minimum (taken over the indefinite functions) on the lower bound on + // each component of the gradient. + std_cxx17::optional> min_lower_abs_grad; + + for (const FunctionBounds &bounds : all_function_bounds) + { + if (is_indefinite(bounds.value)) + { + // For the first indefinite function we find, we write the lower + // bounds on each gradient component to min_lower_abs_grad. + if (!min_lower_abs_grad) + { + min_lower_abs_grad.emplace(); + for (int d = 0; d < dim; ++d) + { + (*min_lower_abs_grad)[d] = + lower_bound_on_abs(bounds.gradient[d]); + } + } + else + { + for (int d = 0; d < dim; ++d) + { + (*min_lower_abs_grad)[d] = + std::min((*min_lower_abs_grad)[d], + lower_bound_on_abs(bounds.gradient[d])); + } + } + } + } + + if (min_lower_abs_grad) + { + const auto max_element = + std::max_element(min_lower_abs_grad->begin(), + min_lower_abs_grad->end()); + + HeightDirectionData data; + data.direction = + std::distance(min_lower_abs_grad->begin(), max_element); + data.min_abs_dfdx = *max_element; + + return data; + } + + return std_cxx17::optional(); + } + + + + /** + * Return true if there are exactly two incoming FunctionBounds and + * they corresponds to one function being positive definite and + * one being negative definite. Return false otherwise. + */ + template + inline bool + one_positive_one_negative_definite( + const std::vector> &all_function_bounds) + { + if (all_function_bounds.size() != 2) + return false; + else + { + const FunctionBounds &bounds0 = all_function_bounds.at(0); + const FunctionBounds &bounds1 = all_function_bounds.at(1); + + if (bounds0.value.first > 0 && bounds1.value.second < 0) + return true; + if (bounds1.value.first > 0 && bounds0.value.second < 0) + return true; + return false; + } + } + + + + /** + * Transform the points and weights of the incoming quadrature, + * unit_quadrature, from unit space to the incoming box and add these to + * quadrature. + * + * Note that unit_quadrature should be a quadrature over [0,1]^dim. + */ + template + void + map_quadrature_to_box(const Quadrature & unit_quadrature, + const BoundingBox & box, + ExtendableQuadrature &quadrature) + { + for (unsigned int i = 0; i < unit_quadrature.size(); i++) + { + const Point point = box.unit_to_real(unit_quadrature.point(i)); + const double weight = unit_quadrature.weight(i) * box.volume(); + + quadrature.push_back(point, weight); + } + } + + + + /** + * For each of the incoming dim-dimensional functions, create the + * restriction to the top and bottom of the incoming BoundingBox and add + * these two (dim-1)-dimensional functions to @p restrictions. Here, top and bottom is + * meant with respect to the incoming @p direction. For each function, the + * "bottom-restriction" will be added before the "top-restriction" + * + * @note @p restrictions will be cleared, so after this function + * restrictions.size() == 2 * functions.size(). + */ + template + void + restrict_to_top_and_bottom( + const std::vector>> + & functions, + const BoundingBox & box, + const unsigned int direction, + std::vector> &restrictions) + { + AssertIndexRange(direction, dim); + + restrictions.clear(); + restrictions.reserve(2 * functions.size()); + + const double bottom = box.lower_bound(direction); + const double top = box.upper_bound(direction); + + for (const auto &function : functions) + { + restrictions.push_back(Functions::CoordinateRestriction( + function, direction, bottom)); + restrictions.push_back(Functions::CoordinateRestriction( + function, direction, top)); + } + } + + + + /** + * Restrict each of the incoming @p functions to @p point, + * while keeping the coordinate direction @p open_direction open, + * and add the restriction to @p restrictions. + * + * @note @p restrictions will be cleared, so after this function + * restrictions.size() == functions.size(). + */ + template + void + restrict_to_point( + const std::vector>> + & functions, + const Point & point, + const unsigned int open_direction, + std::vector> &restrictions) + { + AssertIndexRange(open_direction, dim); + + restrictions.clear(); + restrictions.reserve(functions.size()); + for (const auto &function : functions) + { + restrictions.push_back(Functions::PointRestriction( + function, open_direction, point)); + } + } + + + + /** + * Let $\{ y_0, ..., y_{n+1} \}$ be such that $[y_0, y_{n+1}]$ is the + * @p interval and $\{ y_1, ..., y_n \}$ are the @p roots. In each + * subinterval, $[y_i, y_{i+1}]$, distribute point according to the + * 1D-quadrature rule $\{(x_q, w_q)\}_q$ (@p quadrature1D). + * Take the tensor product with the quadrature point $(x, w)$ + * (@p point, @p weight) to create dim-dimensional quadrature points + * @f[ + * X_q = x_I \times (y_i + (y_{i+1} - y_i) x_q), + * W_q = w_I (y_{i+1} - y_i) w_q, + * @f] + * and add these points to @p q_partitioning. + */ + template + void + distribute_points_between_roots( + const Quadrature<1> & quadrature1D, + const BoundingBox<1> & interval, + const std::vector &roots, + const Point & point, + const double weight, + const unsigned int height_function_direction, + const std::vector>> + & level_sets, + const AdditionalQGeneratorData &additional_data, + QPartitioning & q_partitioning) + { + // Make this int to avoid a warning signed/unsigned comparision. + const int n_roots = roots.size(); + + // The number of intervals are roots.size() + 1 + for (int i = -1; i < n_roots; ++i) + { + // Start and end point of the subinterval. + const double start = i < 0 ? interval.lower_bound(0) : roots[i]; + const double end = + i + 1 < n_roots ? roots[i + 1] : interval.upper_bound(0); + + const double length = end - start; + // It might be that the end points of the subinterval are roots. + // If this is the case then the subinterval has length zero. + // Don't distribute points on the subinterval if it is shorter than + // some tolerance. + if (length > additional_data.min_distance_between_roots) + { + // All points on the interval belong to the same region in + // the QPartitioning. Determine the quadrature we should add + // the points to. + const Point<1> center(start + 0.5 * length); + const Definiteness definiteness = + pointwise_definiteness(level_sets, center); + ExtendableQuadrature &target_quadrature = + q_partitioning.quadrature_by_definiteness(definiteness); + + tensor_point_with_1D_quadrature(point, + weight, + quadrature1D, + start, + end, + height_function_direction, + target_quadrature); + } + } + } + + + + RootFinder::AdditionalData::AdditionalData( + const double tolerance, + const unsigned int max_recursion_depth, + const unsigned int max_iterations) + : tolerance(tolerance) + , max_recursion_depth(max_recursion_depth) + , max_iterations(max_iterations) + {} + + + + RootFinder::RootFinder(const AdditionalData &data) + : additional_data(data) + {} + + + + void + RootFinder::find_roots( + const std::vector>> &functions, + const BoundingBox<1> & interval, + std::vector & roots) + { + for (const Function<1> &function : functions) + { + const unsigned int recursion_depth = 0; + find_roots(function, interval, recursion_depth, roots); + } + // Sort and make sure no roots are duplicated + std::sort(roots.begin(), roots.end()); + + const auto roots_are_equal = [this](const double &a, const double &b) { + return std::abs(a - b) < additional_data.tolerance; + }; + roots.erase(unique(roots.begin(), roots.end(), roots_are_equal), + roots.end()); + } + + + + void + RootFinder::find_roots(const Function<1> & function, + const BoundingBox<1> &interval, + const unsigned int recursion_depth, + std::vector & roots) + { + // Compute function values at end points. + const double left_value = function.value(interval.vertex(0)); + const double right_value = function.value(interval.vertex(1)); + + // If we have a sign change we solve for the root. + if (boost::math::sign(left_value) != boost::math::sign(right_value)) + { + const auto lambda = [&function](const double x) { + return function.value(Point<1>(x)); + }; + + const auto stopping_criteria = [this](const double &a, + const double &b) { + return std::abs(a - b) < additional_data.tolerance; + }; + + boost::uintmax_t iterations = additional_data.max_iterations; + + const std::pair root_bracket = + boost::math::tools::toms748_solve(lambda, + interval.lower_bound(0), + interval.upper_bound(0), + left_value, + right_value, + stopping_criteria, + iterations); + + const double root = .5 * (root_bracket.first + root_bracket.second); + roots.push_back(root); + } + else + { + // Compute bounds on the incoming function to check if there are + // roots. If the function is positive or negative on the whole + // interval we do nothing. + std::pair value_bounds; + std::array, 1> gradient_bounds; + FunctionTools::taylor_estimate_function_bounds<1>(function, + interval, + value_bounds, + gradient_bounds); + + // Since we already know the function values at the interval ends we + // might as well check these for min/max too. + const double function_min = + std::min(std::min(left_value, right_value), value_bounds.first); + + // If the functions is positive there are no roots. + if (function_min > 0) + return; + + const double function_max = + std::max(std::max(left_value, right_value), value_bounds.second); + + // If the functions is negative there are no roots. + if (function_max < 0) + return; + + // If we can't say that the function is strictly positive/negative + // we split the interval in half. We can't split forever, so if we + // have reached the max recursion, we stop looking for roots. + if (recursion_depth < additional_data.max_recursion_depth) + { + find_roots(function, + interval.child(0), + recursion_depth + 1, + roots); + find_roots(function, + interval.child(1), + recursion_depth + 1, + roots); + } + } + } + + + + template + ExtendableQuadrature::ExtendableQuadrature( + const Quadrature &quadrature) + : Quadrature(quadrature) + {} + + + + template + void + ExtendableQuadrature::push_back(const Point &point, + const double weight) + { + this->quadrature_points.push_back(point); + this->weights.push_back(weight); + } + + + + template + ExtendableQuadrature & + QPartitioning::quadrature_by_definiteness( + const Definiteness definiteness) + { + switch (definiteness) + { + case Definiteness::negative: + return negative; + case Definiteness::positive: + return positive; + default: + return indefinite; + } + } + + + + /** + * Takes a (dim-1)-dimensional point from the cross-section (orthogonal + * to direction) of the box. Creates the two dim-dimensional points, which + * are the projections from the cross section to the faces of the box and + * returns the point closest to the zero-contour of the incoming level set + * function. + */ + template + Point + face_projection_closest_zero_contour(const Point & point, + const unsigned int direction, + const BoundingBox &box, + const Function & level_set) + { + const Point bottom_point = + dealii::internal::create_higher_dim_point(point, + direction, + box.lower_bound(direction)); + const double bottom_value = level_set.value(bottom_point); + + const Point top_point = + dealii::internal::create_higher_dim_point(point, + direction, + box.upper_bound(direction)); + const double top_value = level_set.value(top_point); + + // The end point closest to the zero-contour is the one with smallest + // absolute value. + return std::abs(bottom_value) < std::abs(top_value) ? bottom_point : + top_point; + } + + + + template + UpThroughDimensionCreator::UpThroughDimensionCreator( + const hp::QCollection<1> & q_collection1D, + const AdditionalQGeneratorData &additional_data) + : q_collection1D(&q_collection1D) + , additional_data(additional_data) + , root_finder( + RootFinder::AdditionalData(additional_data.root_finder_tolerance, + additional_data.max_root_finder_splits)) + { + q_index = 0; + } + + + + template + void + UpThroughDimensionCreator::generate( + const std::vector>> + & level_sets, + const BoundingBox & box, + const Quadrature &low_dim_quadrature, + const unsigned int height_function_direction, + QPartitioning & q_partitioning) + { + const Quadrature<1> &quadrature1D = (*q_collection1D)[q_index]; + + for (unsigned int q = 0; q < low_dim_quadrature.size(); ++q) + { + const Point &point = low_dim_quadrature.point(q); + const double weight = low_dim_quadrature.weight(q); + restrict_to_point(level_sets, + point, + height_function_direction, + point_restrictions); + + // We need a vector of references to do the recursive call. + const std::vector>> + restrictions(point_restrictions.begin(), + point_restrictions.end()); + + const BoundingBox<1> bounds_in_direction = + box.bounds(height_function_direction); + + roots.clear(); + root_finder.find_roots(restrictions, bounds_in_direction, roots); + + distribute_points_between_roots(quadrature1D, + bounds_in_direction, + roots, + point, + weight, + height_function_direction, + restrictions, + additional_data, + q_partitioning); + + if (dim == spacedim) + create_surface_point(point, + weight, + level_sets, + box, + height_function_direction, + q_partitioning.surface); + } + + point_restrictions.clear(); + } + + + + template + void + UpThroughDimensionCreator::create_surface_point( + const Point &point, + const double weight, + const std::vector>> + & level_sets, + const BoundingBox & box, + const unsigned int height_function_direction, + ImmersedSurfaceQuadrature &surface_quadrature) + { + AssertIndexRange(roots.size(), 2); + Assert(level_sets.size() == 1, ExcInternalError()); + + + const Function &level_set = level_sets.at(0); + + Point surface_point; + if (roots.size() == 1) + { + surface_point = dealii::internal::create_higher_dim_point( + point, height_function_direction, roots[0]); + } + else + { + // If we got here, we have missed roots in the lower dimensional + // algorithm. This is a rare event but can happen if the + // zero-contour has a high curvature. The problem is that the + // incoming point has been incorrectly added to the indefinite + // quadrature in QPartitioning. Since we missed a root on + // this box, we will likely miss it on the neighboring box too. If + // this happens, the point will NOT be in the indefinite quadrature + // on the neighbor. The best thing we can do is to compute the + // surface point by projecting the lower dimensional point to the + // face closest to the zero-contour. We should add a surface point + // because the neighbor will not. + surface_point = face_projection_closest_zero_contour( + point, height_function_direction, box, level_set); + } + + const Tensor<1, dim> gradient = level_set.gradient(surface_point); + Tensor<1, dim> normal = gradient; + normal *= 1. / normal.norm(); + + // Note that gradient[height_function_direction] is non-zero + // because of the implicit function theorem. + const double surface_weight = + weight * gradient.norm() / + std::abs(gradient[height_function_direction]); + surface_quadrature.push_back(surface_point, surface_weight, normal); + } + + + + template + void + UpThroughDimensionCreator::set_1D_quadrature( + unsigned int q_index) + { + AssertIndexRange(q_index, q_collection1D->size()); + this->q_index = q_index; + } + + + + template + QGeneratorBase::QGeneratorBase( + const hp::QCollection<1> & q_collection1D, + const AdditionalQGeneratorData &additional_data) + : additional_data(additional_data) + , q_collection1D(&q_collection1D) + { + q_index = 0; + } + + + + template + QGenerator::QGenerator( + const hp::QCollection<1> & q_collection1D, + const AdditionalQGeneratorData &additional_data) + : QGeneratorBase(q_collection1D, additional_data) + , low_dim_algorithm(q_collection1D, additional_data) + , up_through_dimension_creator(q_collection1D, additional_data) + { + for (unsigned int i = 0; i < q_collection1D.size(); i++) + tensor_products.push_back(Quadrature(q_collection1D[i])); + } + + + + template + void + QGeneratorBase::clear_quadratures() + { + q_partitioning = QPartitioning(); + } + + + + template + const QPartitioning & + QGeneratorBase::get_quadratures() const + { + return q_partitioning; + } + + + + template + void + QGenerator::generate( + const std::vector>> + & level_sets, + const BoundingBox &box, + const unsigned int n_box_splits) + { + std::vector> all_function_bounds; + estimate_function_bounds(level_sets, box, all_function_bounds); + + const std::pair extreme_values = + find_extreme_values(all_function_bounds); + + if (extreme_values.first > this->additional_data.limit_to_be_definite) + { + map_quadrature_to_box(tensor_products[this->q_index], + box, + this->q_partitioning.positive); + } + else if (extreme_values.second < + -(this->additional_data.limit_to_be_definite)) + { + map_quadrature_to_box(tensor_products[this->q_index], + box, + this->q_partitioning.negative); + } + else if (one_positive_one_negative_definite(all_function_bounds)) + { + map_quadrature_to_box(tensor_products[this->q_index], + box, + this->q_partitioning.indefinite); + } + else + { + const std_cxx17::optional data = + find_best_height_direction(all_function_bounds); + + // Check larger than a constant to avoid that min_abs_dfdx is only + // larger by 0 by floating point precision. + if (data && data->min_abs_dfdx > + this->additional_data.lower_bound_implicit_function) + { + create_low_dim_quadratures(data->direction, + level_sets, + box, + n_box_splits); + create_high_dim_quadratures(data->direction, level_sets, box); + } + else if (n_box_splits < this->additional_data.max_box_splits) + { + split_box_and_recurse(level_sets, box, data, n_box_splits); + } + else + { + // We can't split the box recursively forever. Use the midpoint + // method as a last resort. + use_midpoint_method(level_sets, box); + } + } + } + + + + /** + * Return the coordinate direction of the largest side of the box. + * If two or more sides have the same length the returned std::optional + * will be non-set. + */ + template + std_cxx17::optional + direction_of_largest_extent(const BoundingBox &box) + { + // Get the side lengths for each direction and sort them. + std::array, dim> side_lengths; + for (int i = 0; i < dim; i++) + { + side_lengths[i].first = box.side_length(i); + side_lengths[i].second = i; + } + // Sort is lexicographic, so this sorts based on side length first. + std::sort(side_lengths.begin(), side_lengths.end()); + + // Check if the two largest side lengths have the same length. This + // function isn't called in 1D, so the (dim - 2)-element exists. + if (boost::math::epsilon_difference(side_lengths[dim - 1].first, + side_lengths[dim - 2].first) < 100) + return std_cxx17::optional(); + + return side_lengths.back().second; + } + + + + /** + * Return the coordinate direction that the box should be split in, + * assuming that the box should be split it half. + * + * If the box is larger in one coordante direction, this direction is + * returned. If the box have the same extent in all directions, we choose + * the coordinate direction which is closest to being a height-function + * direction. That is, the direction $i$ that has a least negative + * estimate of $|\partial_i \psi_j|$. As a last resort, we choose the + * direction 0, if @p height_direction_data non-set. + */ + template + unsigned int + compute_split_direction( + const BoundingBox & box, + const std_cxx17::optional &height_direction_data) + { + const std_cxx17::optional direction = + direction_of_largest_extent(box); + + if (direction) + return *direction; + + // This direction is closest to being a height direction, so + // we split in this direction. + if (height_direction_data) + return height_direction_data->direction; + + // We have to choose some direction, we might aswell take 0. + return 0; + } + + + + /** + * Split the incoming box in half with respect to the incoming coordinate + * direction and return the left half. + */ + template + inline BoundingBox + left_half(const BoundingBox &box, const unsigned int direction) + { + AssertIndexRange(direction, dim); + + // Move the upper corner half a side-length to the left. + std::pair, Point> corners = box.get_boundary_points(); + corners.second[direction] -= .5 * box.side_length(direction); + + return BoundingBox(corners); + } + + + + /** + * Split the incoming box in half with respect to the incoming coordinate + * direction and return the right half. + */ + template + inline BoundingBox + right_half(const BoundingBox &box, const unsigned int direction) + { + AssertIndexRange(direction, dim); + + // Move the lower corner half a side-length to the right. + std::pair, Point> corners = box.get_boundary_points(); + corners.first[direction] += .5 * box.side_length(direction); + + return BoundingBox(corners); + } + + + + template + void + QGenerator::split_box_and_recurse( + const std::vector>> + & level_sets, + const BoundingBox & box, + const std_cxx17::optional &direction_data, + const unsigned int n_box_splits) + { + if (this->additional_data.split_in_half) + { + const unsigned int direction = + compute_split_direction(box, direction_data); + + const BoundingBox left_box = left_half(box, direction); + const BoundingBox right_box = right_half(box, direction); + + generate(level_sets, left_box, n_box_splits + 1); + generate(level_sets, right_box, n_box_splits + 1); + } + else + { + for (unsigned int i = 0; + i < GeometryInfo::max_children_per_cell; + ++i) + { + generate(level_sets, box.child(i), n_box_splits + 1); + } + } + } + + + + template + void + QGenerator::create_low_dim_quadratures( + const unsigned int height_function_direction, + const std::vector>> + & level_sets, + const BoundingBox &box, + const unsigned int n_box_splits) + { + std::vector> + face_restrictions; + restrict_to_top_and_bottom(level_sets, + box, + height_function_direction, + face_restrictions); + + // We need a vector of references to do the recursive call. + const std::vector>> + restrictions(face_restrictions.begin(), face_restrictions.end()); + + const BoundingBox cross_section = + box.cross_section(height_function_direction); + + low_dim_algorithm.clear_quadratures(); + low_dim_algorithm.generate(restrictions, cross_section, n_box_splits); + } + + + + template + void + QGenerator::create_high_dim_quadratures( + const unsigned int height_function_direction, + const std::vector>> + & level_sets, + const BoundingBox &box) + { + const QPartitioning &low_dim_quadratures = + low_dim_algorithm.get_quadratures(); + + const Quadrature<1> &quadrature1D = + (*this->q_collection1D)[this->q_index]; + + add_tensor_product(low_dim_quadratures.negative, + quadrature1D, + box.lower_bound(height_function_direction), + box.upper_bound(height_function_direction), + height_function_direction, + this->q_partitioning.negative); + + add_tensor_product(low_dim_quadratures.positive, + quadrature1D, + box.lower_bound(height_function_direction), + box.upper_bound(height_function_direction), + height_function_direction, + this->q_partitioning.positive); + + up_through_dimension_creator.generate(level_sets, + box, + low_dim_quadratures.indefinite, + height_function_direction, + this->q_partitioning); + } + + + + template + void + QGenerator::use_midpoint_method( + const std::vector>> + & level_sets, + const BoundingBox &box) + { + const Point center = box.center(); + const Definiteness definiteness = + pointwise_definiteness(level_sets, center); + + ExtendableQuadrature &quadrature = + this->q_partitioning.quadrature_by_definiteness(definiteness); + + quadrature.push_back(center, box.volume()); + } + + + + template + void + QGenerator::set_1D_quadrature(const unsigned int q_index) + { + AssertIndexRange(q_index, this->q_collection1D->size()); + + this->q_index = q_index; + low_dim_algorithm.set_1D_quadrature(q_index); + up_through_dimension_creator.set_1D_quadrature(q_index); + } + + + + template + QGenerator<1, spacedim>::QGenerator( + const hp::QCollection<1> & q_collection1D, + const AdditionalQGeneratorData &additional_data) + : QGeneratorBase<1, spacedim>(q_collection1D, additional_data) + , root_finder( + RootFinder::AdditionalData(additional_data.root_finder_tolerance, + additional_data.max_root_finder_splits)) + { + Assert(q_collection1D.size() > 0, + ExcMessage("Incoming quadrature collection is empty.")); + } + + + + template + void + QGenerator<1, spacedim>::generate( + const std::vector>> + & level_sets, + const BoundingBox<1> &box, + const unsigned int n_box_splits) + { + (void)n_box_splits; + + roots.clear(); + root_finder.find_roots(level_sets, box, roots); + + const Quadrature<1> &quadrature1D = + (*this->q_collection1D)[this->q_index]; + + distribute_points_between_roots(quadrature1D, + box, + roots, + zero_dim_point, + unit_weight, + direction, + level_sets, + this->additional_data, + this->q_partitioning); + + if (spacedim == 1) + this->create_surface_points(level_sets); + } + + + + template + void + QGenerator<1, spacedim>::create_surface_points( + const std::vector>> + &level_sets) + { + Assert(level_sets.size() == 1, ExcInternalError()); + + for (const double root : roots) + { + // A surface integral in 1D is just a point evaluation, + // so the weight is always 1. + const double weight = 1; + const Point<1> point(root); + + Tensor<1, 1> normal = level_sets[0].get().gradient(point); + const double gradient_norm = normal.norm(); + Assert( + gradient_norm > 1e-11, + ExcMessage( + "The level set function has a gradient almost equal to 0.")); + normal *= 1. / gradient_norm; + + this->q_partitioning.surface.push_back(point, weight, normal); + } + } + + + + template + void + QGenerator<1, spacedim>::set_1D_quadrature(const unsigned int q_index) + { + AssertIndexRange(q_index, this->q_collection1D->size()); + this->q_index = q_index; + } + } // namespace QuadratureGeneratorImplementation + } // namespace internal + + using namespace internal::QuadratureGeneratorImplementation; + + + + AdditionalQGeneratorData::AdditionalQGeneratorData( + const unsigned int max_box_splits, + const double lower_bound_implicit_function, + const double min_distance_between_roots, + const double limit_to_be_definite, + const double root_finder_tolerance, + const unsigned int max_root_finder_splits, + bool split_in_half) + : max_box_splits(max_box_splits) + , lower_bound_implicit_function(lower_bound_implicit_function) + , min_distance_between_roots(min_distance_between_roots) + , limit_to_be_definite(limit_to_be_definite) + , root_finder_tolerance(root_finder_tolerance) + , max_root_finder_splits(max_root_finder_splits) + , split_in_half(split_in_half) + {} + + + + template + QuadratureGenerator::QuadratureGenerator( + const hp::QCollection<1> &q_collection, + const AdditionalData & additional_data) + : q_generator(q_collection, additional_data) + { + Assert(q_collection.size() > 0, + ExcMessage("Incoming hp::QCollection<1> is empty.")); + } + + + + template + void + QuadratureGenerator::generate(const Function & level_set, + const BoundingBox &box) + { + Assert(level_set.n_components == 1, + ExcMessage( + "The incoming function should be a scalar level set function," + " it should have one component.")); + Assert(box.volume() > 0, ExcMessage("Incoming box has zero volume.")); + + q_generator.clear_quadratures(); + + std::vector>> level_sets; + level_sets.push_back(level_set); + + const unsigned int n_box_splits = 0; + q_generator.generate(level_sets, box, n_box_splits); + + // With a single level set function, the "indefinite" quadrature should be + // zero. If you call generate() with a ZeroFunction nothing good can be + // done. You will end up here. + Assert( + q_generator.get_quadratures().indefinite.size() == 0, + ExcMessage( + "Generation of quadrature rules failed. This can mean that the level" + "set function is degenerate in some way, e.g. oscillating extremely" + "rapidly.")); + } + + + + template + const Quadrature & + QuadratureGenerator::get_inside_quadrature() const + { + return q_generator.get_quadratures().negative; + } + + + + template + const Quadrature & + QuadratureGenerator::get_outside_quadrature() const + { + return q_generator.get_quadratures().positive; + } + + + + template + const ImmersedSurfaceQuadrature & + QuadratureGenerator::get_surface_quadrature() const + { + return q_generator.get_quadratures().surface; + } + + + template + void + QuadratureGenerator::set_1D_quadrature(const unsigned int q_index) + { + q_generator.set_1D_quadrature(q_index); + } + +} // namespace NonMatching +#include "quadrature_generator.inst" +DEAL_II_NAMESPACE_CLOSE diff --git a/source/non_matching/quadrature_generator.inst.in b/source/non_matching/quadrature_generator.inst.in new file mode 100644 index 0000000000..736ff97ae9 --- /dev/null +++ b/source/non_matching/quadrature_generator.inst.in @@ -0,0 +1,74 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + namespace NonMatching + \{ + template class QuadratureGenerator; + + namespace internal + \{ + namespace QuadratureGeneratorImplementation + \{ + template struct FunctionBounds; + + template std::pair + find_extreme_values( + const std::vector> &); + + template void + estimate_function_bounds( + const std::vector< + std::reference_wrapper>> &, + const BoundingBox &, + std::vector> &); + +// gcc gives a maybe-uninitialized warning in this function when dim = 1, but +// gcc is wrong. We don't need the function when dim = 1, so we avoid +// instantiating it. +#if 1 < deal_II_dimension + template std_cxx17::optional + find_best_height_direction( + const std::vector> &); +#endif + + template void + map_quadrature_to_box(const Quadrature &, + const BoundingBox &, + ExtendableQuadrature &); + + template void + tensor_point_with_1D_quadrature( + const Point &, + const double, + const Quadrature<1> &, + const double, + const double, + const unsigned int, + ExtendableQuadrature &q); + \} + \} + \} + } + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) + { +#if 1 < deal_II_dimension + template class NonMatching::internal::QuadratureGeneratorImplementation:: + UpThroughDimensionCreator; +#endif + } diff --git a/tests/non_matching/find_best_height_direction.cc b/tests/non_matching/find_best_height_direction.cc new file mode 100644 index 0000000000..109f221ea4 --- /dev/null +++ b/tests/non_matching/find_best_height_direction.cc @@ -0,0 +1,106 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test the function find_best_height_direction in + * NonMatching::internal::QuadratureGeneratorImplementation. + */ + +#include + +#include "../tests.h" + + +using namespace dealii; +using namespace NonMatching::internal::QuadratureGeneratorImplementation; + + +// Return a pair with both entries equal to value. +std::pair +pair_with_equal_entries(const double value) +{ + return std::pair(value, value); +} + + + +/* + * Test that find_best_height_direction returns an unset optional if the + * incoming bounds correspond to negative/positive definite functions. + */ +void +test_ignores_definite_functions() +{ + const int dim = 2; + deallog << "test_ignores_definite_functions" << std::endl; + + // Bounds corresponding to one negative and one positive definite function. + std::vector> bounds(2); + bounds[0].value = pair_with_equal_entries(-1); + bounds[1].value = pair_with_equal_entries(1); + + const std_cxx17::optional data = + find_best_height_direction(bounds); + + if (!data) + deallog << "OK" << std::endl; +} + + + +/** + * Create a vector containing two FunctionBounds, set them up so that + * there is one height function direction that is the best. Test that this is + * the direction returned from find_best_height_direction(). + */ +void +test_find_best_height_direction() +{ + deallog << "test_find_best_height_direction" << std::endl; + + const int dim = 2; + + std::vector> bounds(2); + // Set up so that the bounds correspond to indefinite functions. + for (unsigned int i = 0; i < bounds.size(); i++) + { + bounds[i].value.first = -1; + bounds[i].value.second = 1; + } + + // Set up the bounds so that the componenetwise min (over function bounds) + // of the gradient is [3, 5]. This makes 1 the best direction. + bounds[0].gradient[0] = pair_with_equal_entries(3); + bounds[1].gradient[0] = pair_with_equal_entries(4); + bounds[0].gradient[1] = pair_with_equal_entries(6); + bounds[1].gradient[1] = pair_with_equal_entries(5); + + const std_cxx17::optional data = + find_best_height_direction(bounds); + + deallog << "height direction = " << data->direction << std::endl; + deallog << "min_abs_dfdx = " << data->min_abs_dfdx << std::endl; +} + + + +int +main() +{ + initlog(); + test_ignores_definite_functions(); + deallog << std::endl; + test_find_best_height_direction(); +} diff --git a/tests/non_matching/find_best_height_direction.output b/tests/non_matching/find_best_height_direction.output new file mode 100644 index 0000000000..43736c2003 --- /dev/null +++ b/tests/non_matching/find_best_height_direction.output @@ -0,0 +1,7 @@ + +DEAL::test_ignores_definite_functions +DEAL::OK +DEAL:: +DEAL::test_find_best_height_direction +DEAL::height direction = 1 +DEAL::min_abs_dfdx = 5.00000 diff --git a/tests/non_matching/find_extreme_values.cc b/tests/non_matching/find_extreme_values.cc new file mode 100644 index 0000000000..39fea777bd --- /dev/null +++ b/tests/non_matching/find_extreme_values.cc @@ -0,0 +1,85 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test the function find_extreme_values in + * NonMatching::internal::QuadratureGeneratorImplementation. + */ + +#include + +#include "../tests.h" + + +using namespace dealii; +using namespace NonMatching::internal::QuadratureGeneratorImplementation; + + +/** + * Send in a vector with two different function bounds to + * find_extreme_values(). Check that what we get back are the actual extreme + * values. + */ +template +void +test_extreme_values_are_found() +{ + deallog << "test_extreme_values_are_found" << std::endl; + + std::vector> bounds(2); + bounds[0].value.first = 1; + bounds[0].value.second = 2; + bounds[1].value.first = -1; + bounds[1].value.second = 3; + + const std::pair extremes = find_extreme_values(bounds); + + deallog << "min = " << extremes.first << std::endl; + deallog << "max = " << extremes.second << std::endl; +} + + + +/** + * Since the implementation of find_extreme_values() treats the 0th entry + * differently, we check that we get the same entry back if we send in a + * vector with only one entry. + */ +template +void +test_extreme_values_initialized_to_first() +{ + deallog << "test_extreme_values_initialized_to_first" << std::endl; + + std::vector> bounds(1); + bounds[0].value.first = 1; + bounds[0].value.second = 2; + + const std::pair extremes = find_extreme_values(bounds); + + deallog << "min = " << extremes.first << std::endl; + deallog << "max = " << extremes.second << std::endl; +} + + + +int +main() +{ + initlog(); + test_extreme_values_initialized_to_first<1>(); + deallog << std::endl; + test_extreme_values_are_found<1>(); +} diff --git a/tests/non_matching/find_extreme_values.output b/tests/non_matching/find_extreme_values.output new file mode 100644 index 0000000000..3c68e29a14 --- /dev/null +++ b/tests/non_matching/find_extreme_values.output @@ -0,0 +1,8 @@ + +DEAL::test_extreme_values_initialized_to_first +DEAL::min = 1.00000 +DEAL::max = 2.00000 +DEAL:: +DEAL::test_extreme_values_are_found +DEAL::min = -1.00000 +DEAL::max = 3.00000 diff --git a/tests/non_matching/pointwise_definiteness.cc b/tests/non_matching/pointwise_definiteness.cc new file mode 100644 index 0000000000..a8123fa027 --- /dev/null +++ b/tests/non_matching/pointwise_definiteness.cc @@ -0,0 +1,146 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test the function pointwise_definiteness in + * NonMatching::internal::QuadratureGeneratorImplementation. + */ + +#include +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; +using namespace NonMatching::internal::QuadratureGeneratorImplementation; + + +/** + * Call pointwise_definiteness with two positive Functions, + * check that it returns Definiteness::positive. + */ +template +void +test_with_positive_functions() +{ + std::vector> functions; + functions.push_back(Functions::ConstantFunction(1)); + functions.push_back(Functions::ConstantFunction(1)); + + const std::vector>> function_refs( + functions.begin(), functions.end()); + + const Definiteness definiteness = + pointwise_definiteness(function_refs, Point()); + + AssertThrow(definiteness == Definiteness::positive, ExcInternalError()); +} + + + +/** + * Call pointwise_definiteness with two negative Functions, + * check that it returns Definiteness::negative. + */ +template +void +test_with_negative_functions() +{ + std::vector> functions; + functions.push_back(Functions::ConstantFunction(-1)); + functions.push_back(Functions::ConstantFunction(-1)); + + const std::vector>> function_refs( + functions.begin(), functions.end()); + + const Definiteness definiteness = + pointwise_definiteness(function_refs, Point()); + + AssertThrow(definiteness == Definiteness::negative, ExcInternalError()); +} + + + +/** + * Call pointwise_definiteness with with one positive and one negative Function, + * check that it returns Definiteness::indefinite. + */ +template +void +test_with_functions_of_different_sign() +{ + std::vector> functions; + functions.push_back(Functions::ConstantFunction(-1)); + functions.push_back(Functions::ConstantFunction(1)); + + const std::vector>> function_refs( + functions.begin(), functions.end()); + + const Definiteness definiteness = + pointwise_definiteness(function_refs, Point()); + + AssertThrow(definiteness == Definiteness::indefinite, ExcInternalError()); +} + + + +/** + * Call pointwise_definiteness with a single Function which is zero, + * check that it returns Definiteness::indefinite. + * + * This is a special case in the implementation. + */ +template +void +test_first_function_zero() +{ + Functions::ZeroFunction zero_function; + + std::vector>> function_refs; + function_refs.push_back(zero_function); + + const Definiteness definiteness = + pointwise_definiteness(function_refs, Point()); + + AssertThrow(definiteness == Definiteness::indefinite, ExcInternalError()); +} + + + +template +void +run_test() +{ + test_with_positive_functions(); + test_with_negative_functions(); + test_with_functions_of_different_sign(); + test_first_function_zero(); +} + + + +int +main() +{ + initlog(); + run_test<1>(); + run_test<2>(); + run_test<3>(); + deallog << "OK" << std::endl; +} diff --git a/tests/non_matching/pointwise_definiteness.output b/tests/non_matching/pointwise_definiteness.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/non_matching/pointwise_definiteness.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/non_matching/quadrature_generator.cc b/tests/non_matching/quadrature_generator.cc new file mode 100644 index 0000000000..e98a3cc902 --- /dev/null +++ b/tests/non_matching/quadrature_generator.cc @@ -0,0 +1,302 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test the QuadratureGenerator class, by setting up a few simple cuts over the + * unit box and writing the generated quadrature rules to the output file. + * + * Each function beginning with "test_" sets up a level set function and then + * calls the function create_and_print_quadratures() to generate the + * quadratures. + */ + +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + +#include "quadrature_printing.h" + +using NonMatching::QuadratureGenerator; + +/* + * Create immersed quadrature rules over a unit box intersected by the + * incoming level set function. Use a 1D-Gauss quadrature with n_1D_points + * points as a base. Print the constructed quadrature rules to deallog. + */ +template +void +create_and_print_quadratures( + const Function & level_set, + const unsigned int n_1D_points = 2, + const typename QuadratureGenerator::AdditionalData &additional_data = + typename QuadratureGenerator::AdditionalData()) +{ + deallog << "dim=" << dim << std::endl; + + hp::QCollection<1> q_collection; + q_collection.push_back(QGauss<1>(n_1D_points)); + + QuadratureGenerator quadrature_generator(q_collection, additional_data); + + const BoundingBox box = create_unit_bounding_box(); + quadrature_generator.generate(level_set, box); + + deallog << "Inside quadrature" << std::endl; + print_quadrature(quadrature_generator.get_inside_quadrature()); + deallog << "Outside quadrature" << std::endl; + print_quadrature(quadrature_generator.get_outside_quadrature()); + deallog << "Surface quadrature" << std::endl; + print_surface_quadrature(quadrature_generator.get_surface_quadrature()); +} + + + +/* + * Construct level set with a zero contour as a plane cutting straight + * through the unit box. Create and print the constructed quadratures. Do this + * for all unit normals aligned with the coordinate directions. We expect that + * the constructed quadrature has equal number of points in the inside/outside + * region and that they are tensor products. + */ +template +void +test_vertical_cuts_through_center() +{ + deallog << "test_vertical_cuts_through_center" << std::endl; + + Point center; + for (int d = 0; d < dim; ++d) + { + center(d) = .5; + } + for (int direction = 0; direction < dim; ++direction) + { + deallog << "direction=" << direction << std::endl; + const Tensor<1, dim> normal = Point::unit_vector(direction); + const Functions::LevelSet::Plane level_set(center, normal); + create_and_print_quadratures(level_set); + } +} + + + +/* + * Set up a constant negative/positive level set function. Check that + * the constructed quadratures only the inside/outside have points and that + * this quadrature is a tensor product. + */ +template +void +test_constant_level_sets_both_signs() +{ + const Functions::ConstantFunction constant_positive(1); + const Functions::ConstantFunction constant_negative(-1); + deallog << std::endl; + + deallog << "constant_positive" << std::endl; + create_and_print_quadratures(constant_positive); + + deallog << std::endl; + + deallog << "constant_negative" << std::endl; + create_and_print_quadratures(constant_negative); +} + + + +// Set up a level set function corresponding to a plane with normal (1,1) in 2D +// and (1,1,1) in 3D. This makes the inside region a simplex, with vertices +// (0, 0), (0, l), (l, 0), in 2D, +// (0, 0 ,0), (0, 0, l), (0, l, 0), (l, 0, 0), in 3D. +// where l is the edge length. +// This is a good test because we know that the inside weights should sum to the +// area/volume: $V = l^{dim}/dim!$, and that the the surface weights should sum +// to $S = \sqrt(2) l$ in 2D and $S = \sqrt(3)/2 l^2$ in 3D. +template +void +test_simplex_cut() +{ + deallog << "test_simplex_cut" << std::endl; + + const double edge_length = 1. / std::sqrt(2); + + Tensor<1, dim> normal; + for (int i = 0; i < dim; ++i) + normal[i] = 1; + + Point point_in_plane; + point_in_plane[0] = edge_length; + + const Functions::LevelSet::Plane level_set(point_in_plane, normal); + + create_and_print_quadratures(level_set); +} + + + +// Set up a level set function with a zero contour being a plane in the +// direction (1,1) in 2D and (1,1,1) in 3D, such that it cuts the bottom corner +// of the reference cell with a cut of size epsilon. Test that the epsilon cut +// is ignored and we get a tensor product quadrature over the outside region. +template +void +test_epsilon_cut_at_bottom_corner() +{ + deallog << "test_epsilon_cut_at_bottom_corner" << std::endl; + const double epsilon = 1e-15; + Tensor<1, dim> normal; + Point center; + for (int i = 0; i < dim; ++i) + { + normal[i] = 1; + center(i) += epsilon; + } + const Functions::LevelSet::Plane level_set(center, normal); + + create_and_print_quadratures(level_set); +} + + + +/* + * Set up a spherical level set with radius R centered in (0, R) in 2D and + * (0, 0, R) in 3D. The result of this is that the zero contour of the level set + * function cuts exactly through vertex 0 of the unit box + * + * When we choose to split the cell in 4. This test case is difficult for the + * algorithm in 3D. We first get dim - 1 as the first height-direction. But + * after restricting once we get L_a = 0 (where L_a is defined by |\partial_i + * psi| > L_a) for both i = 1,2. Thus we can not choose a second height function + * direction. The results is that the cell is split several times until the + * maximum recursion is reached. When this happens the algorithm uses the + * midpoint method as a fallback. + */ +template +void +test_sphere_cutting_corner_exactly() +{ + deallog << "test_sphere_cutting_corner_exactly" << std::endl; + const double radius = 4; + Point center; + center[dim - 1] = radius; + const Functions::LevelSet::Sphere level_set(center, radius); + + typename QuadratureGenerator::AdditionalData data; + data.split_in_half = false; + data.max_box_splits = 2; + + const unsigned int n_1D_points = 2; + + create_and_print_quadratures(level_set, n_1D_points, data); +} + + + +// A "fake" function used in test_splitting(). This function is constant 1, +// except close to the unit box center, x_i = 0.5, where it has a very large +// Hessian. +template +class ConstantOneButLargeHessianInCenter + : public Functions::ConstantFunction +{ +public: + ConstantOneButLargeHessianInCenter() + : Functions::ConstantFunction(1) + { + for (int d = 0; d < dim; ++d) + unit_box_center(d) = .5; + } + + SymmetricTensor<2, dim> + hessian(const Point &point, const unsigned int) const override + { + SymmetricTensor<2, dim> hessian; + + const double max_distance = 1e-3; + const double diagonal_value = + point.distance(unit_box_center) < max_distance ? 1E3 : 0; + + for (int d = 0; d < dim; ++d) + hessian[d][d] = diagonal_value; + + return hessian; + } + +private: + Point unit_box_center; +}; + + + +// Test the box splitting. Call QuadratureGenerator with a function that is +// constant 1, but has a large Hessian close to the center of the box. This +// should make the algorithm split the box, since the function bounds will be +// large. +template +void +test_splitting() +{ + deallog << "test_splitting" << std::endl; + + const ConstantOneButLargeHessianInCenter level_set; + create_and_print_quadratures(level_set); +} + + + +// Some of the test cases only make sense for a given dimension, +// so we list the cases for each dimension. +int +main() +{ + initlog(); + // 1D + test_vertical_cuts_through_center<1>(); + deallog << std::endl; + // 2D + test_vertical_cuts_through_center<2>(); + deallog << std::endl; + test_constant_level_sets_both_signs<2>(); + deallog << std::endl; + test_simplex_cut<2>(); + deallog << std::endl; + test_epsilon_cut_at_bottom_corner<2>(); + deallog << std::endl; + test_sphere_cutting_corner_exactly<2>(); + deallog << std::endl; + test_splitting<2>(); + deallog << std::endl; + // 3D + test_vertical_cuts_through_center<3>(); + deallog << std::endl; + test_simplex_cut<3>(); + deallog << std::endl; + test_epsilon_cut_at_bottom_corner<3>(); + deallog << std::endl; + test_sphere_cutting_corner_exactly<3>(); +} diff --git a/tests/non_matching/quadrature_generator.output b/tests/non_matching/quadrature_generator.output new file mode 100644 index 0000000000..20f513046b --- /dev/null +++ b/tests/non_matching/quadrature_generator.output @@ -0,0 +1,389 @@ + +DEAL::test_vertical_cuts_through_center +DEAL::direction=0 +DEAL::dim=1 +DEAL::Inside quadrature +DEAL::0.105662, 0.250000 +DEAL::0.394338, 0.250000 +DEAL::Outside quadrature +DEAL::0.605662, 0.250000 +DEAL::0.894338, 0.250000 +DEAL::Surface quadrature +DEAL::0.500000, 1.00000, 1.00000 +DEAL:: +DEAL::test_vertical_cuts_through_center +DEAL::direction=0 +DEAL::dim=2 +DEAL::Inside quadrature +DEAL::0.105662, 0.211325, 0.125000 +DEAL::0.394338, 0.211325, 0.125000 +DEAL::0.105662, 0.788675, 0.125000 +DEAL::0.394338, 0.788675, 0.125000 +DEAL::Outside quadrature +DEAL::0.605662, 0.211325, 0.125000 +DEAL::0.894338, 0.211325, 0.125000 +DEAL::0.605662, 0.788675, 0.125000 +DEAL::0.894338, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL::0.500000, 0.211325, 0.500000, 1.00000, 0.00000 +DEAL::0.500000, 0.788675, 0.500000, 1.00000, 0.00000 +DEAL::direction=1 +DEAL::dim=2 +DEAL::Inside quadrature +DEAL::0.211325, 0.105662, 0.125000 +DEAL::0.211325, 0.394338, 0.125000 +DEAL::0.788675, 0.105662, 0.125000 +DEAL::0.788675, 0.394338, 0.125000 +DEAL::Outside quadrature +DEAL::0.211325, 0.605662, 0.125000 +DEAL::0.211325, 0.894338, 0.125000 +DEAL::0.788675, 0.605662, 0.125000 +DEAL::0.788675, 0.894338, 0.125000 +DEAL::Surface quadrature +DEAL::0.211325, 0.500000, 0.500000, 0.00000, 1.00000 +DEAL::0.788675, 0.500000, 0.500000, 0.00000, 1.00000 +DEAL:: +DEAL:: +DEAL::constant_positive +DEAL::dim=2 +DEAL::Inside quadrature +DEAL::Outside quadrature +DEAL::0.211325, 0.211325, 0.250000 +DEAL::0.788675, 0.211325, 0.250000 +DEAL::0.211325, 0.788675, 0.250000 +DEAL::0.788675, 0.788675, 0.250000 +DEAL::Surface quadrature +DEAL:: +DEAL::constant_negative +DEAL::dim=2 +DEAL::Inside quadrature +DEAL::0.211325, 0.211325, 0.250000 +DEAL::0.788675, 0.211325, 0.250000 +DEAL::0.211325, 0.788675, 0.250000 +DEAL::0.788675, 0.788675, 0.250000 +DEAL::Outside quadrature +DEAL::Surface quadrature +DEAL:: +DEAL::test_simplex_cut +DEAL::dim=2 +DEAL::Inside quadrature +DEAL::0.117851, 0.149429, 0.0985844 +DEAL::0.439826, 0.149429, 0.0985844 +DEAL::0.0315781, 0.557678, 0.0264156 +DEAL::0.117851, 0.557678, 0.0264156 +DEAL::Outside quadrature +DEAL::0.211325, 0.769002, 0.0732233 +DEAL::0.788675, 0.769002, 0.0732233 +DEAL::0.211325, 0.938104, 0.0732233 +DEAL::0.788675, 0.938104, 0.0732233 +DEAL::0.651151, 0.149429, 0.0781923 +DEAL::0.906526, 0.149429, 0.0781923 +DEAL::0.329176, 0.557678, 0.150361 +DEAL::0.820253, 0.557678, 0.150361 +DEAL::Surface quadrature +DEAL::0.557678, 0.149429, 0.500000, 0.707107, 0.707107 +DEAL::0.149429, 0.557678, 0.500000, 0.707107, 0.707107 +DEAL:: +DEAL::test_epsilon_cut_at_bottom_corner +DEAL::dim=2 +DEAL::Inside quadrature +DEAL::Outside quadrature +DEAL::0.211325, 0.211325, 0.250000 +DEAL::0.788675, 0.211325, 0.250000 +DEAL::0.211325, 0.788675, 0.250000 +DEAL::0.788675, 0.788675, 0.250000 +DEAL::Surface quadrature +DEAL:: +DEAL::test_sphere_cutting_corner_exactly +DEAL::dim=2 +DEAL::Inside quadrature +DEAL::0.211325, 0.215731, 0.248603 +DEAL::0.211325, 0.789856, 0.248603 +DEAL::0.788675, 0.273253, 0.230370 +DEAL::0.788675, 0.805269, 0.230370 +DEAL::Outside quadrature +DEAL::0.211325, 0.00118050, 0.00139654 +DEAL::0.211325, 0.00440568, 0.00139654 +DEAL::0.788675, 0.0165936, 0.0196304 +DEAL::0.788675, 0.0619282, 0.0196304 +DEAL::Surface quadrature +DEAL::0.211325, 0.00558618, 0.500699, 0.0528312, -0.998603 +DEAL::0.788675, 0.0785218, 0.510012, 0.197169, -0.980370 +DEAL:: +DEAL::test_splitting +DEAL::dim=2 +DEAL::Inside quadrature +DEAL::Outside quadrature +DEAL::0.105662, 0.211325, 0.125000 +DEAL::0.394338, 0.211325, 0.125000 +DEAL::0.105662, 0.788675, 0.125000 +DEAL::0.394338, 0.788675, 0.125000 +DEAL::0.605662, 0.211325, 0.125000 +DEAL::0.894338, 0.211325, 0.125000 +DEAL::0.605662, 0.788675, 0.125000 +DEAL::0.894338, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL:: +DEAL::test_vertical_cuts_through_center +DEAL::direction=0 +DEAL::dim=3 +DEAL::Inside quadrature +DEAL::0.105662, 0.211325, 0.211325, 0.0625000 +DEAL::0.394338, 0.211325, 0.211325, 0.0625000 +DEAL::0.105662, 0.788675, 0.211325, 0.0625000 +DEAL::0.394338, 0.788675, 0.211325, 0.0625000 +DEAL::0.105662, 0.211325, 0.788675, 0.0625000 +DEAL::0.394338, 0.211325, 0.788675, 0.0625000 +DEAL::0.105662, 0.788675, 0.788675, 0.0625000 +DEAL::0.394338, 0.788675, 0.788675, 0.0625000 +DEAL::Outside quadrature +DEAL::0.605662, 0.211325, 0.211325, 0.0625000 +DEAL::0.894338, 0.211325, 0.211325, 0.0625000 +DEAL::0.605662, 0.788675, 0.211325, 0.0625000 +DEAL::0.894338, 0.788675, 0.211325, 0.0625000 +DEAL::0.605662, 0.211325, 0.788675, 0.0625000 +DEAL::0.894338, 0.211325, 0.788675, 0.0625000 +DEAL::0.605662, 0.788675, 0.788675, 0.0625000 +DEAL::0.894338, 0.788675, 0.788675, 0.0625000 +DEAL::Surface quadrature +DEAL::0.500000, 0.211325, 0.211325, 0.250000, 1.00000, 0.00000, 0.00000 +DEAL::0.500000, 0.788675, 0.211325, 0.250000, 1.00000, 0.00000, 0.00000 +DEAL::0.500000, 0.211325, 0.788675, 0.250000, 1.00000, 0.00000, 0.00000 +DEAL::0.500000, 0.788675, 0.788675, 0.250000, 1.00000, 0.00000, 0.00000 +DEAL::direction=1 +DEAL::dim=3 +DEAL::Inside quadrature +DEAL::0.211325, 0.105662, 0.211325, 0.0625000 +DEAL::0.211325, 0.394338, 0.211325, 0.0625000 +DEAL::0.211325, 0.105662, 0.788675, 0.0625000 +DEAL::0.211325, 0.394338, 0.788675, 0.0625000 +DEAL::0.788675, 0.105662, 0.211325, 0.0625000 +DEAL::0.788675, 0.394338, 0.211325, 0.0625000 +DEAL::0.788675, 0.105662, 0.788675, 0.0625000 +DEAL::0.788675, 0.394338, 0.788675, 0.0625000 +DEAL::Outside quadrature +DEAL::0.211325, 0.605662, 0.211325, 0.0625000 +DEAL::0.211325, 0.894338, 0.211325, 0.0625000 +DEAL::0.211325, 0.605662, 0.788675, 0.0625000 +DEAL::0.211325, 0.894338, 0.788675, 0.0625000 +DEAL::0.788675, 0.605662, 0.211325, 0.0625000 +DEAL::0.788675, 0.894338, 0.211325, 0.0625000 +DEAL::0.788675, 0.605662, 0.788675, 0.0625000 +DEAL::0.788675, 0.894338, 0.788675, 0.0625000 +DEAL::Surface quadrature +DEAL::0.211325, 0.500000, 0.211325, 0.250000, 0.00000, 1.00000, 0.00000 +DEAL::0.211325, 0.500000, 0.788675, 0.250000, 0.00000, 1.00000, 0.00000 +DEAL::0.788675, 0.500000, 0.211325, 0.250000, 0.00000, 1.00000, 0.00000 +DEAL::0.788675, 0.500000, 0.788675, 0.250000, 0.00000, 1.00000, 0.00000 +DEAL::direction=2 +DEAL::dim=3 +DEAL::Inside quadrature +DEAL::0.211325, 0.211325, 0.105662, 0.0625000 +DEAL::0.211325, 0.211325, 0.394338, 0.0625000 +DEAL::0.788675, 0.211325, 0.105662, 0.0625000 +DEAL::0.788675, 0.211325, 0.394338, 0.0625000 +DEAL::0.211325, 0.788675, 0.105662, 0.0625000 +DEAL::0.211325, 0.788675, 0.394338, 0.0625000 +DEAL::0.788675, 0.788675, 0.105662, 0.0625000 +DEAL::0.788675, 0.788675, 0.394338, 0.0625000 +DEAL::Outside quadrature +DEAL::0.211325, 0.211325, 0.605662, 0.0625000 +DEAL::0.211325, 0.211325, 0.894338, 0.0625000 +DEAL::0.788675, 0.211325, 0.605662, 0.0625000 +DEAL::0.788675, 0.211325, 0.894338, 0.0625000 +DEAL::0.211325, 0.788675, 0.605662, 0.0625000 +DEAL::0.211325, 0.788675, 0.894338, 0.0625000 +DEAL::0.788675, 0.788675, 0.605662, 0.0625000 +DEAL::0.788675, 0.788675, 0.894338, 0.0625000 +DEAL::Surface quadrature +DEAL::0.211325, 0.211325, 0.500000, 0.250000, 0.00000, 0.00000, 1.00000 +DEAL::0.788675, 0.211325, 0.500000, 0.250000, 0.00000, 0.00000, 1.00000 +DEAL::0.211325, 0.788675, 0.500000, 0.250000, 0.00000, 0.00000, 1.00000 +DEAL::0.788675, 0.788675, 0.500000, 0.250000, 0.00000, 0.00000, 1.00000 +DEAL:: +DEAL::test_simplex_cut +DEAL::dim=3 +DEAL::Inside quadrature +DEAL::0.0929463, 0.117851, 0.149429, 0.0216800 +DEAL::0.346880, 0.117851, 0.149429, 0.0216800 +DEAL::0.0249049, 0.439826, 0.149429, 0.00580914 +DEAL::0.0929463, 0.439826, 0.149429, 0.00580914 +DEAL::0.0249049, 0.0315781, 0.557678, 0.00155655 +DEAL::0.0929463, 0.0315781, 0.557678, 0.00155655 +DEAL::0.00667324, 0.117851, 0.557678, 0.000417078 +DEAL::0.0249049, 0.117851, 0.557678, 0.000417078 +DEAL::Outside quadrature +DEAL::0.211325, 0.211325, 0.769002, 0.0366117 +DEAL::0.788675, 0.211325, 0.769002, 0.0366117 +DEAL::0.211325, 0.788675, 0.769002, 0.0366117 +DEAL::0.788675, 0.788675, 0.769002, 0.0366117 +DEAL::0.211325, 0.211325, 0.938104, 0.0366117 +DEAL::0.788675, 0.211325, 0.938104, 0.0366117 +DEAL::0.211325, 0.788675, 0.938104, 0.0366117 +DEAL::0.788675, 0.788675, 0.938104, 0.0366117 +DEAL::0.211325, 0.651151, 0.149429, 0.0390962 +DEAL::0.788675, 0.651151, 0.149429, 0.0390962 +DEAL::0.211325, 0.906526, 0.149429, 0.0390962 +DEAL::0.788675, 0.906526, 0.149429, 0.0390962 +DEAL::0.211325, 0.329176, 0.557678, 0.0751805 +DEAL::0.788675, 0.329176, 0.557678, 0.0751805 +DEAL::0.211325, 0.820253, 0.557678, 0.0751805 +DEAL::0.788675, 0.820253, 0.557678, 0.0751805 +DEAL::0.558205, 0.117851, 0.149429, 0.0276122 +DEAL::0.881621, 0.117851, 0.149429, 0.0276122 +DEAL::0.304271, 0.439826, 0.149429, 0.0434831 +DEAL::0.813580, 0.439826, 0.149429, 0.0434831 +DEAL::0.304271, 0.0315781, 0.557678, 0.0116512 +DEAL::0.813580, 0.0315781, 0.557678, 0.0116512 +DEAL::0.236230, 0.117851, 0.557678, 0.0127907 +DEAL::0.795348, 0.117851, 0.557678, 0.0127907 +DEAL::Surface quadrature +DEAL::0.439826, 0.117851, 0.149429, 0.170753, 0.577350, 0.577350, 0.577350 +DEAL::0.117851, 0.439826, 0.149429, 0.170753, 0.577350, 0.577350, 0.577350 +DEAL::0.117851, 0.0315781, 0.557678, 0.0457532, 0.577350, 0.577350, 0.577350 +DEAL::0.0315781, 0.117851, 0.557678, 0.0457532, 0.577350, 0.577350, 0.577350 +DEAL:: +DEAL::test_epsilon_cut_at_bottom_corner +DEAL::dim=3 +DEAL::Inside quadrature +DEAL::Outside quadrature +DEAL::0.211325, 0.211325, 0.211325, 0.125000 +DEAL::0.788675, 0.211325, 0.211325, 0.125000 +DEAL::0.211325, 0.788675, 0.211325, 0.125000 +DEAL::0.788675, 0.788675, 0.211325, 0.125000 +DEAL::0.211325, 0.211325, 0.788675, 0.125000 +DEAL::0.788675, 0.211325, 0.788675, 0.125000 +DEAL::0.211325, 0.788675, 0.788675, 0.125000 +DEAL::0.788675, 0.788675, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL:: +DEAL::test_sphere_cutting_corner_exactly +DEAL::dim=3 +DEAL::Inside quadrature +DEAL::0.125000, 0.125000, 0.214407, 0.0311279 +DEAL::0.125000, 0.125000, 0.789501, 0.0311279 +DEAL::0.302831, 0.0528312, 0.220655, 0.00772008 +DEAL::0.302831, 0.0528312, 0.791175, 0.00772008 +DEAL::0.447169, 0.0528312, 0.231377, 0.00761387 +DEAL::0.447169, 0.0528312, 0.794048, 0.00761387 +DEAL::0.302831, 0.197169, 0.224225, 0.00768472 +DEAL::0.302831, 0.197169, 0.792132, 0.00768472 +DEAL::0.447169, 0.197169, 0.234959, 0.00757839 +DEAL::0.447169, 0.197169, 0.795008, 0.00757839 +DEAL::0.0528312, 0.302831, 0.220655, 0.00772008 +DEAL::0.0528312, 0.302831, 0.791175, 0.00772008 +DEAL::0.197169, 0.302831, 0.224225, 0.00768472 +DEAL::0.197169, 0.302831, 0.792132, 0.00768472 +DEAL::0.0528312, 0.447169, 0.231377, 0.00761387 +DEAL::0.0528312, 0.447169, 0.794048, 0.00761387 +DEAL::0.197169, 0.447169, 0.234959, 0.00757839 +DEAL::0.197169, 0.447169, 0.795008, 0.00757839 +DEAL::0.302831, 0.302831, 0.229459, 0.00763287 +DEAL::0.302831, 0.302831, 0.793534, 0.00763287 +DEAL::0.447169, 0.302831, 0.240211, 0.00752636 +DEAL::0.447169, 0.302831, 0.796415, 0.00752636 +DEAL::0.302831, 0.447169, 0.240211, 0.00752636 +DEAL::0.302831, 0.447169, 0.796415, 0.00752636 +DEAL::0.447169, 0.447169, 0.251000, 0.00741948 +DEAL::0.447169, 0.447169, 0.799306, 0.00741948 +DEAL::0.605662, 0.105662, 0.248812, 0.0297646 +DEAL::0.605662, 0.105662, 0.798720, 0.0297646 +DEAL::0.894338, 0.105662, 0.292317, 0.0280408 +DEAL::0.894338, 0.105662, 0.810377, 0.0280408 +DEAL::0.605662, 0.394338, 0.263246, 0.0291927 +DEAL::0.605662, 0.394338, 0.802587, 0.0291927 +DEAL::0.894338, 0.394338, 0.306956, 0.0274608 +DEAL::0.894338, 0.394338, 0.814299, 0.0274608 +DEAL::0.105662, 0.605662, 0.248812, 0.0297646 +DEAL::0.105662, 0.605662, 0.798720, 0.0297646 +DEAL::0.394338, 0.605662, 0.263246, 0.0291927 +DEAL::0.394338, 0.605662, 0.802587, 0.0291927 +DEAL::0.105662, 0.894338, 0.292317, 0.0280408 +DEAL::0.105662, 0.894338, 0.810377, 0.0280408 +DEAL::0.394338, 0.894338, 0.306956, 0.0274608 +DEAL::0.394338, 0.894338, 0.814299, 0.0274608 +DEAL::0.605662, 0.605662, 0.284500, 0.0283505 +DEAL::0.605662, 0.605662, 0.808282, 0.0283505 +DEAL::0.894338, 0.605662, 0.328517, 0.0266065 +DEAL::0.894338, 0.605662, 0.820077, 0.0266065 +DEAL::0.605662, 0.894338, 0.328517, 0.0266065 +DEAL::0.605662, 0.894338, 0.820077, 0.0266065 +DEAL::0.894338, 0.894338, 0.373180, 0.0248367 +DEAL::0.894338, 0.894338, 0.832044, 0.0248367 +DEAL::Outside quadrature +DEAL::0.125000, 0.125000, 0.000825891, 0.000122130 +DEAL::0.125000, 0.125000, 0.00308227, 0.000122130 +DEAL::0.302831, 0.0528312, 0.00249992, 9.24198e-05 +DEAL::0.302831, 0.0528312, 0.00932981, 9.24198e-05 +DEAL::0.447169, 0.0528312, 0.00537287, 0.000198630 +DEAL::0.447169, 0.0528312, 0.0200518, 0.000198630 +DEAL::0.302831, 0.197169, 0.00345648, 0.000127783 +DEAL::0.302831, 0.197169, 0.0128997, 0.000127783 +DEAL::0.447169, 0.197169, 0.00633270, 0.000234115 +DEAL::0.447169, 0.197169, 0.0236340, 0.000234115 +DEAL::0.0528312, 0.302831, 0.00249992, 9.24198e-05 +DEAL::0.0528312, 0.302831, 0.00932981, 9.24198e-05 +DEAL::0.197169, 0.302831, 0.00345648, 0.000127783 +DEAL::0.197169, 0.302831, 0.0128997, 0.000127783 +DEAL::0.0528312, 0.447169, 0.00537287, 0.000198630 +DEAL::0.0528312, 0.447169, 0.0200518, 0.000198630 +DEAL::0.197169, 0.447169, 0.00633270, 0.000234115 +DEAL::0.197169, 0.447169, 0.0236340, 0.000234115 +DEAL::0.302831, 0.302831, 0.00485894, 0.000179631 +DEAL::0.302831, 0.302831, 0.0181338, 0.000179631 +DEAL::0.447169, 0.302831, 0.00773999, 0.000286141 +DEAL::0.447169, 0.302831, 0.0288860, 0.000286141 +DEAL::0.302831, 0.447169, 0.00773999, 0.000286141 +DEAL::0.302831, 0.447169, 0.0288860, 0.000286141 +DEAL::0.447169, 0.447169, 0.0106310, 0.000393018 +DEAL::0.447169, 0.447169, 0.0396753, 0.000393018 +DEAL::0.605662, 0.105662, 0.0100446, 0.00148535 +DEAL::0.605662, 0.105662, 0.0374868, 0.00148535 +DEAL::0.894338, 0.105662, 0.0217017, 0.00320918 +DEAL::0.894338, 0.105662, 0.0809920, 0.00320918 +DEAL::0.605662, 0.394338, 0.0139121, 0.00205728 +DEAL::0.605662, 0.394338, 0.0519208, 0.00205728 +DEAL::0.894338, 0.394338, 0.0256243, 0.00378924 +DEAL::0.894338, 0.394338, 0.0956313, 0.00378924 +DEAL::0.105662, 0.605662, 0.0100446, 0.00148535 +DEAL::0.105662, 0.605662, 0.0374868, 0.00148535 +DEAL::0.394338, 0.605662, 0.0139121, 0.00205728 +DEAL::0.394338, 0.605662, 0.0519208, 0.00205728 +DEAL::0.105662, 0.894338, 0.0217017, 0.00320918 +DEAL::0.105662, 0.894338, 0.0809920, 0.00320918 +DEAL::0.394338, 0.894338, 0.0256243, 0.00378924 +DEAL::0.394338, 0.894338, 0.0956313, 0.00378924 +DEAL::0.605662, 0.605662, 0.0196073, 0.00289946 +DEAL::0.605662, 0.605662, 0.0731755, 0.00289946 +DEAL::0.894338, 0.605662, 0.0314015, 0.00464354 +DEAL::0.894338, 0.605662, 0.117192, 0.00464354 +DEAL::0.605662, 0.894338, 0.0314015, 0.00464354 +DEAL::0.605662, 0.894338, 0.117192, 0.00464354 +DEAL::0.894338, 0.894338, 0.0433691, 0.00641327 +DEAL::0.894338, 0.894338, 0.161856, 0.00641327 +DEAL::Surface quadrature +DEAL::0.125000, 0.125000, 0.00390816, 0.0625611, 0.0312500, 0.0312500, -0.999023 +DEAL::0.302831, 0.0528312, 0.0118297, 0.0156713, 0.0757078, 0.0132078, -0.997043 +DEAL::0.447169, 0.0528312, 0.0254247, 0.0157250, 0.111792, 0.0132078, -0.993644 +DEAL::0.302831, 0.197169, 0.0163562, 0.0156892, 0.0757078, 0.0492922, -0.995911 +DEAL::0.447169, 0.197169, 0.0299667, 0.0157429, 0.111792, 0.0492922, -0.992508 +DEAL::0.0528312, 0.302831, 0.0118297, 0.0156713, 0.0132078, 0.0757078, -0.997043 +DEAL::0.197169, 0.302831, 0.0163562, 0.0156892, 0.0492922, 0.0757078, -0.995911 +DEAL::0.0528312, 0.447169, 0.0254247, 0.0157250, 0.0132078, 0.111792, -0.993644 +DEAL::0.197169, 0.447169, 0.0299667, 0.0157429, 0.0492922, 0.111792, -0.992508 +DEAL::0.302831, 0.302831, 0.0229928, 0.0157153, 0.0757078, 0.0757078, -0.994252 +DEAL::0.447169, 0.302831, 0.0366260, 0.0157694, 0.111792, 0.0757078, -0.990843 +DEAL::0.302831, 0.447169, 0.0366260, 0.0157694, 0.0757078, 0.111792, -0.990843 +DEAL::0.447169, 0.447169, 0.0503063, 0.0158240, 0.111792, 0.111792, -0.987423 +DEAL::0.605662, 0.105662, 0.0475313, 0.0632516, 0.151416, 0.0264156, -0.988117 +DEAL::0.894338, 0.105662, 0.102694, 0.0641469, 0.223584, 0.0264156, -0.974327 +DEAL::0.605662, 0.394338, 0.0658329, 0.0635459, 0.151416, 0.0985844, -0.983542 +DEAL::0.894338, 0.394338, 0.121256, 0.0644538, 0.223584, 0.0985844, -0.969686 +DEAL::0.105662, 0.605662, 0.0475313, 0.0632516, 0.0264156, 0.151416, -0.988117 +DEAL::0.394338, 0.605662, 0.0658329, 0.0635459, 0.0985844, 0.151416, -0.983542 +DEAL::0.105662, 0.894338, 0.102694, 0.0641469, 0.0264156, 0.223584, -0.974327 +DEAL::0.394338, 0.894338, 0.121256, 0.0644538, 0.0985844, 0.223584, -0.969686 +DEAL::0.605662, 0.605662, 0.0927828, 0.0639842, 0.151416, 0.151416, -0.976804 +DEAL::0.894338, 0.605662, 0.148593, 0.0649113, 0.223584, 0.151416, -0.962852 +DEAL::0.605662, 0.894338, 0.148593, 0.0649113, 0.151416, 0.223584, -0.962852 +DEAL::0.894338, 0.894338, 0.205225, 0.0658801, 0.223584, 0.223584, -0.948694 diff --git a/tests/non_matching/quadrature_generator_clears_between_calls.cc b/tests/non_matching/quadrature_generator_clears_between_calls.cc new file mode 100644 index 0000000000..ee45b232a3 --- /dev/null +++ b/tests/non_matching/quadrature_generator_clears_between_calls.cc @@ -0,0 +1,88 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test that QuadratureGenerator clears its previously created quadratures + * when we call generate() again. + */ + +#include +#include + +#include + +#include "deal.II/non_matching/quadrature_generator.h" + +#include "../tests.h" + + +// Print the sizes of all the quadratures that QuadratureGenerator creates to +// deallog. +template +void +print_n_quadrature_points( + const NonMatching::QuadratureGenerator &quadrature_generator) +{ + deallog << "inside " << quadrature_generator.get_inside_quadrature().size() + << std::endl; + + deallog << "outside " << quadrature_generator.get_outside_quadrature().size() + << std::endl; + + deallog << "surface " << quadrature_generator.get_surface_quadrature().size() + << std::endl; + + deallog << std::endl; +} + + + +// Call the QuadratureGenerator::generate with the same level set function +// twice. Make sure that the sizes of the constructed quadratures are the same +// both times. The purpose is to make sure that the previously created +// quadratures have been cleared before we create the new ones. +template +void +test() +{ + deallog << "dim = " << dim << std::endl; + + const hp::QCollection<1> q_collection(QGauss<1>(1)); + NonMatching::QuadratureGenerator quadrature_generator(q_collection); + + const Functions::ConstantFunction level_set(1); + + const BoundingBox box = create_unit_bounding_box(); + + deallog << "quadrature sizes first call" << std::endl; + quadrature_generator.generate(level_set, box); + print_n_quadrature_points(quadrature_generator); + + deallog << "quadrature sizes second call" << std::endl; + quadrature_generator.generate(level_set, box); + print_n_quadrature_points(quadrature_generator); +} + + + +int +main() +{ + initlog(); + + test<1>(); + test<2>(); + test<3>(); +} diff --git a/tests/non_matching/quadrature_generator_clears_between_calls.output b/tests/non_matching/quadrature_generator_clears_between_calls.output new file mode 100644 index 0000000000..b12c146d76 --- /dev/null +++ b/tests/non_matching/quadrature_generator_clears_between_calls.output @@ -0,0 +1,34 @@ + +DEAL::dim = 1 +DEAL::quadrature sizes first call +DEAL::inside 0 +DEAL::outside 1 +DEAL::surface 0 +DEAL:: +DEAL::quadrature sizes second call +DEAL::inside 0 +DEAL::outside 1 +DEAL::surface 0 +DEAL:: +DEAL::dim = 2 +DEAL::quadrature sizes first call +DEAL::inside 0 +DEAL::outside 1 +DEAL::surface 0 +DEAL:: +DEAL::quadrature sizes second call +DEAL::inside 0 +DEAL::outside 1 +DEAL::surface 0 +DEAL:: +DEAL::dim = 3 +DEAL::quadrature sizes first call +DEAL::inside 0 +DEAL::outside 1 +DEAL::surface 0 +DEAL:: +DEAL::quadrature sizes second call +DEAL::inside 0 +DEAL::outside 1 +DEAL::surface 0 +DEAL:: diff --git a/tests/non_matching/quadrature_generator_sphere.cc b/tests/non_matching/quadrature_generator_sphere.cc new file mode 100644 index 0000000000..4ed0361cbe --- /dev/null +++ b/tests/non_matching/quadrature_generator_sphere.cc @@ -0,0 +1,101 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +/** + * Compute the volume and surface area of a ball/sphere by setting up a + * level set function immersed in a background mesh, generating the + * quadrature rules, and summing the weights. + */ +template +void +calculate_volume_and_surface_area() +{ + // Set up a background mesh + Triangulation triangulation; + const int n_subdivisions = 12; + const double gridsize = 2.07; + GridGenerator::subdivided_hyper_cube(triangulation, + n_subdivisions, + -gridsize / 2, + gridsize / 2); + + + // Description of the immersed domain. + const Functions::LevelSet::Sphere level_set; + + // Create a quadrature generator. + const hp::QCollection<1> q_collection1D(QGauss<1>(2)); + NonMatching::QuadratureGenerator quadrature_generator(q_collection1D); + + // Go over all cells and compute the volume and surface area. + double surface_area = 0, volume = 0; + for (const auto cell : triangulation.active_cell_iterators()) + { + // Create a box corresponding to the cell. + std::pair, Point> lower_upper_corner; + lower_upper_corner.first = cell->vertex(0); + lower_upper_corner.second = + cell->vertex(GeometryInfo::vertices_per_cell - 1); + const BoundingBox box(lower_upper_corner); + + // Generate immersed quadrature rules. + quadrature_generator.generate(level_set, box); + + // Get the quadrature rules. + const Quadrature &inside_quadrature = + quadrature_generator.get_inside_quadrature(); + const NonMatching::ImmersedSurfaceQuadrature &surface_quadrature = + quadrature_generator.get_surface_quadrature(); + + // Sum the weights to get the area/volume of the sphere. + for (unsigned int i = 0; i < inside_quadrature.size(); ++i) + volume += inside_quadrature.weight(i); + + // Sum the weights to get the circumference/surface area of the sphere. + for (unsigned int i = 0; i < surface_quadrature.size(); ++i) + surface_area += surface_quadrature.weight(i); + } + + deallog << "dim = " << dim << std::endl; + deallog << (2 == dim ? "area = " : "volume = "); + deallog << volume / M_PI << " * pi" << std::endl; + deallog << (2 == dim ? "circumference = " : "surface area = "); + deallog << surface_area / M_PI << " * pi" << std::endl; + deallog << std::endl; +} + + +int +main() +{ + initlog(); + calculate_volume_and_surface_area<2>(); + calculate_volume_and_surface_area<3>(); +} diff --git a/tests/non_matching/quadrature_generator_sphere.output b/tests/non_matching/quadrature_generator_sphere.output new file mode 100644 index 0000000000..caf1f3fa93 --- /dev/null +++ b/tests/non_matching/quadrature_generator_sphere.output @@ -0,0 +1,9 @@ + +DEAL::dim = 2 +DEAL::area = 1.00000 * pi +DEAL::circumference = 1.99997 * pi +DEAL:: +DEAL::dim = 3 +DEAL::volume = 1.33334 * pi +DEAL::surface area = 3.99996 * pi +DEAL:: diff --git a/tests/non_matching/quadrature_printing.h b/tests/non_matching/quadrature_printing.h new file mode 100644 index 0000000000..004c8fccbf --- /dev/null +++ b/tests/non_matching/quadrature_printing.h @@ -0,0 +1,70 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_quadrature_printing_h_ +#define dealii_quadrature_printing_h_ + +#include +#include + +#include + +using namespace dealii; + +/* + * Print the incoming quadrature to deallog as comma separated values: + * point[0], ..., point[dim-1], weight + */ +template +void +print_quadrature(const Quadrature &quadrature) +{ + for (unsigned int i = 0; i < quadrature.size(); ++i) + { + const Point &point = quadrature.point(i); + for (int d = 0; d < dim; d++) + deallog << point[d] << ", "; + + deallog << quadrature.weight(i) << std::endl; + } +} + + + +/* + * Print the incoming surface quadrature to deallog as comma separated values: + * p[0], ..., p[dim-1], weight, normal[0], ..., normal[dim-1] + */ +template +void +print_surface_quadrature( + const NonMatching::ImmersedSurfaceQuadrature &quadrature) +{ + for (unsigned int i = 0; i < quadrature.size(); ++i) + { + const Point &point = quadrature.point(i); + for (int d = 0; d < dim; d++) + deallog << point[d] << ", "; + + deallog << quadrature.weight(i); + + const Tensor<1, dim> &normal = quadrature.normal_vector(i); + for (int d = 0; d < dim; d++) + deallog << ", " << normal[d]; + deallog << std::endl; + } +} + +#endif diff --git a/tests/non_matching/root_finder.cc b/tests/non_matching/root_finder.cc new file mode 100644 index 0000000000..5b1f98cda0 --- /dev/null +++ b/tests/non_matching/root_finder.cc @@ -0,0 +1,160 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test the RootFinder class in internal::QuadratureGeneratorImplementation. + */ + +#include +#include +#include +#include +#include +#include + +#include + +#include "../tests.h" + + +using namespace dealii; +using namespace NonMatching::internal::QuadratureGeneratorImplementation; + +// Use RootFinder to find the roots of the incoming functions over the interval +// [0, 1]. Print the roots to deallog. +void +find_and_print_roots( + const std::vector>> &functions) +{ + const BoundingBox<1> interval = create_unit_bounding_box<1>(); + + std::vector roots; + RootFinder root_finder; + root_finder.find_roots(functions, interval, roots); + + for (unsigned int i = 0; i < roots.size(); i++) + { + deallog << roots[i]; + if (i < roots.size() - 1) + deallog << ", "; + } + deallog << std::endl; +} + + + +// Test that the roots we get back from RootFinder are sorted and +// not duplicated. +// +// Call find_roots with 3 linear functions f_i(x) = x - x_i, +// where x_0 = 0.75, x_1 = 0.25, x_2 = 0.25 +// and check that RootFinder gives back the vector {0.25, 0.75}. +void +test_roots_sorted_not_duplicated() +{ + deallog << "test_roots_sorted_not_duplicated" << std::endl; + + std::vector> linear_functions; + + const std::vector roots = {.75, .25, .25}; + for (unsigned int i = 0; i < roots.size(); ++i) + { + Tensor<1, 1> normal; + normal[0] = 1; + const Point<1> point(roots.at(i)); + linear_functions.push_back(Functions::LevelSet::Plane<1>(point, normal)); + } + + const std::vector>> functions( + linear_functions.begin(), linear_functions.end()); + + find_and_print_roots(functions); +} + + + +/* + * The function: + * f(x) = C(x - x_0)^2 + y_0 + */ +class QuadraticFunction : public Function<1> +{ +public: + QuadraticFunction(const double C, const double x_0, const double y_0) + : C(C) + , x_0(x_0) + , y_0(y_0) + {} + + double + value(const Point<1> &point, const unsigned int component = 0) const override + { + return C * std::pow(point(0) - x_0, 2) + y_0; + }; + + Tensor<1, 1> + gradient(const Point<1> & point, + const unsigned int component = 0) const override + { + Tensor<1, 1> grad; + grad[0] = 2 * C * (point(0) - x_0); + + return grad; + }; + + SymmetricTensor<2, 1> + hessian(const Point<1> & point, + const unsigned int component = 0) const override + { + SymmetricTensor<2, 1> grad; + grad[0][0] = 2 * C; + + return grad; + }; + +private: + const double C; + const double x_0; + const double y_0; +}; + + + +// Test that RootFinder can find both roots of the function +// f(x) = 4(x-0.5)^2 - 0.25 +// which are x_0 = 0.25 and x_1 = 0.75. +void +test_find_both_roots() +{ + deallog << "test_find_both_roots" << std::endl; + + const QuadraticFunction function(4, 0.5, -0.25); + + std::vector>> functions; + functions.push_back(function); + + find_and_print_roots(functions); +} + + + +int +main() +{ + initlog(); + test_roots_sorted_not_duplicated(); + deallog << std::endl; + test_find_both_roots(); +} diff --git a/tests/non_matching/root_finder.output b/tests/non_matching/root_finder.output new file mode 100644 index 0000000000..e4202fe256 --- /dev/null +++ b/tests/non_matching/root_finder.output @@ -0,0 +1,6 @@ + +DEAL::test_roots_sorted_not_duplicated +DEAL::0.250000, 0.750000 +DEAL:: +DEAL::test_find_both_roots +DEAL::0.250000, 0.750000 diff --git a/tests/non_matching/tensor_point_with_1D_quadrature.cc b/tests/non_matching/tensor_point_with_1D_quadrature.cc new file mode 100644 index 0000000000..3f2d3a1c4d --- /dev/null +++ b/tests/non_matching/tensor_point_with_1D_quadrature.cc @@ -0,0 +1,79 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/** + * Test the function tensor_point_with_1D_quadrature() + * in NonMatching::internal::QuadratureGeneratorImplementation. + */ + +#include + +#include + +#include "../tests.h" + +#include "quadrature_printing.h" + +using namespace dealii; +using namespace NonMatching::internal::QuadratureGeneratorImplementation; + + +/* + * Set up a (dim-1)-dimensional point and a 1D-quadrature. Call + * tensor_point_with_1D_quadrature for each possible coordinate direction + * and print the resulting dim-dimensional quadrature. + */ +template +void +create_and_output_quadrature_for_each_direction() +{ + deallog << "dim=" << dim << std::endl; + + const unsigned int n_quadrature_points = 2; + const QGaussLobatto<1> quadrature1D(n_quadrature_points); + // Choose the points coordinates to something + // easily distinguished. + Point point; + for (int i = 0; i < dim - 1; ++i) + { + point(i) = 10 * (i + 1); + } + // Both points in the 1D-quadrature have weight 1/2 so + // this should also be the weight of the points in the final + // quadrature. + const double weight = 5; + const double start = -1, end = 1; + + for (int direction = 0; direction < dim; ++direction) + { + deallog << "direction=" << direction << std::endl; + ExtendableQuadrature result; + tensor_point_with_1D_quadrature( + point, weight, quadrature1D, start, end, direction, result); + print_quadrature(result); + deallog << std::endl; + } +} + + + +int +main() +{ + initlog(); + create_and_output_quadrature_for_each_direction<1>(); + create_and_output_quadrature_for_each_direction<2>(); + create_and_output_quadrature_for_each_direction<3>(); +} diff --git a/tests/non_matching/tensor_point_with_1D_quadrature.output b/tests/non_matching/tensor_point_with_1D_quadrature.output new file mode 100644 index 0000000000..72077f8672 --- /dev/null +++ b/tests/non_matching/tensor_point_with_1D_quadrature.output @@ -0,0 +1,28 @@ + +DEAL::dim=1 +DEAL::direction=0 +DEAL::-1.00000, 5.00000 +DEAL::1.00000, 5.00000 +DEAL:: +DEAL::dim=2 +DEAL::direction=0 +DEAL::-1.00000, 10.0000, 5.00000 +DEAL::1.00000, 10.0000, 5.00000 +DEAL:: +DEAL::direction=1 +DEAL::10.0000, -1.00000, 5.00000 +DEAL::10.0000, 1.00000, 5.00000 +DEAL:: +DEAL::dim=3 +DEAL::direction=0 +DEAL::-1.00000, 10.0000, 20.0000, 5.00000 +DEAL::1.00000, 10.0000, 20.0000, 5.00000 +DEAL:: +DEAL::direction=1 +DEAL::20.0000, -1.00000, 10.0000, 5.00000 +DEAL::20.0000, 1.00000, 10.0000, 5.00000 +DEAL:: +DEAL::direction=2 +DEAL::10.0000, 20.0000, -1.00000, 5.00000 +DEAL::10.0000, 20.0000, 1.00000, 5.00000 +DEAL:: diff --git a/tests/non_matching/up_through_dimension_creator.cc b/tests/non_matching/up_through_dimension_creator.cc new file mode 100644 index 0000000000..653be417bb --- /dev/null +++ b/tests/non_matching/up_through_dimension_creator.cc @@ -0,0 +1,166 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/** + * Test the class UpThroughDimensionCreator + * in NonMatching::internal::QuadratureGeneratorImplementation. + */ + +#include +#include + +#include + +#include + +#include "../tests.h" + +#include "quadrature_printing.h" + +using namespace dealii; +using namespace NonMatching::internal::QuadratureGeneratorImplementation; + + +// Print each quadrature in the incoming QPartitioning. +template +void +print(const QPartitioning &q_partitioning) +{ + deallog << "Negative" << std::endl; + print_quadrature(q_partitioning.negative); + deallog << "Positive" << std::endl; + print_quadrature(q_partitioning.positive); + deallog << "Indefinite" << std::endl; + print_quadrature(q_partitioning.indefinite); + deallog << "Surface" << std::endl; + print_surface_quadrature(q_partitioning.surface); +} + + + +// Let the height function direction be dim - 1 and let the lower dimensional +// quadrature contain a single point. Call UpThroughDimensionCreator with the +// incoming level set function over the unit box to generate a QPartitioning. +// Print the quadratures in the partitioning to make sure it is correct. +template +void +create_and_print_partitioning(const Function &level_set) +{ + const hp::QCollection<1> q_collection1D(QGauss<1>(2)); + const NonMatching::AdditionalQGeneratorData additional_data; + + UpThroughDimensionCreator up_through_dimension_creator( + q_collection1D, additional_data); + + std::vector>> level_sets; + level_sets.push_back(level_set); + const BoundingBox box = create_unit_bounding_box(); + const QMidpoint low_dim_quadrature; + const unsigned int height_function_direction = dim - 1; + + QPartitioning q_partitioning; + up_through_dimension_creator.generate(level_sets, + box, + low_dim_quadrature, + height_function_direction, + q_partitioning); + print(q_partitioning); +} + + + +// Set up a level set function with the zero contour along, x_{dim-1} = 0.5 +// Call create_and_print_partitioning to test that the points are added as +// expected: +// "negative" points should have x_{dim-1} \in (0, 0.5) +// "positive" points should have x_{dim-1} \in (0.5, 1) +// surface points should have x_{dim-1} = 0.5 +template +void +test_cut_through_center() +{ + deallog << "test_cut_through_center" << std::endl; + deallog << std::endl; + + Point point_through_plane = .5 * Point::unit_vector(dim - 1); + Tensor<1, dim> plane_normal = Point::unit_vector(dim - 1); + const Functions::LevelSet::Plane level_set(point_through_plane, + plane_normal); + + create_and_print_partitioning(level_set); +} + + + +// Fabricate the case when we have missed roots when creating the quadrature in +// the lower dimensions. See the comment in the implementation of +// UpThroughDimensionCreator::create_surface_point(..). +// +// In this test, the zero contour goes outside the cell but close to the "bottom +// face" at x_{dim-1} = 0. Check that the surface quadrature points gets placed +// on x_{dim-1} = 0. This is the "least bad" option we have. +template +void +test_missed_roots_on_bottom_face() +{ + deallog << "test_missed_roots_on_bottom_face" << std::endl; + deallog << std::endl; + + const Tensor<1, dim> plane_normal = Point::unit_vector(dim - 1); + Point point_in_plane; + point_in_plane[dim - 1] = -.1; + const Functions::LevelSet::Plane level_set(point_in_plane, plane_normal); + + create_and_print_partitioning(level_set); +} + + + +// Same test as above, but the zero contour just outside the "top face" at +// x_{dim-1} = 1. Check that the surface quadrature points are placed on +// x_{dim-1} = 1. +template +void +test_missed_roots_on_top_face() +{ + deallog << "test_missed_roots_on_top_face" << std::endl; + deallog << std::endl; + + const Tensor<1, dim> plane_normal = Point::unit_vector(dim - 1); + Point point_in_plane; + point_in_plane[dim - 1] = 1.1; + const Functions::LevelSet::Plane level_set(point_in_plane, plane_normal); + + create_and_print_partitioning(level_set); +} + + + +int +main() +{ + initlog(); + + const int dim = 2; + + test_cut_through_center(); + deallog << std::endl; + + test_missed_roots_on_bottom_face(); + deallog << std::endl; + + test_missed_roots_on_top_face(); + deallog << std::endl; +} diff --git a/tests/non_matching/up_through_dimension_creator.output b/tests/non_matching/up_through_dimension_creator.output new file mode 100644 index 0000000000..336cbb1fc3 --- /dev/null +++ b/tests/non_matching/up_through_dimension_creator.output @@ -0,0 +1,33 @@ + +DEAL::test_cut_through_center +DEAL:: +DEAL::Negative +DEAL::0.500000, 0.105662, 0.250000 +DEAL::0.500000, 0.394338, 0.250000 +DEAL::Positive +DEAL::0.500000, 0.605662, 0.250000 +DEAL::0.500000, 0.894338, 0.250000 +DEAL::Indefinite +DEAL::Surface +DEAL::0.500000, 0.500000, 1.00000, 0.00000, 1.00000 +DEAL:: +DEAL::test_missed_roots_on_bottom_face +DEAL:: +DEAL::Negative +DEAL::Positive +DEAL::0.500000, 0.211325, 0.500000 +DEAL::0.500000, 0.788675, 0.500000 +DEAL::Indefinite +DEAL::Surface +DEAL::0.500000, 0.00000, 1.00000, 0.00000, 1.00000 +DEAL:: +DEAL::test_missed_roots_on_top_face +DEAL:: +DEAL::Negative +DEAL::0.500000, 0.211325, 0.500000 +DEAL::0.500000, 0.788675, 0.500000 +DEAL::Positive +DEAL::Indefinite +DEAL::Surface +DEAL::0.500000, 1.00000, 1.00000, 0.00000, 1.00000 +DEAL::