You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

969 lines
41 KiB

#ifndef ALGOIM_QUADRATURE_MULTIPOLY_HPP
#define ALGOIM_QUADRATURE_MULTIPOLY_HPP
// High-order accurate quadrature algorithms for multi-component domains implicitly
// defined by (one or more) multivariate Bernstein polynomials, based on the algorithms
// developed in the paper
// R. I. Saye, High-order quadrature on multi-component domains implicitly defined
// by multivariate polynomials, Journal of Computational Physics, 448, 110720 (2022),
// https://doi.org/10.1016/j.jcp.2021.110720
//
// See examples/examples_quad_multipoly.cpp, as well as the short tutorial on the GitHub
// page https://algoim.github.io/ for examples of usage.
// The algorithms make use of various "masking" operations. A mask divides the reference
// cube [0,1]^N into a regular grid of M x ... x M subcells; on each subcell, a mask has
// binary value 0 or 1 and indicates whether its accompanying polynomial is provably
// nonzero on that subcell, or if its roots can be ignored. Typically, M = 4 or 8 is
// a good choice, and this is specified by the following macro def.
#include <cassert>
#include <cmath>
#include <iterator>
#include <vector>
#include "organizer/primitive.hpp"
#define ALGOIM_M 8
#include <algorithm>
#include "real.hpp"
#include <cmath>
#include "uvector.hpp"
#include "booluarray.hpp"
#include "multiloop.hpp"
#include "xarray.hpp"
#include "polyset.hpp"
#include "sparkstack.hpp"
#include "gaussquad.hpp"
#include "tanhsinh.hpp"
#include "bernstein.hpp"
#define STOP_WHEN_BLOCKED true
#if STOP_WHEN_BLOCKED
#include <chrono>
static thread_local auto timerStart = std::chrono::high_resolution_clock::now();
static thread_local bool stopCurrentQuadrature = false;
const double MAX_DURATION = 3; // seconds
#endif
// #define
namespace algoim
{
inline namespace v1
{
namespace detail
{
/**
* g 和 gmask之所以用ptr,是因为其不一定存在,方便传入空
*/
template <int N>
booluarray<N, ALGOIM_M> mask_driver(const xarray<real, N>& f,
const booluarray<N, ALGOIM_M>& fmask,
const xarray<real, N>* g,
const booluarray<N, ALGOIM_M>* gmask)
{
booluarray<N, ALGOIM_M> mask(false);
auto helper = [&](auto&& self, uvector<int, N> a, uvector<int, N> b) {
bool overlap = false;
for (MultiLoop<N> i(a, b); ~i; ++i)
if (fmask(i()) && (!gmask || (*gmask)(i()))) overlap = true;
if (!overlap) return;
real eps = 0.05 / ALGOIM_M;
uvector<real, N> xa, xb;
for (int dim = 0; dim < N; ++dim) {
xa(dim) = real(a(dim)) / ALGOIM_M - eps;
xb(dim) = real(b(dim)) / ALGOIM_M + eps;
}
if (g) {
xarray<real, N> fab(nullptr, f.ext()), gab(nullptr, g->ext());
algoim_spark_alloc(real, fab, gab);
bernstein::deCasteljau(f, xa, xb, fab);
bernstein::deCasteljau(*g, xa, xb, gab);
if (bernstein::orthantTest(fab, gab)) return;
} else {
xarray<real, N> fab(nullptr, f.ext());
algoim_spark_alloc(real, fab);
bernstein::deCasteljau(f, xa, xb, fab);
if (bernstein::uniformSign(fab) != 0) return;
}
if (b(0) - a(0) == 1) {
assert(all(b - a == 1));
assert(fmask(a) && (!gmask || (*gmask)(a)));
mask(a) = true;
return;
}
assert(all(b - a > 1) && all((b - a) % 2 == 0));
uvector<int, N> delta = (b - a) / 2;
for (MultiLoop<N> i(0, 2); ~i; ++i) self(self, a + i() * delta, a + (i() + 1) * delta);
};
helper(helper, 0, ALGOIM_M);
return mask;
}
// Using orthant tests, compute a mask for the possible subrectangles for which f and g share
// common zeros, i.e., intersecting zero level sets; if the mask is false somewhere, then it is
// guaranteed that f and/or g were originally masked off at the same place, or that they
// definitively do not share common zeros in that subrectangle; if the returned mask is true,
// then shared zeros may exist (and with high likelihood)
template <int N>
booluarray<N, ALGOIM_M> intersectionMask(const xarray<real, N>& f,
const booluarray<N, ALGOIM_M>& fmask,
const xarray<real, N>& g,
const booluarray<N, ALGOIM_M>& gmask)
{
return mask_driver(f, fmask, &g, &gmask);
}
// Using orthant tests, compute a mask for the possible subrectangles for which f has a zero
// level set; if the mask is false somewhere, then it is guaranteed that f was originally
// masked off at the same place, or that it definitively does not have any zeros in that
// subrectangle; if the returned mask is true, then shared zeros may exist (and with high likelihood)
template <int N>
booluarray<N, ALGOIM_M> nonzeroMask(const xarray<real, N>& f, const booluarray<N, ALGOIM_M>& fmask)
{
return mask_driver<N>(f, fmask, nullptr, nullptr);
}
// Collapse a mask along dimension k by bitwise-OR-ing along columns
template <int N>
booluarray<N - 1, ALGOIM_M> collapseMask(const booluarray<N, ALGOIM_M>& mask, int k)
{
booluarray<N - 1, ALGOIM_M> r(false);
for (MultiLoop<N> i(0, ALGOIM_M); ~i; ++i)
if (mask(i())) r(remove_component(i(), k)) = true;
return r;
}
// Test if a mask is empty, i.e., all entries are false
template <int N>
bool maskEmpty(const booluarray<N, ALGOIM_M>& mask)
{
return mask.none();
}
// Test if a point x \in [0,1]^N is in a true subrectangle of a mask; if x is exactly
// on the border between two subrectangles, the left subrectangle shall be used
template <int N>
bool pointWithinMask(const booluarray<N, ALGOIM_M>& mask, const uvector<real, N>& x)
{
using std::floor;
uvector<int, N> cell;
for (int dim = 0; dim < N; ++dim)
cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast<int>(floor(x(dim) * ALGOIM_M))));
return mask(cell);
}
// Test if a point {x + alpha e_k} is in a true subrectangle of a mask for some alpha \in [0,1]
template <int N>
bool lineIntersectsMask(const booluarray<N, ALGOIM_M>& mask, const uvector<real, N - 1>& x, int k)
{
using std::floor;
if constexpr (N > 1) {
uvector<int, N> cell;
for (int dim = 0; dim < N; ++dim)
if (dim < k)
cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast<int>(floor(x(dim) * ALGOIM_M))));
else if (dim > k)
cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast<int>(floor(x(dim - 1) * ALGOIM_M))));
for (int i = 0; i < ALGOIM_M; ++i) {
cell(k) = i;
if (mask(cell)) return true;
}
return false;
} else
return !maskEmpty(mask);
}
template <int N>
void restrictToFace(const xarray<real, N>& a, int k, int side, xarray<real, N - 1>& out)
{
assert(0 <= k && k < N && (side == 0 || side == 1));
assert(all(out.ext() == remove_component(a.ext(), k)));
int P = a.ext(k);
for (auto i = out.loop(); ~i; ++i) {
uvector<int, N> j;
for (int dim = 0; dim < N; ++dim) j(dim) = (dim < k) ? i(dim) : ((dim == k) ? side * (P - 1) : i(dim - 1));
out.l(i) = a.m(j);
}
}
template <int N>
booluarray<N - 1, ALGOIM_M> restrictToFace(const booluarray<N, ALGOIM_M>& a, int k, int side)
{
assert(0 <= k && k < N && (side == 0 || side == 1));
booluarray<N - 1, ALGOIM_M> r;
for (MultiLoop<N - 1> i(0, ALGOIM_M); ~i; ++i) {
uvector<int, N> j;
for (int dim = 0; dim < N; ++dim) j(dim) = (dim < k) ? i(dim) : ((dim == k) ? side * (ALGOIM_M - 1) : i(dim - 1));
r(i()) = a(j);
}
return r;
}
// Compute determinant of the given matrix using QR + Givens rotations + column
// pivoting, along with approximated rank
// in: square matrix A, which will be overwritten
template <typename T>
T det_qr(xarray<T, 2>& A, int& rank, T tol = 10.0)
{
assert(A.ext(0) == A.ext(1) && A.ext(0) > 0);
using std::abs;
using std::max;
T det = 1.0;
int n = A.ext(0);
T max_diag_r = 0.0;
for (int j = 0; j < n; ++j) {
T m = -1;
int k = -1;
for (int i = j; i < n; ++i) {
T mag = 0;
for (int a = 0; a < n; ++a) mag += util::sqr(A(a, i));
if (k == -1 || mag >= m) {
m = mag;
k = i;
}
}
assert(j <= k && k < n);
if (k != j) {
for (int a = 0; a < n; ++a) std::swap(A(a, j), A(a, k));
det *= -1.0;
}
for (int i = n - 1; i >= j + 1; --i) {
T c, s;
util::givens_get(A(i - 1, j), A(i, j), c, s);
for (int k = j; k < n; ++k) util::givens_rotate(A(i - 1, k), A(i, k), c, s);
}
det *= A(j, j);
max_diag_r = max(max_diag_r, abs(A(j, j)));
}
tol *= max_diag_r * n * std::numeric_limits<T>::epsilon();
rank = 0;
for (int i = 0; i < n; ++i)
if (abs(A(i, i)) > tol) ++rank;
return det;
}
// Determine the largest possible degree of the resultant of two general polynomials
template <int N>
uvector<int, N - 1> resultantExtent(const uvector<int, N>& p, const uvector<int, N>& q, int dim)
{
uvector<int, N - 1> ext;
for (int i = 0; i < N - 1; ++i) {
int ii = (i < dim) ? i : i + 1;
ext(i) = (p(dim) - 1) * (q(ii) - 1) + (q(dim) - 1) * (p(ii) - 1) + 1;
}
return ext;
}
// Determine the largest possible degree of the discriminant of a polynomial
template <int N>
uvector<int, N - 1> discriminantExtent(const uvector<int, N>& p, int dim)
{
uvector<int, N - 1> ext;
for (int i = 0; i < N - 1; ++i) {
int ii = (i < dim) ? i : i + 1;
ext(i) = (2 * p(dim) - 3) * (p(ii) - 1) + 1;
}
return ext;
}
// Compute the resultant of p and q and store the result in out
// ========================================= NOTE =========================================
// This is a heavily simplified method and does not handle rank deficiency caused by, e.g.,
// common polynomial factors, nor does it handle ill-conditioning caused by extreme values,
// among various other aspects. If your application requires this kind of special handling,
// consider contacting the author of this code for suggestions.
// ========================================= NOTE =========================================
template <int N>
bool resultant_core(const xarray<real, N>& p, const xarray<real, N>* q, int k, xarray<real, N - 1>& out)
{
assert(0 <= k && k < N);
int P = p.ext(k);
int Q = q ? q->ext(k) : P - 1;
int M = (P == Q) ? P - 1 : P + Q - 2;
xarray<real, 2> mat(nullptr, uvector<int, 2>{M, M});
/**gjj***/
// 当Q == 0,应该是对p求P类点,且p完全垂直降维方向的情形,这样的R类点可以直接丢弃
// 因为产生的R类点和Q类点是重合的(平行降维方向的多项式一定和上下边界交于同一位置,其降维后也是R类点降维后的结果)
// 当P和Q都<=1,要么是Q == 0,要么是两个垂直降维方向的多项式求R类点,此时交线退化为点,也丢弃
if (Q == 0 || (P <= 1 && Q <= 1)) { return false; }
if (!(P >= 2 && Q >= 1 && M >= 1)) {
int aaa = 1;
int bbb = 1;
}
assert(Q >= 1 && M >= 1);
/***gjj***/
// assert(P >= 2 && Q >= 1 && M >= 1);
xarray<real, N - 1> f(nullptr, out.ext());
real *pk, *qk;
algoim_spark_alloc(real, f, mat);
algoim_spark_alloc(real, &pk, P, &qk, Q);
for (auto i = f.loop(); ~i; ++i) {
uvector<real, N - 1> x;
for (int dim = 0; dim < N - 1; ++dim) x(dim) = bernstein::modifiedChebyshevNode(i(dim), f.ext(dim));
bernstein::collapseAlongAxis(p, x, k, pk);
if (q)
bernstein::collapseAlongAxis(*q, x, k, qk);
else
bernstein::bernsteinDerivative(pk, P, qk);
if (P == Q)
bernstein::bezoutMatrix(pk, qk, P, mat);
else
bernstein::sylvesterMatrix(pk, P, qk, Q, mat);
int rank;
f.l(i) = det_qr(mat, rank);
}
// Interpolate the resultant on the tensor-product grid
bernstein::normalise(f);
bernstein::bernsteinInterpolate(f, std::pow(100.0 * std::numeric_limits<real>::epsilon(), 1.0 / (N - 1)), out);
// Try for polynomial degree reduction
bool b = bernstein::autoReduction(out, 1e4 * std::numeric_limits<real>::epsilon());
// If able to reduce the degree, recompute the resultant on the lower-degree poly space
// which is expected to have better conditioning
if (b) resultant_core(p, q, k, out);
return true;
}
// Compute the pseudo-resultant R(p,q) along dimension k
template <int N>
bool resultant(const xarray<real, N>& p, const xarray<real, N>& q, int k, xarray<real, N - 1>& out)
{
return resultant_core(p, &q, k, out);
}
// Compute the (intentially unnormalised) pseudo-discriminant R(p,p') along dimension k
template <int N>
bool discriminant(const xarray<real, N>& p, int k, xarray<real, N - 1>& out)
{
/**gjj ***/
if (p.ext(k) <= 1) {
return false; // p和降维方向完全平行。这样的p可以不要,因为一模一样的N-1维psi会由Q类点等价刻画 }
} /**gjj ***/
xarray<real, N> prime(nullptr, inc_component(p.ext(), k, -1));
algoim_spark_alloc(real, prime);
bernstein::bernsteinDerivative(p, k, prime);
return resultant_core(p, &prime, k, out);
}
// Using the polynomials in phi, eliminate the axis k by restricting to faces, computing
// discriminants and resultants, and storing the computed polynomials in psi
template <int N>
void eliminate_axis(PolySet<N, ALGOIM_M>& phi, int k, PolySet<N - 1, ALGOIM_M>& psi)
{
static_assert(N >= 2, "N >= 2 required to eliminate axis");
assert(0 <= k && k < N);
assert(psi.count() == 0);
// For every phi(i) ...
for (int i = 0; i < phi.count(); ++i) {
#if STOP_WHEN_BLOCKED
auto timerEnd = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> duration = timerEnd - timerStart;
// auto duration = std::chrono::duration_cast<std::chrono::seconds>(timerEnd - timerStart);
if (duration.count() > MAX_DURATION) {
stopCurrentQuadrature = true;
return;
}
#endif
const auto& p = phi.poly(i);
const auto& mask = phi.mask(i);
// Examine bottom and top faces in the k'th dimension
// Q
for (int side = 0; side <= 1; ++side) {
xarray<real, N - 1> p_face(nullptr, remove_component(p.ext(), k));
algoim_spark_alloc(real, p_face);
restrictToFace(p, k, side, p_face);
auto p_face_mask = nonzeroMask(p_face, restrictToFace(mask, k, side));
if (!maskEmpty(p_face_mask)) {
bernstein::autoReduction(p_face);
bernstein::normalise(p_face);
psi.push_back(p_face, p_face_mask);
}
}
// Consider discriminant
// P
xarray<real, N> p_k(nullptr, p.ext());
algoim_spark_alloc(real, p_k);
bernstein::elevatedDerivative(p, k, p_k);
auto disc_mask = intersectionMask(p, mask, p_k, mask);
if (!maskEmpty(disc_mask)) {
// note: computed disc might have lower degree than the following
uvector<int, N - 1> R = discriminantExtent(p.ext(), k);
xarray<real, N - 1> disc(nullptr, R);
algoim_spark_alloc(real, disc);
if (discriminant(p, k, disc)) {
bernstein::normalise(disc);
psi.push_back(disc, collapseMask(disc_mask, k));
}
}
}
// Consider every pairwise combination of resultants ...
// R
for (int i = 0; i < phi.count(); ++i)
for (int j = i + 1; j < phi.count(); ++j) {
#if STOP_WHEN_BLOCKED
auto timerEnd = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> duration = timerEnd - timerStart;
// auto duration = std::chrono::duration_cast<std::chrono::seconds>(timerEnd - timerStart);
if (duration.count() > MAX_DURATION) {
stopCurrentQuadrature = true;
return;
}
#endif
const auto& p = phi.poly(i);
const auto& pmask = phi.mask(i);
const auto& q = phi.poly(j);
const auto& qmask = phi.mask(j);
auto mask = intersectionMask(p, pmask, q, qmask);
if (!maskEmpty(mask)) {
// note: computed resultant might have lower degree than the following
uvector<int, N - 1> R = resultantExtent(p.ext(), q.ext(), k);
xarray<real, N - 1> res(nullptr, R);
algoim_spark_alloc(real, res);
if (resultant(p, q, k, res)) {
bernstein::normalise(res);
psi.push_back(res, collapseMask(mask, k));
}
}
};
}
/**
@param k 降维方向
*/
template <int N>
void refWNVCompress(const PolySet<N, ALGOIM_M>& originPhis,
const std::vector<std::vector<int>>& refWNV,
std::vector<std::vector<int>>& refWNVCompressed,
int k)
{
assert(refWNV.size() >= 2);
assert(refWNV[0].size() == originPhis.count());
refWNVCompressed = std::vector<std::vector<int>>(refWNV.size() / 2, std::vector<int>(refWNV[0].size(), 0));
int base = pow(2, k);
for (int i = 0; i < refWNVCompressed.size(); ++i) {
uvector<int, N - 1> xBaseVec;
for (int j = 0; j < N - 1; ++j) { xBaseVec(j) = (i & (1 << j)) >> j; }
int remainder = i % base;
int x0 = ((i - remainder) << 1) + remainder;
int x1 = ((i - remainder + 1) << 1) + remainder;
for (int j = 0; j < refWNVCompressed[0].size(); ++j) {
if (refWNV[x0][j] == 2 || refWNV[x1][j] == 2) {
refWNVCompressed[i][j] = 2; // both
continue;
}
if (refWNV[x0][j] != refWNV[x1][j]) {
// 一个是1一个是0
refWNVCompressed[i][j] = 2; // both
continue;
}
if (!detail::lineIntersectsMask(originPhis.mask(i), xBaseVec, k)) {
// 没交
if (refWNV[x0][j] == 1 && refWNV[x1][j] == 1) {
refWNVCompressed[i][j] = 1; // in
} else if (refWNV[x0][j] == 0 && refWNV[x1][j] == 0) {
refWNVCompressed[i][j] = 0; // out
}
} else {
// 交了
refWNVCompressed[i][j] = 2; // both
}
}
}
}
// Compute the 'score' across all dimensions
template <int N>
uvector<real, N> score_estimate(PolySet<N, ALGOIM_M>& phi, uvector<bool, N>& has_disc)
{
static_assert(N > 1, "score_estimate of practical use only with N > 1");
using std::abs;
uvector<real, N> s = 0;
has_disc = false;
// For every phi(i) ...
for (int i = 0; i < phi.count(); ++i) {
const auto& p = phi.poly(i);
const auto& mask = phi.mask(i);
// Accumulate onto score by sampling at midpoint of every subcell of mask
for (MultiLoop<N> j(0, ALGOIM_M); ~j; ++j)
if (mask(j())) {
uvector<real, N> x = (j() + 0.5) / real(ALGOIM_M);
uvector<real, N> g = bernstein::evalBernsteinPolyGradient(p, x);
real sum = 0;
for (int dim = 0; dim < N; ++dim) {
g(dim) = abs(g(dim));
sum += g(dim);
}
if (sum > 0) s += g / sum;
}
// Consider discriminant
// 这里是再用结式(intersection)来求判别式,即,求导数对应的多项式与原多项式的交集
xarray<real, N> p_k(nullptr, p.ext());
algoim_spark_alloc(real, p_k);
for (int k = 0; k < N; ++k) {
bernstein::elevatedDerivative(p, k, p_k);
auto disc_mask = intersectionMask(p, mask, p_k, mask);
has_disc(k) |= !maskEmpty(disc_mask);
}
}
return s;
}
} // namespace detail
// Main engine for generating high-order accurature quadrature schemes on multi-component domains
// implicitly defined by multivariate Bernstein polynomials in the unit hyperrectangle [0,1]^N.
// See examples_quad_multipoly.cpp for examples demonstrating the usage of these methods.
// See also additional examples provided on the GitHub documentation page,
// https://algoim.github.io/
// After building the quadrature hierarchy (via dimension reduction), the specific kind
// of quadrature scheme applied to the one-dimensional line integrals is chosen by the
// user via a QuadStrategy parameter:
// AlwaysGL: applies Gauss-Legendre quadrature on every interval, regardless of
// the level; this strategy is generally good for small-to-medium q
// AlwaysTS: applies tanh-sinh quadrature on every interval, regardless of the level;
// this strategy is generally the worst of all three and is provided mainly
// for exploratory purposes
// AutoMixed: apply the automated strategy discussed in detail in the paper
// https://doi.org/10.1016/j.jcp.2021.110720; briefly, the method: (i) applies
// tanh-sinh on base integrals for which vertical tangents are detected in
// the height-function-based representation of their implicitly-defined
// geometry; (ii) applies Gauss-Legendre quadrature on the inner-most
// integral, and on all base integrals whose geometry is the graph of a
// multi-valued height function devoid of vertical tangents/branching.
enum QuadStrategy { AlwaysGL, AlwaysTS, AutoMixed };
template <int N>
struct ImplicitPolyQuadrature {
enum IntegralType { Inner, OuterSingle, OuterAggregate };
PolySet<N, ALGOIM_M> phi; // Given N-dimensional polynomials
int k; // Elimination axis/height direction; k = N if there are no interfaces
// 这种struct的递归定义一般是不允许的,因为会内存爆炸,但是模板类可以通过特化定义递归的终止条件
// 所以有template<> struct ImplicitPolyQuadrature<0> {};
ImplicitPolyQuadrature<N - 1> base; // Base polynomials corresponding to removal of axis k
bool auto_apply_TS; // If quad method is auto chosen, indicates whether TS is applied // what is ts ? tanh-sinh?
IntegralType type; // Whether an inner integral, or outer of two kinds
std::array<std::tuple<int, ImplicitPolyQuadrature<N - 1>>, N - 1>
base_other; // Stores other base cases, besides k, when in aggregate mode
// Default ctor sets to an uninitialised state
ImplicitPolyQuadrature() : k(-1) {}
// Build quadrature hierarchy for a domain implicitly defined by a single polynomial
ImplicitPolyQuadrature(const xarray<real, N>& p)
{
#if STOP_WHEN_BLOCKED
timerStart = std::chrono::high_resolution_clock::now();
#endif
auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true));
if (!detail::maskEmpty(mask)) {
phi.push_back(p, mask);
int a = 1;
}
build(true, false);
}
// Build quadrature hierarchy for a domain implicitly defined by two polynomials
ImplicitPolyQuadrature(const xarray<real, N>& p, const xarray<real, N>& q)
{
#if STOP_WHEN_BLOCKED
timerStart = std::chrono::high_resolution_clock::now();
#endif
{
// auto tmp = booluarray<N, ALGOIM_M>(false);
auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true));
// tmp(0) = true;
// auto res = !detail::maskEmpty(tmp); // 似乎是测试代码?
if (!detail::maskEmpty(mask)) phi.push_back(p, mask);
}
{
auto mask = detail::nonzeroMask(q, booluarray<N, ALGOIM_M>(true));
if (!detail::maskEmpty(mask)) phi.push_back(q, mask);
}
build(true, false);
}
// Build quadrature hierarchy for a domain implicitly defined by multiple polynomials
ImplicitPolyQuadrature(const std::vector<xarray<real, N>>& ps)
{
#if STOP_WHEN_BLOCKED
timerStart = std::chrono::high_resolution_clock::now();
#endif
for (const auto& p : ps) {
auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true));
if (!detail::maskEmpty(mask)) phi.push_back(p, mask);
}
build(true, false);
}
// ImplicitPolyQuadrature(const std::vector<organizer::CompleteTensorRep>& completeTensorReps)
// {
// for (const auto& ctr : completeTensorReps) {
// auto mask = detail::nonzeroMask(ctr.compositedBernstein, booluarray<N, ALGOIM_M>(true));
// if (!detail::maskEmpty(mask)) phi.push_back(ctr.compositedBernstein, mask);
// }
// build(true, false);
// }
// ImplicitPolyQuadrature(const std::vector<organizer::CompleteTensorRep>& completeTensorReps,
// const std::vector<int>& polyIndices)
// {
// for (auto i : polyIndices) {
// auto mask = detail::nonzeroMask(completeTensorReps[i].compositedBernstein, booluarray<N, ALGOIM_M>(true));
// if (!detail::maskEmpty(mask)) phi.push_back(completeTensorReps[i].compositedBernstein, mask);
// }
// build(true, false);
// }
ImplicitPolyQuadrature(const std::vector<tensor3>& tensors, const std::vector<int>& polyIndices)
{
#if STOP_WHEN_BLOCKED
timerStart = std::chrono::high_resolution_clock::now();
#endif
for (auto i : polyIndices) {
auto mask = detail::nonzeroMask(tensors[i], booluarray<N, ALGOIM_M>(true));
if (!detail::maskEmpty(mask)) phi.push_back(tensors[i], mask);
}
build(true, false);
}
ImplicitPolyQuadrature(const std::vector<organizer::MinimalPrimitiveRep>& minimalReps, const std::vector<int>& polyIndices)
{
#if STOP_WHEN_BLOCKED
timerStart = std::chrono::high_resolution_clock::now();
#endif
for (auto i : polyIndices) {
auto mask = detail::nonzeroMask(minimalReps[i].tensor, booluarray<N, ALGOIM_M>(true));
if (!detail::maskEmpty(mask)) phi.push_back(minimalReps[i].tensor, mask);
}
build(true, false);
}
// Build quadrature hierarchy for a given domain implicitly defined by two polynomials with user-defined masks
ImplicitPolyQuadrature(const xarray<real, N>& p,
const booluarray<N, ALGOIM_M>& pmask,
const xarray<real, N>& q,
const booluarray<N, ALGOIM_M>& qmask)
{
#if STOP_WHEN_BLOCKED
timerStart = std::chrono::high_resolution_clock::now();
#endif
{
auto mask = detail::nonzeroMask(p, pmask);
if (!maskEmpty(mask)) phi.push_back(p, mask);
}
{
auto mask = detail::nonzeroMask(q, qmask);
if (!maskEmpty(mask)) phi.push_back(q, mask);
}
build(true, false);
}
// ssss
// Assuming phi has been instantiated, determine elimination axis and build base
void build(bool outer, bool auto_apply_TS)
{
#if STOP_WHEN_BLOCKED
if (stopCurrentQuadrature) { return; }
#endif
type = outer ? OuterSingle : Inner;
this->auto_apply_TS = auto_apply_TS;
// If phi is empty, apply a tensor-product Gaussian quadrature
if (phi.count() == 0) {
k = N;
this->auto_apply_TS = false;
return;
}
if constexpr (N == 1) {
// If in one dimension, there is only one choice of height direction and
// the recursive process halts
k = 0;
return;
} else {
// Compute score; penalise any directions which likely contain vertical tangents
uvector<bool, N> has_disc;
uvector<real, N> score = detail::score_estimate(phi, has_disc);
/**gjj */
if (max(abs(score)) == 0) score(0) = 1.0;
/**gjj */
assert(max(abs(score)) > 0);
score /= 2 * max(abs(score));
for (int i = 0; i < N; ++i)
if (!has_disc(i)) score(i) += 1.0; // 一个与k相切的多项式都没有时,才+1
// Choose height direction and form base polynomials; if tanh-sinh is being used at this
// level, suggest the same all the way down; moreover, suggest tanh-sinh if a non-empty
// discriminant mask has been found
k = argmax(score);
detail::eliminate_axis(phi, k, base.phi);
base.build(false, this->auto_apply_TS || has_disc(k));
// If this is the outer integral, and surface quadrature schemes are required, apply
// the dimension-aggregated scheme when necessary
// if (outer && has_disc(k)) {
// type = OuterAggregate;
// for (int i = 0; i < N; ++i)
// if (i != k) {
// auto& [kother, base] = base_other[i < k ? i : i - 1];
// kother = i;
// detail::eliminate_axis(phi, kother, base.phi);
// // In aggregate mode, triggered by non-empty discriminant mask,
// // base integrals always have T-S suggested
// base.build(false, this->auto_apply_TS || true);
// }
// }
}
}
// sssss
// Integrate a functional via quadrature of the base integral, adding the dimension k
template <typename F>
void integrate(QuadStrategy strategy, int q, const F& f)
{
assert(0 <= k && k <= N);
// If there are no interfaces, apply rectangular tensor-product Gauss-Legendre quadrature
// 所以这个是积分域包裹了整个box时的做法?
if (k == N) {
assert(!auto_apply_TS);
for (MultiLoop<N> i(0, q); ~i; ++i) {
uvector<real, N> x;
real w = 1.0;
for (int dim = 0; dim < N; ++dim) {
x(dim) = GaussQuad::x(q, i(dim)); // 0 <= i(dim) < q
w *= GaussQuad::w(q, i(dim));
}
f(x, w);
}
return;
}
// Determine maximum possible number of roots; used to allocate a small buffer
// ???
int max_count = 2;
for (size_t i = 0; i < phi.count(); ++i) max_count += phi.poly(i).ext(k) - 1;
// Base integral invokes the following integrand
auto integrand = [&](const uvector<real, N - 1>& xbase, real w) {
#if STOP_WHEN_BLOCKED
auto timerEnd = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> duration = timerEnd - timerStart;
// auto duration = std::chrono::duration_cast<std::chrono::seconds>(timerEnd - timerStart);
if (duration.count() > MAX_DURATION) {
stopCurrentQuadrature = true;
return;
}
#endif
// Allocate node buffer of sufficient size and initialise with {0, 1}
real* nodes;
algoim_spark_alloc(real, &nodes, max_count);
nodes[0] = 0.0;
nodes[1] = 1.0;
int count = 2;
// For every phi(i) ...
for (size_t i = 0; i < phi.count(); ++i) {
#if STOP_WHEN_BLOCKED
auto timerEnd = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> duration = timerEnd - timerStart;
// auto duration = std::chrono::duration_cast<std::chrono::seconds>(timerEnd - timerStart);
if (duration.count() > MAX_DURATION) {
stopCurrentQuadrature = true;
return;
}
#endif
const auto& p = phi.poly(i);
const auto& mask = phi.mask(i);
int P = p.ext(k);
// Ignore phi if its mask is void everywhere above the base point
if (!detail::lineIntersectsMask(mask, xbase, k)) continue;
// Restrict polynomial to axis-aligned line and compute its roots ????
real *pline, *roots;
algoim_spark_alloc(real, &pline, P, &roots, P - 1);
bernstein::collapseAlongAxis(p, xbase, k, pline);
int rcount = bernstein::bernsteinUnitIntervalRealRoots(pline, P, roots);
// Add all real roots in [0,1] which are also within masked region of phi
// 为什么[0,1]的real roots还必须要在mask为true的subregion?
for (int j = 0; j < rcount; ++j) {
auto x = add_component(xbase, k, roots[j]);
if (std::isnan(roots[j])) {
int aaa = 1;
int bbb = 1;
}
if (detail::pointWithinMask(mask, x)) nodes[count++] = roots[j];
}
};
/**gjj */
auto newEnd = std::remove_if(nodes, nodes + count, [](real x) { return std::isnan(x); });
count = std::distance(nodes, newEnd);
// for (int i = 0; i < count; ++i) {
// if (std::isnan(nodes[i])) { nanCount++; }
// }
// count -= nanCount;
/**gjj */
// Sort the nodes
std::sort(nodes, nodes + count);
if (!(nodes[0] == real(0) && nodes[count - 1] == real(1))) {
int aaa = 1;
int bbb = 1;
}
assert(nodes[0] == real(0) && nodes[count - 1] == real(1));
// Force nearly-degenerate sub-intervals to be exactly degenerate
// real tol = 10.0 * std::numeric_limits<real>::epsilon();
real tol = .01 * std::numeric_limits<real>::epsilon();
using std::abs;
for (int i = 1; i < count - 1; ++i)
if (abs(nodes[i]) < tol)
nodes[i] = 0.0;
else if (abs(nodes[i] - 1) < tol)
nodes[i] = 1.0;
else if (abs(nodes[i] - nodes[i + 1]) < tol)
nodes[i + 1] = nodes[i];
assert(nodes[0] == real(0) && nodes[count - 1] == real(1));
// Apply quadrature to non-degenerate sub-intervals
for (int i = 0; i < count - 1; ++i) {
real x0 = nodes[i];
real x1 = nodes[i + 1];
if (x0 == x1) continue;
// Choose between Gauss-Legendre and tanh-sinh according to the user-defined strategy
bool GL = true;
if (strategy == AlwaysTS) GL = false;
if (strategy == AutoMixed) GL = !this->auto_apply_TS;
if (GL)
for (int j = 0; j < q; ++j)
f(add_component(xbase, k, x0 + (x1 - x0) * GaussQuad::x(q, j)), w * (x1 - x0) * GaussQuad::w(q, j));
else
for (int j = 0; j < q; ++j)
f(add_component(xbase, k, TanhSinhQuadrature::x(q, j, x0, x1)),
w * TanhSinhQuadrature::w(q, j, x0, x1));
}
};
// When N = 1, the base case of recursion is invoked on a zero-dimensional point with unit weight
if constexpr (N > 1)
base.integrate(strategy, q, integrand);
else
integrand(uvector<real, 0>(), real(1));
}
// Surface-integrate a functional via quadrature of the base integral, adding the dimension k
template <typename F>
void integrate_surf(QuadStrategy strategy, int q, const F& f)
{
static_assert(N > 1, "surface integral only implemented in N > 1 dimensions");
assert(type == OuterSingle || type == OuterAggregate);
// If there is no interface, there is no surface integral
if (k == N) return;
// Base integral invokes the following integrand which operates in the height direction k_active
int k_active = -1;
auto integrand = [&](const uvector<real, N - 1>& xbase, real w) {
assert(0 <= k_active && k_active < N);
// For every phi(i) ...
for (size_t i = 0; i < phi.count(); ++i) {
const auto& p = phi.poly(i);
const auto& mask = phi.mask(i);
int P = p.ext(k_active);
// Ignore phi if its mask is void everywhere above the base point
if (!detail::lineIntersectsMask(mask, xbase, k_active)) continue;
// Compute roots of { x \mapsto phi(xbase + x e_k) }
real *pline, *roots;
algoim_spark_alloc(real, &pline, P, &roots, P - 1);
bernstein::collapseAlongAxis(p, xbase, k_active, pline);
int rcount = bernstein::bernsteinUnitIntervalRealRoots(pline, P, roots);
// Consider all real roots in (0,1) which are also within masked region of phi; evaluate
// integrand at interfacial points, multiplying weights by the effective surface Jacobian
for (int j = 0; j < rcount; ++j) {
auto x = add_component(xbase, k_active, roots[j]);
if (detail::pointWithinMask(mask, x)) {
using std::abs;
uvector<real, N> g = bernstein::evalBernsteinPolyGradient(p, x);
if (type == OuterAggregate) {
// When in aggregate mode, the scalar surf integral multiplies f by |n_k|^2, whose net effect
// is multiply weight by |n_k|; the flux surf integral multiplies f by sign(n_k) = sign(g(k))
real alpha = max(abs(g));
if (alpha > 0) {
g /= alpha;
alpha = abs(g(k_active)) / norm(g);
}
// Simplistic method to compute sign(n_k). NOTE: This method relies on a reasonable consistency
// between the gradient calculation of original polynomial, and that of the roots computed from
// pline; when near high-multiplicity roots, this simple method can break down; other, more
// sophisticated methods should be used in such cases, but these are not implemented here
f(x, w * alpha, set_component<real, N>(real(0.0), k_active, w * util::sign(g(k_active))));
} else {
// When in non-aggregate mode, the scalar surf integral multiples f by 1, whose net effect
// is multiply weight by 1/|n_k|; the flux surf integral multiplies f by n
uvector<real, N> n = g;
if (norm(n) > 0) n *= real(1.0) / norm(n);
real alpha = w * norm(g) / abs(g(k_active));
f(x, alpha, alpha * n);
}
}
}
};
};
// Apply primary base integral
k_active = k;
base.integrate(strategy, q, integrand);
// situation the
// If in aggregate mode, apply to other dimensions as well
if (type == OuterAggregate) {
for (int i = 0; i < N - 1; ++i) {
auto& [k, base] = base_other[i];
k_active = k;
base.integrate(strategy, q, integrand);
}
}
}
};
template <>
struct ImplicitPolyQuadrature<0> {
};
} // namespace v1
} // namespace algoim
#endif