GJJ
/
HeteQuadrature
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
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.
49 Commits
1 Branch
0 Tags
2.4 MiB
 
 
Tree: 4914ecdfc1
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '4914ecdfc1'
${ noResults }
HeteQuadrature/gjj/myDebug.hpp

1228 lines
49 KiB
Raw Blame History

#include <array>
#include <bitset>
#include <iostream>
#include <booluarray.hpp>
#include <cstddef>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include "bernstein.hpp"
#include "multiloop.hpp"
#include "quadrature_multipoly.hpp"
#include "binomial.hpp"
#include "real.hpp"
#include "sparkstack.hpp"
#include "uvector.hpp"
#include "vector"
#include "xarray.hpp"
#include <chrono>
using namespace algoim;
// Driver method which takes a functor phi defining a single polynomial in the reference
// rectangle [xmin, xmax]^N, of Bernstein degree P, along with an integrand function,
// and performances a q-refinement convergence study, comparing the computed integral
// with the given exact answers, for 1 <= q <= qMax.
template <int N, typename Phi, typename F>
void qConv1(const Phi& phi,
real xmin,
real xmax,
uvector<int, N> P,
const F& integrand,
int qMax,
real volume_exact,
real surf_exact)
{
// Construct Bernstein polynomial by mapping [0,1] onto bounding box [xmin,xmax]
xarray<real, N> phipoly(nullptr, P);
algoim_spark_alloc(real, phipoly);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return phi(xmin + x * (xmax - xmin)); }, phipoly);
// Build quadrature hierarchy
ImplicitPolyQuadrature<N> ipquad(phipoly);
// Functional to evaluate volume and surface integrals of given integrand
real volume, surf;
auto compute = [&](int q) {
volume = 0.0;
surf = 0.0;
// compute volume integral over {phi < 0} using AutoMixed strategy
ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
if (bernstein::evalBernsteinPoly(phipoly, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
});
// compute surface integral over {phi == 0} using AutoMixed strategy
ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, N>& x, real w, const uvector<real, N>& wn) {
surf += w * integrand(xmin + x * (xmax - xmin));
});
// scale appropriately
volume *= pow(xmax - xmin, N);
surf *= pow(xmax - xmin, N - 1);
};
// Compute results for all q and output in a convergence table
for (int q = 1; q <= qMax; ++q) {
compute(q);
std::cout << q << ' ' << volume << ' ' << surf << ' ' << std::abs(volume - volume_exact) / volume_exact << ' '
<< std::abs(surf - surf_exact) / surf_exact << std::endl;
}
}
enum BoolOp { Union, Intersection, Difference };
// Driver method which takes two phi functors defining two polynomials in the reference
// rectangle [xmin, xmax]^N, each of of Bernstein degree P, builds a quadrature scheme with the
// given q, and outputs it for visualisation in a set of VTP XML files
template <int N, typename F1, typename F2, typename F>
void qConvMultiPloy(const F1& fphi1,
const F2& fphi2,
real xmin,
real xmax,
const uvector<int, N>& P,
const F& integrand,
int q,
BoolOp op)
{
// Construct phi by mapping [0,1] onto bounding box [xmin,xmax]
xarray<real, N> phi1(nullptr, P), phi2(nullptr, P);
algoim_spark_alloc(real, phi1, phi2);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi1(xmin + x * (xmax - xmin)); }, phi1);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2);
// Build quadrature hierarchy
auto t = std::chrono::system_clock::now();
ImplicitPolyQuadrature<N> ipquad(phi1, phi2);
// ImplicitPolyQuadrature<N> ipquad(phi1);
// Functional to evaluate volume and surface integrals of given integrand
real volume, surf;
auto compute = [&](int q) {
volume = 0.0;
surf = 0.0;
// compute volume integral over {phi < 0} using AutoMixed strategy
ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
if (op == Union) {
if (bernstein::evalBernsteinPoly(phi1, x) < 0 || bernstein::evalBernsteinPoly(phi2, x) < 0)
volume += w * integrand(xmin + x * (xmax - xmin));
} else if (op == Intersection) {
if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0)
volume += w * integrand(xmin + x * (xmax - xmin));
} else if (op == Difference) {
if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) > 0)
volume += w * integrand(xmin + x * (xmax - xmin));
}
// if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
});
// compute surface integral over {phi == 0} using AutoMixed strategy
ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, N>& x, real w, const uvector<real, N>& wn) {
surf += w * integrand(xmin + x * (xmax - xmin));
});
// scale appropriately
volume *= pow(xmax - xmin, N);
surf *= pow(xmax - xmin, N - 1);
};
compute(q);
auto tAfter = std::chrono::system_clock::now();
std::chrono::duration<double> elapsed_seconds = tAfter - t;
std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
std::cout << "q volume: " << volume << std::endl;
}
template <int N, typename F>
void qConvMultiPloy2(real xmin, real xmax, const std::vector<xarray<real, N>>& phis, const F& integrand, int q)
{
auto t = std::chrono::system_clock::now();
ImplicitPolyQuadrature<N> ipquad(phis);
// Functional to evaluate volume and surface integrals of given integrand
real volume, surf;
auto compute = [&](int q) {
volume = 0.0;
surf = 0.0;
// compute volume integral over {phi < 0} using AutoMixed strategy
ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
// if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0) volume += w *
// integrand(xmin + x * (xmax - xmin)); if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin
// + x * (xmax - xmin));
bool in = true;
for (int i = 0; i < phis.size(); ++i) {
if (bernstein::evalBernsteinPoly(phis[i], x) >= 0) {
in = false;
break;
}
}
if (in) volume += w * integrand(xmin + x * (xmax - xmin));
});
// compute surface integral over {phi == 0} using AutoMixed strategy
// ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, N>& x, real w, const uvector<real, N>& wn) {
// surf += w * integrand(xmin + x * (xmax - xmin));
// });
// scale appropriately
volume *= pow(xmax - xmin, N);
surf *= pow(xmax - xmin, N - 1);
};
compute(q);
auto tAfter = std::chrono::system_clock::now();
std::chrono::duration<double> elapsed_seconds = tAfter - t;
std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
std::cout << "volume: " << volume << std::endl;
}
void fromCommon2Bernstein() {}
void testMultiPolys()
{
const int N = 3;
const uvector<int, N> P = 3;
real xmin = -1;
real xmax = 1;
std::vector<xarray<real, N>> phis; // 大量sphere
// xarray<real, N>* phis = new xarray<real, N>[100];
for (int i = 0; i < 100; i++) {
real centerX = 1.0 - rand() % 20 / 10.0, centerY = 1.0 - rand() % 20 / 10.0, centerZ = 1.0 - rand() % 20 / 10.0;
real r = 0.01 + rand() % 10 / 100.0;
auto fphi = [&](const uvector<real, 3>& xx) {
real x = xx(0) - centerX;
real y = xx(1) - centerY;
real z = xx(2) - centerZ;
return x * x + y * y + z * z - r * r;
};
xarray<real, 3> phi(nullptr, P);
algoim_spark_alloc(real, phi);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi(xmin + x * (xmax - xmin)); }, phi);
phis.emplace_back(phi);
}
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
// qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
std::cout << "testMultiPolys" << std::endl;
qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}
template <int N>
void power2BernsteinTensorProduct(const xarray<real, N>& phiPower, xarray<real, N>& phiBernstein)
{
xarrayInit(phiBernstein);
for (auto i = phiPower.loop(); ~i; ++i) {
// phi.l(i) = powerFactors.l(i);
real factorBase = phiPower.l(i);
if (factorBase == 0) continue;
auto traverseRange = phiPower.ext() - i();
std::vector<std::vector<real>> decompFactors(N, std::vector<real>(max(traverseRange), 0.));
for (int dim = 0; dim < N; ++dim) {
// Sigma
size_t nDim = phiPower.ext()(dim) - 1;
const real* binomNDim = Binomial::row(nDim);
for (int j = i(dim); j <= nDim; ++j) {
const real* binomJ = Binomial::row(j);
decompFactors[dim][j - i(dim)] = binomJ[i(dim)] / binomNDim[i(dim)];
}
}
xarray<real, N> subgrid(nullptr, traverseRange);
algoim_spark_alloc(real, subgrid);
for (auto ii = subgrid.loop(); ~ii; ++ii) {
real factor = factorBase;
for (int dim = 0; dim < N; ++dim) { factor *= decompFactors[dim][ii(dim)]; }
phiBernstein.m(i() + ii()) += factor;
}
}
}
template <int N>
real powerEvaluation(const xarray<real, N>& phi, const uvector<real, N>& x)
{
real res = 0;
for (auto i = phi.loop(); ~i; ++i) {
real item = phi.l(i);
auto exps = i();
for (int dim = 0; dim < N; ++dim) { item *= pow(x(dim), i(dim)); }
res += item;
}
return res;
}
template <int N, typename F>
void qConvBernstein(const xarray<real, N>& phi,
real xmin,
real xmax,
const F& integrand,
int q,
const std::vector<std::array<real, 4>>& halfFaces)
{
auto t = std::chrono::system_clock::now();
uvector<real, 3> testX(0., 0.75, 0.);
real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
auto vec1 = xarray2StdVector(phi);
std::cout << "eval bernstein without interpolation:" << testEvalBernstein << std::endl;
// Build quadrature hierarchy
ImplicitPolyQuadrature<N> ipquad(phi);
// ImplicitPolyQuadrature<N> ipquad(phi1);
// Functional to evaluate volume and surface integrals of given integrand
real volume;
auto compute = [&](int q) {
volume = 0.0;
// compute volume integral over {phi < 0} using AutoMixed strategy
if (!halfFaces.empty()) {
ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
if (bernstein::evalBernsteinPoly(phi, x) <= 0) return;
uvector<real, N> trueX = xmin + x * (xmax - xmin);
// uvector<real, N> trueX = x;
// std::cout << "trueX: " << trueX << std::endl;
bool in = true;
for (int i = 0; i < halfFaces.size(); ++i) {
uvector<real, 3> factors(halfFaces[i][0], halfFaces[i][1], halfFaces[i][2]);
if (sum(factors * trueX) + halfFaces[i][3] >= 0) {
in = false;
break;
}
}
if (in) volume += w * integrand(xmin + x * (xmax - xmin));
});
} else
ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
if (bernstein::evalBernsteinPoly(phi, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
});
// scale appropriately
volume *= pow(xmax - xmin, N);
};
compute(q);
auto tAfter = std::chrono::system_clock::now();
std::chrono::duration<double> elapsed_seconds = tAfter - t;
std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
std::cout << "q volume: " << volume << std::endl;
}
// template <typename T, int N>
// void xarrayInit(xarray<T, N>& arr)
// {
// for (auto i = arr.loop(); ~i; ++i) { arr.l(i) = 0; }
// }
void testSpherePowerDirectly()
{
// a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
uvector<real, 3> xmin = -1, xmax = 1;
uvector<real, 3> range = xmax - xmin;
assert(all(range != 0));
// uvector<real, 3> k = 1 / range;
// uvector<real, 3> bias = -k * xmin;
uvector<real, 3> k = range;
uvector<real, 3> bias = xmin;
real a[3] = {1, 1, 1};
real c[3] = {0, 0, 0};
real r = 1;
uvector<int, 3> ext = 3;
xarray<real, 3> phiPower(nullptr, ext), phiBernstein(nullptr, ext);
algoim_spark_alloc(real, phiPower, phiBernstein);
xarrayInit(phiPower);
auto v = xarray2StdVector(phiPower);
auto vv = xarray2StdVector(phiBernstein);
for (int dim = 0; dim < 3; ++dim) {
uvector<int, 3> idx = 0;
phiPower.m(idx) += a[dim] * (bias(dim) - c[dim]) * (bias(dim) - c[dim]);
idx(dim) = 1;
phiPower.m(idx) = 2 * a[dim] * k(dim) * (bias(dim) - c[dim]);
idx(dim) = 2;
phiPower.m(idx) = a[dim] * k(dim) * k(dim);
}
phiPower.m(0) -= r * r;
auto integrand = [](const uvector<real, 3>&) { return 1.0; };
power2BernsteinTensorProduct(phiPower, phiBernstein);
v = xarray2StdVector(phiPower);
uvector<real, 3> testX(0., 0., 0.5);
real testEval = powerEvaluation(phiPower, testX);
std::cout << "eval power:" << testEval << std::endl;
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10, {});
}
void powerTransformation(const uvector<real, 3>& scale, const uvector<real, 3>& bias, xarray<real, 3>& phiPower)
{
std::vector<std::vector<std::vector<real>>> dimOrderExpansion;
const auto& ext = phiPower.ext();
for (int dim = 0; dim < 3; ++dim) {
dimOrderExpansion.push_back(std::vector<std::vector<real>>(ext(dim)));
for (int degree = 0; degree < ext(dim); ++degree) {
const real* binomDegree = Binomial::row(degree);
dimOrderExpansion[dim][degree].reserve(degree + 1);
// 根据二项定理展开
for (int i = 0; i <= degree; ++i) {
dimOrderExpansion[dim][degree].push_back(binomDegree[i] * pow(scale(dim), i) * pow(bias(dim), degree - i));
}
}
}
for (auto i = phiPower.loop(); ~i; ++i) {
// 迭代器必须按照坐标的升序进行访问，即，访问ijk时，(i-1)jk，i(j-1)k，ij(k-1)必须已经访问过
real base = phiPower.l(i);
phiPower.l(i) = 0;
auto exps = i();
for (MultiLoop<3> j(0, exps + 1); ~j; ++j) {
real item = base;
for (int dim = 0; dim < 3; ++dim) { item *= dimOrderExpansion[dim][exps(dim)][j(dim)]; }
phiPower.m(j()) += item;
}
}
}
void testSpherePowerDirectlyByTransformation()
{
// a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
uvector<real, 3> xmin = -1, xmax = 1;
uvector<real, 3> range = xmax - xmin;
assert(all(range != 0));
// uvector<real, 3> k = 1 / range;
// uvector<real, 3> bias = -k * xmin;
real a[3] = {1, 1, 1};
real c[3] = {0, 0, 0};
real r = 0.2;
uvector<int, 3> ext = 3;
xarray<real, 3> phiPower(nullptr, ext), phiBernstein(nullptr, ext);
algoim_spark_alloc(real, phiPower, phiBernstein);
xarrayInit(phiPower);
auto v = xarray2StdVector(phiPower);
auto vv = xarray2StdVector(phiBernstein);
for (int dim = 0; dim < 3; ++dim) {
uvector<int, 3> idx = 0;
phiPower.m(idx) += a[dim] * (c[dim]) * (c[dim]);
idx(dim) = 1;
phiPower.m(idx) = 2 * a[dim] * (-c[dim]);
idx(dim) = 2;
phiPower.m(idx) = a[dim];
}
phiPower.m(0) -= r * r;
powerTransformation(range, xmin, phiPower);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
power2BernsteinTensorProduct(phiPower, phiBernstein);
v = xarray2StdVector(phiPower);
uvector<real, 3> testX(0., 0., 0.5);
real testEval = powerEvaluation(phiPower, testX);
std::cout << "eval power:" << testEval << std::endl;
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10, {});
}
void testCylinderPowerDirectly()
{
// a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
uvector<real, 3> xmin = -1, xmax = 1;
uvector<real, 3> range = xmax - xmin;
assert(all(range != 0));
uvector<real, 3> k = range;
uvector<real, 3> bias = xmin;
real c[3] = {0, 0, 0};
real r = 1, h = 0.8;
uvector<int, 3> ext = 3;
xarray<real, 3> phiPower(nullptr, ext), phiBernstein(nullptr, ext);
algoim_spark_alloc(real, phiPower, phiBernstein);
xarrayInit(phiPower);
auto v = xarray2StdVector(phiPower);
auto vv = xarray2StdVector(phiBernstein);
real top = c[2] + h * 0.5, bottom = c[2] - h * 0.5;
phiPower.m(uvector<int, 3>(2, 0, 2)) = -1;
phiPower.m(uvector<int, 3>(0, 2, 2)) = -1;
phiPower.m(uvector<int, 3>(2, 0, 1)) = top + bottom;
phiPower.m(uvector<int, 3>(0, 2, 1)) = top + bottom;
phiPower.m(uvector<int, 3>(2, 0, 0)) = -bottom * top;
phiPower.m(uvector<int, 3>(0, 2, 0)) = -bottom * top;
phiPower.m(uvector<int, 3>(1, 0, 2)) = 2 * c[0];
phiPower.m(uvector<int, 3>(0, 1, 2)) = 2 * c[1];
phiPower.m(uvector<int, 3>(1, 0, 1)) = -2 * c[0] * (top + bottom);
phiPower.m(uvector<int, 3>(0, 1, 1)) = -2 * c[1] * (top + bottom);
phiPower.m(uvector<int, 3>(1, 0, 0)) = 2 * c[0] * (top * bottom);
phiPower.m(uvector<int, 3>(0, 1, 0)) = 2 * c[0] * (top * bottom);
phiPower.m(uvector<int, 3>(0, 0, 2)) = -(c[0] * c[0] + c[1] * c[1] - r * r);
phiPower.m(uvector<int, 3>(0, 0, 1)) = (bottom + top) * (c[0] * c[0] + c[1] * c[1] - r * r);
phiPower.m(uvector<int, 3>(0, 0, 0)) = -bottom * top * (c[0] * c[0] + c[1] * c[1] - r * r);
powerTransformation(range, xmin, phiPower);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
power2BernsteinTensorProduct(phiPower, phiBernstein);
v = xarray2StdVector(phiPower);
uvector<real, 3> testX(0., 0., 0.75);
real testEval = powerEvaluation(phiPower, testX);
std::cout << "eval power:" << testEval << std::endl;
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
std::cout << "testCylinderPowerDirectly" << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10,
{
{0, 0, 1, -top },
{0, 0, -1, bottom}
});
}
void testCylinderAlgoim()
{
const int N = 3;
const uvector<int, N> P = 3;
real xmin = -1, xmax = 1;
std::vector<xarray<real, N>> phis(3, xarray<real, N>(nullptr, P));
auto fphiCircle = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return x * x + y * y - 1;
};
auto fphiUp = [](const uvector<real, 3>& xx) {
real z = xx(2);
return z - 0.5;
};
auto fphiFloor = [](const uvector<real, 3>& xx) {
real z = xx(2);
return -z - 0.5;
};
algoim_spark_alloc(real, phis[0]);
algoim_spark_alloc(real, phis[1]);
algoim_spark_alloc(real, phis[2]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiCircle(xmin + x * (xmax - xmin)); },
phis[0]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiUp(xmin + x * (xmax - xmin)); }, phis[1]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiFloor(xmin + x * (xmax - xmin)); }, phis[2]);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
// // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
std::cout << "testCylinderAlgoim" << std::endl;
qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}
void test6PlanesAlgoim()
{
const int N = 3;
const uvector<int, N> P = 3;
real xmin = -1, xmax = 1;
std::vector<xarray<real, N>> phis(6, xarray<real, N>(nullptr, P));
// auto fphiPlane1 = [](const uvector<real, 3>& xx) {
// real x = xx(0);
// real y = xx(1);
// real z = xx(2);
// return x + y + z - 1;
// };
// auto fphiPlane2 = [](const uvector<real, 3>& xx) {
// real x = xx(0);
// real y = xx(1);
// real z = xx(2);
// return x - y + z - 1;
// };
// auto fphiPlane3 = [](const uvector<real, 3>& xx) {
// real x = xx(0);
// real y = xx(1);
// real z = xx(2);
// return -x + y + z - 1;
// };
auto fphiPlane1 = [](const uvector<real, 3>& xx) {
real x = xx(0);
return x - 0.5;
};
auto fphiPlane2 = [](const uvector<real, 3>& xx) {
real x = xx(0);
return -x - 0.5;
};
auto fphiPlane3 = [](const uvector<real, 3>& xx) {
real y = xx(1);
return y - 0.5;
};
auto fphiPlane4 = [](const uvector<real, 3>& xx) {
real y = xx(1);
return -y - 0.5;
};
auto fphiPlane5 = [](const uvector<real, 3>& xx) {
real z = xx(2);
return z - 0.5;
};
auto fphiPlane6 = [](const uvector<real, 3>& xx) {
real z = xx(2);
return -z - 0.5;
};
algoim_spark_alloc(real, phis[0]);
algoim_spark_alloc(real, phis[1]);
algoim_spark_alloc(real, phis[2]);
algoim_spark_alloc(real, phis[3]);
algoim_spark_alloc(real, phis[4]);
algoim_spark_alloc(real, phis[5]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
phis[0]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
phis[1]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
phis[2]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane4(xmin + x * (xmax - xmin)); },
phis[3]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane5(xmin + x * (xmax - xmin)); },
phis[4]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane6(xmin + x * (xmax - xmin)); },
phis[5]);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
// // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
std::cout << "testPlanesAlgoim" << std::endl;
qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}
void test8PlanesAlgoim()
{
const int N = 3;
const uvector<int, N> P = 3;
real xmin = -1, xmax = 1;
std::vector<xarray<real, N>> phis(8, xarray<real, N>(nullptr, P));
auto fphiPlane1 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return x + y + z - 1;
};
auto fphiPlane2 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return x - y + z - 1;
};
auto fphiPlane3 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return -x + y + z - 1;
};
auto fphiPlane4 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return -x - y + z - 1;
};
auto fphiPlane5 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return x + y - z - 1;
};
auto fphiPlane6 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return x - y - z - 1;
};
auto fphiPlane7 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return -x + y - z - 1;
};
auto fphiPlane8 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return -x - y - z - 1;
};
algoim_spark_alloc(real, phis[0]);
algoim_spark_alloc(real, phis[1]);
algoim_spark_alloc(real, phis[2]);
algoim_spark_alloc(real, phis[3]);
algoim_spark_alloc(real, phis[4]);
algoim_spark_alloc(real, phis[5]);
algoim_spark_alloc(real, phis[6]);
algoim_spark_alloc(real, phis[7]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
phis[0]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
phis[1]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
phis[2]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane4(xmin + x * (xmax - xmin)); },
phis[3]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane5(xmin + x * (xmax - xmin)); },
phis[4]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane6(xmin + x * (xmax - xmin)); },
phis[5]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane7(xmin + x * (xmax - xmin)); },
phis[6]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane8(xmin + x * (xmax - xmin)); },
phis[7]);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
// // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
std::cout << "testPlanesAlgoim" << std::endl;
qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}
void compositePower(const std::vector<xarray<real, 3>>& powers,
int powerIdx,
uvector<int, 3> powerSum,
real factor,
xarray<real, 3>& res)
{
if (powerIdx == 0) {
{
uvector3 ext(1, 1, 1);
for (auto& t : powers) { ext += t.ext() - 1; }
assert(all(ext == res.ext()));
}
xarrayInit(res);
}
if (powerIdx == powers.size()) {
res.m(powerSum) += factor;
return;
}
auto& power = powers[powerIdx];
for (auto i = power.loop(); ~i; ++i) {
if (power.l(i) == 0) {
// factor = 0;
continue;
}
compositePower(powers, powerIdx + 1, powerSum + i(), factor * power.l(i), res);
}
// int a = 1;
}
void test8PlanesPowerDirectly()
{
// a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
uvector<real, 3> xmin = -1, xmax = 1;
uvector<real, 3> range = xmax - xmin;
assert(all(range != 0));
uvector<real, 3> k = range;
uvector<real, 3> bias = xmin;
uvector<int, 3> ext = 2;
const int planeNum = 8;
std::vector<xarray<real, 3>> phiPowers(planeNum, xarray<real, 3>(nullptr, ext));
algoim_spark_alloc(real, phiPowers[0], phiPowers[1], phiPowers[2], phiPowers[3], phiPowers[4], phiPowers[5], phiPowers[6],
phiPowers[7]);
for (int i = 0; i < planeNum; i++) {
xarrayInit(phiPowers[i]);
if (i & 1) {
phiPowers[i].m(uvector<int, 3>(0, 0, 1)) = 1;
} else {
phiPowers[i].m(uvector<int, 3>(0, 0, 1)) = -1;
}
if (i & 2) {
phiPowers[i].m(uvector<int, 3>(0, 1, 0)) = 1;
} else {
phiPowers[i].m(uvector<int, 3>(0, 1, 0)) = -1;
}
if (i & 4) {
phiPowers[i].m(uvector<int, 3>(1, 0, 0)) = 1;
} else {
phiPowers[i].m(uvector<int, 3>(1, 0, 0)) = -1;
}
phiPowers[i].m(uvector<int, 3>(0, 0, 0)) = -1;
}
uvector<int, 3> resExt = 1 + planeNum;
xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
algoim_spark_alloc(real, phiPowerAll, phiBernstein);
compositePower(phiPowers, 0, 0, 1, phiPowerAll);
powerTransformation(range, xmin, phiPowerAll);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
power2BernsteinTensorProduct(phiPowerAll, phiBernstein);
uvector<real, 3> testX(0., 0., 0.75);
real testEval = powerEvaluation(phiPowerAll, testX);
std::cout << "eval power:" << testEval << std::endl;
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
std::cout << "test8PlanesPowerDirectly" << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10,
{
{1, 1, 1, -1},
{1, 1, -1, -1},
{1, -1, 1, -1},
{1, -1, -1, -1},
{-1, 1, 1, -1},
{-1, 1, -1, -1},
{-1, -1, 1, -1},
{-1, -1, -1, -1}
});
}
void test3PlanesAlgoim()
{
const int N = 3;
const uvector<int, N> P = 3;
real xmin = -1, xmax = 1;
std::vector<xarray<real, N>> phis(3, xarray<real, N>(nullptr, P));
auto fphiPlane1 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
return x + y;
};
auto fphiPlane2 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
return x - y;
};
auto fphiPlane3 = [](const uvector<real, 3>& xx) {
real x = xx(0);
return -x - 1;
};
algoim_spark_alloc(real, phis[0]);
algoim_spark_alloc(real, phis[1]);
algoim_spark_alloc(real, phis[2]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
phis[0]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
phis[1]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
phis[2]);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
// // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
std::cout << "testPlanesAlgoim" << std::endl;
qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}
void test3PlanesPowerDirectly()
{
// a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
uvector<real, 3> xmin = -1, xmax = 1;
uvector<real, 3> range = xmax - xmin;
assert(all(range != 0));
uvector<real, 3> k = range;
uvector<real, 3> bias = xmin;
uvector<int, 3> ext = 2;
const int planeNum = 3;
std::vector<xarray<real, 3>> phiPowers(planeNum, xarray<real, 3>(nullptr, ext));
algoim_spark_alloc(real, phiPowers[0]);
algoim_spark_alloc(real, phiPowers[1]);
algoim_spark_alloc(real, phiPowers[2]);
for (int i = 0; i < planeNum; i++) { xarrayInit(phiPowers[i]); }
phiPowers[0].m(uvector<int, 3>(1, 0, 0)) = 1;
phiPowers[0].m(uvector<int, 3>(0, 1, 0)) = 1;
phiPowers[1].m(uvector<int, 3>(1, 0, 0)) = 1;
phiPowers[1].m(uvector<int, 3>(0, 1, 0)) = -1;
phiPowers[2].m(uvector<int, 3>(1, 0, 0)) = -1;
phiPowers[2].m(uvector<int, 3>(0, 0, 0)) = -1;
uvector<int, 3> resExt = 1 + planeNum;
xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
algoim_spark_alloc(real, phiPowerAll);
algoim_spark_alloc(real, phiBernstein);
compositePower(phiPowers, 0, 0, 1, phiPowerAll);
auto v = xarray2StdVector(phiPowerAll);
uvector<real, 3> testX(0., 0.75, 0);
real testEval = powerEvaluation(phiPowerAll, testX);
std::cout << "eval power before trans:" << testEval << std::endl;
powerTransformation(range, xmin, phiPowerAll);
auto vAfterTrans = xarray2StdVector(phiPowerAll);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
power2BernsteinTensorProduct(phiPowerAll, phiBernstein);
// bernstein::normalise(phiBernstein);
real testEval0 = powerEvaluation(phiPowers[0], testX);
std::cout << "eval power0:" << testEval0 << std::endl;
real testEval1 = powerEvaluation(phiPowers[1], testX);
std::cout << "eval power1:" << testEval1 << std::endl;
real testEval2 = powerEvaluation(phiPowers[2], testX);
std::cout << "eval power2:" << testEval2 << std::endl;
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
std::cout << "test8PlanesPowerDirectly" << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10,
{
{1, 1, 0, 0 },
{1, -1, 0, 0 },
{-1, 0, 0, -1},
});
}
void testConicAlgoim()
{
const int N = 3;
const uvector<int, N> P = 3;
real xmin = -1, xmax = 1;
std::vector<xarray<real, N>> phis(3, xarray<real, N>(nullptr, P));
auto fphiPlane1 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return x * x + y * y - z * z;
};
auto fphiPlane2 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
return x * z - y;
};
auto fphiPlane3 = [](const uvector<real, 3>& xx) {
real x = xx(0);
return -x - 1;
};
algoim_spark_alloc(real, phis[0]);
algoim_spark_alloc(real, phis[1]);
algoim_spark_alloc(real, phis[2]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
phis[0]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
phis[1]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
phis[2]);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
// // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
std::cout << "testPlanesAlgoim" << std::endl;
qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}
void testQuarterSphere()
{
uvector<real, 3> xmin = -1, xmax = 1;
uvector<real, 3> range = xmax - xmin;
assert(all(range != 0));
uvector<real, 3> k = range;
uvector<real, 3> bias = xmin;
std::vector<xarray<real, 3>> phiPowers;
phiPowers.push_back(xarray<real, 3>(nullptr, 3));
phiPowers.push_back(xarray<real, 3>(nullptr, 2));
phiPowers.push_back(xarray<real, 3>(nullptr, 2));
algoim_spark_alloc(real, phiPowers[0]);
algoim_spark_alloc(real, phiPowers[1]);
algoim_spark_alloc(real, phiPowers[2]);
for (int i = 0; i < 3; i++) { xarrayInit(phiPowers[i]); }
phiPowers[0].m(uvector<int, 3>(2, 0, 0)) = 1;
phiPowers[0].m(uvector<int, 3>(0, 2, 0)) = 1;
phiPowers[0].m(uvector<int, 3>(0, 0, 2)) = 1;
phiPowers[0].m(uvector<int, 3>(0, 0, 0)) = -1;
phiPowers[1].m(uvector<int, 3>(1, 0, 0)) = 1;
phiPowers[1].m(uvector<int, 3>(0, 1, 0)) = 1;
// phiPowers[1].m(uvector<int, 3>(0, 0, 1)) = 1;
phiPowers[2].m(uvector<int, 3>(1, 0, 0)) = 1;
phiPowers[2].m(uvector<int, 3>(0, 1, 0)) = -1;
// phiPowers[2].m(uvector<int, 3>(0, 0, 1)) = -1;
uvector<int, 3> resExt = 1 + 2 + 1 + 1;
xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
algoim_spark_alloc(real, phiPowerAll);
algoim_spark_alloc(real, phiBernstein);
compositePower(phiPowers, 0, 0, 1, phiPowerAll);
auto v = xarray2StdVector(phiPowerAll);
uvector<real, 3> testX(0., 0.75, 0);
real testEval = powerEvaluation(phiPowerAll, testX);
std::cout << "eval power before trans:" << testEval << std::endl;
powerTransformation(range, xmin, phiPowerAll);
auto vAfterTrans = xarray2StdVector(phiPowerAll);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
power2BernsteinTensorProduct(phiPowerAll, phiBernstein);
// bernstein::normalise(phiBernstein);
real testEval0 = powerEvaluation(phiPowers[0], testX);
std::cout << "eval power0:" << testEval0 << std::endl;
real testEval1 = powerEvaluation(phiPowers[1], testX);
std::cout << "eval power1:" << testEval1 << std::endl;
real testEval2 = powerEvaluation(phiPowers[2], testX);
std::cout << "eval power2:" << testEval2 << std::endl;
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
std::cout << "testSphereQuarterDirectly" << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 30,
{
{1, 1, 0, 0},
{1, -1, 0, 0},
});
}
void test4PlanesAlgoim()
{
const int N = 3;
const uvector<int, N> P = 2;
real xmin = -1, xmax = 1;
std::vector<xarray<real, N>> phis(4, xarray<real, N>(nullptr, P));
auto fphiPlane1 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
return x + y;
};
auto fphiPlane2 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
return x - y;
};
auto fphiPlane3 = [](const uvector<real, 3>& xx) {
real x = xx(0);
return -x - 1;
};
auto fphiPlane4 = [](const uvector<real, 3>& xx) {
real x = xx(0);
real y = xx(1);
real z = xx(2);
// return x + y - z;
return y - 1;
};
algoim_spark_alloc(real, phis[0]);
algoim_spark_alloc(real, phis[1]);
algoim_spark_alloc(real, phis[2]);
algoim_spark_alloc(real, phis[3]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
phis[0]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
phis[1]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
phis[2]);
bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane4(xmin + x * (xmax - xmin)); },
phis[3]);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
// // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
std::cout << "test4PlanesAlgoim" << std::endl;
qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}
void test4PlanesPowerDirectly()
{
// a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
uvector<real, 3> xmin = -1, xmax = 1;
uvector<real, 3> range = xmax - xmin;
assert(all(range != 0));
uvector<real, 3> k = range;
uvector<real, 3> bias = xmin;
uvector<int, 3> ext = 2;
const int planeNum = 4;
std::vector<xarray<real, 3>> phiPowers(planeNum, xarray<real, 3>(nullptr, ext));
algoim_spark_alloc(real, phiPowers[0]);
algoim_spark_alloc(real, phiPowers[1]);
algoim_spark_alloc(real, phiPowers[2]);
algoim_spark_alloc(real, phiPowers[3]);
for (int i = 0; i < planeNum; i++) { xarrayInit(phiPowers[i]); }
phiPowers[0].m(uvector<int, 3>(1, 0, 0)) = 1;
phiPowers[0].m(uvector<int, 3>(0, 1, 0)) = 1;
phiPowers[1].m(uvector<int, 3>(1, 0, 0)) = 1;
phiPowers[1].m(uvector<int, 3>(0, 1, 0)) = -1;
phiPowers[2].m(uvector<int, 3>(1, 0, 0)) = -1;
phiPowers[2].m(uvector<int, 3>(0, 0, 0)) = -1;
phiPowers[3].m(uvector<int, 3>(1, 0, 0)) = 1;
phiPowers[3].m(uvector<int, 3>(0, 1, 0)) = 1;
phiPowers[3].m(uvector<int, 3>(0, 0, 1)) = -1;
uvector<int, 3> resExt = 1 + planeNum;
xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
algoim_spark_alloc(real, phiPowerAll);
algoim_spark_alloc(real, phiBernstein);
compositePower(phiPowers, 0, 0, 1, phiPowerAll);
auto v = xarray2StdVector(phiPowerAll);
uvector<real, 3> testX(0., 0.75, 0.2);
real testEval = powerEvaluation(phiPowerAll, testX);
std::cout << "eval power before trans:" << testEval << std::endl;
powerTransformation(range, xmin, phiPowerAll);
real testEvalAfterTrans = powerEvaluation(phiPowerAll, testX);
std::cout << "eval power after trans:" << testEvalAfterTrans << std::endl;
auto vAfterTrans = xarray2StdVector(phiPowerAll);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
power2BernsteinTensorProduct(phiPowerAll, phiBernstein);
// bernstein::normalise(phiBernstein);
real testEval0 = powerEvaluation(phiPowers[0], testX);
std::cout << "eval power0:" << testEval0 << std::endl;
real testEval1 = powerEvaluation(phiPowers[1], testX);
std::cout << "eval power1:" << testEval1 << std::endl;
real testEval2 = powerEvaluation(phiPowers[2], testX);
std::cout << "eval power2:" << testEval2 << std::endl;
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
std::cout << "test4PlanesPowerDirectly" << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10,
{
{1, 1, 0, 0 },
{1, -1, 0, 0 },
{-1, 0, 0, -1},
{1, 1, -1, 0 },
});
}
void test2PlanesPowerDirectly()
{
// a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
uvector<real, 3> xmin = -1, xmax = 1;
uvector<real, 3> range = xmax - xmin;
assert(all(range != 0));
uvector<real, 3> k = range;
uvector<real, 3> bias = xmin;
uvector<int, 3> ext = 2;
const int planeNum = 2;
std::vector<xarray<real, 3>> phiPowers(planeNum, xarray<real, 3>(nullptr, ext));
algoim_spark_alloc(real, phiPowers[0]);
algoim_spark_alloc(real, phiPowers[1]);
for (int i = 0; i < planeNum; i++) { xarrayInit(phiPowers[i]); }
phiPowers[0].m(uvector<int, 3>(1, 0, 0)) = -1;
phiPowers[0].m(uvector<int, 3>(0, 0, 0)) = -0.8;
phiPowers[1].m(uvector<int, 3>(0, 1, 0)) = 1;
phiPowers[1].m(uvector<int, 3>(0, 0, 0)) = -0.8;
uvector<int, 3> resExt = 1 + planeNum;
xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
algoim_spark_alloc(real, phiPowerAll);
algoim_spark_alloc(real, phiBernstein);
compositePower(phiPowers, 0, 0, 1, phiPowerAll);
auto v = xarray2StdVector(phiPowerAll);
uvector<real, 3> testX(0., 0.75, 0.2);
real testEval = powerEvaluation(phiPowerAll, testX);
std::cout << "eval power before trans:" << testEval << std::endl;
powerTransformation(range, xmin, phiPowerAll);
real testEvalAfterTrans = powerEvaluation(phiPowerAll, testX);
std::cout << "eval power after trans:" << testEvalAfterTrans << std::endl;
auto vAfterTrans = xarray2StdVector(phiPowerAll);
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
power2BernsteinTensorProduct(phiPowerAll, phiBernstein);
// bernstein::normalise(phiBernstein);
real testEval0 = powerEvaluation(phiPowers[0], testX);
std::cout << "eval power0:" << testEval0 << std::endl;
real testEval1 = powerEvaluation(phiPowers[1], testX);
std::cout << "eval power1:" << testEval1 << std::endl;
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
std::cout << "test2PlanesPowerDirectly" << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10,
{
{1, 0, 0, -0.8},
{0, 1, 0, -0.8},
});
}
void testMultiScale()
{
auto phi0 = [](const uvector<real, 3>& xx) {
real x = xx(0) + 100000.;
real y = xx(1);
real z = xx(2) - 100000.;
real r = 800000.;
return x * x + y * y + z * z - r * r;
};
auto phi1 = [](const uvector<real, 3>& xx) {
// real x = xx(0);
// real y = xx(1);
// real z = xx(2);
// return x * x + y * y + z * z - 1;
real x = xx(0) + 100000.;
real y = xx(1) - 800000.;
real z = xx(2) - 100000.;
real r = 1000.;
return x * x + y * y + z * z - r * r;
};
auto phi2 = [](const uvector<real, 3>& xx) {
real z = xx(2);
return -z - 100000.;
};
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
qConvMultiPloy<3>(phi0, phi2, -1000000., 1000000., 3, integrand, 20, Difference);
}
void testBitSet()
{
std::bitset<10> set(128);
set.set();
std::cout << set << std::endl;
}
void testBooluarray()
{
algoim::booluarray<2, 3> tmp;
tmp(0) = true;
tmp(1) = true;
tmp(2) = true;
std::cout << tmp.bits << std::endl;
}
void testPower2Bernstein()
{
uvector<real, 3> ext(3, 3, 3);
xarray<real, 3> power(nullptr, ext), bernstein(nullptr, ext);
algoim_spark_alloc(real, power, bernstein);
xarrayInit(power);
power.m(uvector<int, 3>(1, 2, 0)) = 3;
power.m(uvector<int, 3>(0, 1, 2)) = 4;
power.m(uvector<int, 3>(1, 2, 2)) = 5;
power.m(uvector<int, 3>(0, 0, 0)) = 3;
uvector3 X(0.8, 0.7, 0.6);
auto evalPower = powerEvaluation(power, X);
power2BernsteinTensorProduct(power, bernstein);
auto evalBernstein = bernstein::evalBernsteinPoly(bernstein, X);
std::cout << "eval power: " << evalPower << std::endl;
std::cout << "eval bernstein: " << evalBernstein << std::endl;
}
void testMain()
{
// testBooluarray();
// testMultiPolys();
// testMultiScale();
// testSpherePowerDirectly();
// testSpherePowerDirectlyByTransformation();
// testCylinderPowerDirectly();
// testCylinderAlgoim();
// test8PlanesAlgoim();
// test8PlanesPowerDirectly();
// test4PlanesPowerDirectly();
test4PlanesAlgoim();
// test2PlanesPowerDirectly();
// testQuarterSphere();
// testPower2Bernstein();
}