GJJ
/
HeteQuadrature
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

plane tests

master
gjj 7 months ago
parent
3e3056f3eb
commit
01623b69e6
4 changed files with 642 additions and 40 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 38
      algoim/organizer/primitive.hpp
  2. 18
      algoim/sparkstack.hpp
  3. 6
      examples/examples_quad_multipoly.cpp
  4. 620
      gjj/myDebug.hpp

38
algoim/organizer/primitive.hpp
View File

@ -160,8 +160,8 @@ public:
~PowerTensor()
{
if (sparkStackPtr) {
sparkStackPtr = nullptr;
algoim_spark_release_heap(sparkStackPtr);
sparkStackPtr = nullptr;
}
}
};
@ -178,14 +178,14 @@ public:
// PowerTensorComplex() {}
PowerTensorComplex(const std::vector<xarray<real, 3>>& ts_) : tensors(ts_)
{
uvector3 ext(1, 1, 1);
for (auto& t : ts_) { ext += t.ext() - 1; }
compositeTensor.ext_ = ext;
sparkStackPtr = algoim_spark_alloc_heap(real, compositeTensor);
detail::compositePower(tensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
}
// PowerTensorComplex(const std::vector<xarray<real, 3>>& ts_) : tensors(ts_)
// {
// uvector3 ext(1, 1, 1);
// for (auto& t : ts_) { ext += t.ext() - 1; }
// compositeTensor.ext_ = ext;
// sparkStackPtr = algoim_spark_alloc_heap(real, compositeTensor);
// // detail::compositePower(tensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
// }
PowerTensorComplex(const std::vector<PowerTensor>& pts_)
{
@ -193,7 +193,7 @@ public:
for (auto& pt : pts_) { ext += pt.tensor.ext() - 1; }
compositeTensor.ext_ = ext;
sparkStackPtr = algoim_spark_alloc_heap(real, compositeTensor);
detail::compositePower(tensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
// detail::compositePower(tensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
}
PowerTensorComplex(const std::vector<PowerTensorComplex>& pcs_)
@ -224,7 +224,7 @@ public:
{
std::vector<PowerTensorComplex> pcs;
pcs.emplace_back(*this);
pcs.emplace_back(PowerTensorComplex({pt.tensor}));
// pcs.emplace_back(PowerTensorComplex({pt.tensor}));
return PowerTensorComplex(pcs);
}
@ -233,7 +233,7 @@ public:
bool isInside(uvector3 p)
{
for (auto& t : tensors) {
if (evalPower(t, p) >= 0) { return false; }
// if (evalPower(t, p) >= 0) { return false; }
}
return true;
}
@ -264,7 +264,7 @@ PowerTensorComplex makeMesh(const std::vector<uvector3>& vertices, const std::ve
PowerTensorComplex pc(pts);
pc.compositeTensor.ext_ = ext;
algoim_spark_alloc(real, pc.compositeTensor);
detail::compositePower(pc.tensors, 0, 0, 1, pc.compositeTensor);
// detail::compositePower(pc.tensors, 0, 0, 1, pc.compositeTensor);
return pc;
}
@ -335,12 +335,12 @@ public:
real eval(uvector3 p) override { return evalBernstein(compositeTensor, p); }
BernsteinTensorComplex(const PowerTensorComplex& pc_) : fromPower(true), tensors(pc_.tensors)
{
compositeTensor.ext_ = pc_.compositeTensor.ext();
algoim_spark_alloc(real, compositeTensor);
detail::power2BernsteinTensor(pc_.compositeTensor, compositeTensor);
};
// BernsteinTensorComplex(const PowerTensorComplex& pc_) : fromPower(true), tensors(pc_.tensors)
// {
// compositeTensor.ext_ = pc_.compositeTensor.ext();
// algoim_spark_alloc(real, compositeTensor);
// detail::power2BernsteinTensor(pc_.compositeTensor, compositeTensor);
// };
bool isInside(uvector3 p)
{

18
algoim/sparkstack.hpp
View File

@ -29,15 +29,15 @@ class SparkStack
<< __FILE__ << '\n';
throw std::bad_alloc();
}
std::cout << "Before alloc, the **ptr pointing to nullptr is not nullptr: " << ptr << std::endl;
*ptr = base() + pos();
std::cout << "pos before += len:" << pos() << std::endl;
// std::cout << "Before alloc, the **ptr pointing to nullptr is not nullptr: " << ptr << std::endl;
*ptr = base() + pos();
// std::cout << "pos before += len:" << pos() << std::endl;
pos() += len;
std::cout << "pos after += len:" << pos() << std::endl;
std::cout << "base: " << base() << std::endl;
std::cout << "ptr: " << ptr << std::endl;
std::cout << "*ptr: " << *ptr << std::endl;
std::cout << "==========" << std::endl;
// std::cout << "pos after += len:" << pos() << std::endl;
// std::cout << "base: " << base() << std::endl;
// std::cout << "ptr: " << ptr << std::endl;
// std::cout << "*ptr: " << *ptr << std::endl;
// std::cout << "==========" << std::endl;
if constexpr (sizeof...(rest) == 0)
return len;
else
@ -100,7 +100,7 @@ public:
~SparkStack()
{
pos() -= len_;
std::cout << "Here!" << std::endl;
// std::cout << "Here!" << std::endl;
}
};

6
examples/examples_quad_multipoly.cpp
View File

@ -16,7 +16,7 @@
#include "xarray.hpp"
#include "myDebug.hpp"
#include "primitiveDebug.hpp"
// #include "primitiveDebug.hpp"
using namespace algoim;
@ -397,8 +397,8 @@ int main(int argc, char* argv[])
}
// module_test();
// testMain();
testPrimitive();
testMain();
// testPrimitive();
return 0;
}
#endif

620
gjj/myDebug.hpp
View File

@ -15,6 +15,7 @@
#include "binomial.hpp"
#include "real.hpp"
#include "sparkstack.hpp"
#include "uvector.hpp"
#include "vector"
#include "xarray.hpp"
@ -150,7 +151,7 @@ void qConvMultiPloy2(real xmin, real xmax, const std::vector<xarray<real, N>>& p
// + x * (xmax - xmin));
bool in = true;
for (int i = 0; i < phis.size(); ++i) {
if (!(bernstein::evalBernsteinPoly(phis[i], x) < 0)) {
if (bernstein::evalBernsteinPoly(phis[i], x) >= 0) {
in = false;
break;
}
@ -222,7 +223,7 @@ void power2BernsteinTensorProduct(const xarray<real, N>& phiPower, xarray<real,
}
}
xarray<real, N> subgrid(nullptr, traverseRange);
// algoim_spark_alloc(real, subgrid);
algoim_spark_alloc(real, subgrid);
for (auto ii = subgrid.loop(); ~ii; ++ii) {
real factor = factorBase;
@ -239,7 +240,7 @@ real powerEvaluation(const xarray<real, N>& phi, const uvector<real, N>& x)
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), exps(dim)); }
for (int dim = 0; dim < N; ++dim) { item *= pow(x(dim), i(dim)); }
res += item;
}
return res;
@ -255,7 +256,7 @@ void qConvBernstein(const xarray<real, N>& phi,
{
auto t = std::chrono::system_clock::now();
uvector<real, 3> testX(0., 0., 0.25);
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;
@ -273,6 +274,7 @@ void qConvBernstein(const xarray<real, N>& phi,
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) {
@ -344,7 +346,7 @@ void testSpherePowerDirectly()
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEval << std::endl;
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10, {});
}
@ -415,7 +417,7 @@ void testSpherePowerDirectlyByTransformation()
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEval << std::endl;
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
qConvBernstein(phiBernstein, -1, 1, integrand, 10, {});
}
@ -428,7 +430,7 @@ void testCylinderPowerDirectly()
uvector<real, 3> k = range;
uvector<real, 3> bias = xmin;
real c[3] = {0, 0, 0};
real r = 1, h = 1;
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);
@ -465,7 +467,7 @@ void testCylinderPowerDirectly()
real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
auto vec = xarray2StdVector(phiBernstein);
std::cout << "eval bernstein:" << testEval << std::endl;
std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
std::cout << "testCylinderPowerDirectly" << std::endl;
qConvBernstein(phiBernstein,
-1,
@ -512,6 +514,595 @@ void testCylinderAlgoim()
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;
std::vector<xarray<real, 3>> phiPowers(8, 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 < 8; 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 + phiPowers.size();
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;
std::vector<xarray<real, 3>> phiPowers(3, 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 < 3; 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 + phiPowers.size();
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,
30,
{
{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 = 3;
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;
};
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;
std::vector<xarray<real, 3>> phiPowers(4, 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 < 4; 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 + phiPowers.size();
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 << "test4PlanesPowerDirectly" << std::endl;
qConvBernstein(phiBernstein,
-1,
1,
integrand,
20,
{
{1, 1, 0, 0 },
{1, -1, 0, 0 },
{-1, 0, 0, -1},
{1, 1, -1, 0 },
});
}
void testMultiScale()
{
auto phi0 = [](const uvector<real, 3>& xx) {
@ -564,7 +1155,18 @@ void testMain()
// testMultiPolys();
// testMultiScale();
// testSpherePowerDirectly();
testSpherePowerDirectlyByTransformation();
// testSpherePowerDirectlyByTransformation();
// testCylinderPowerDirectly();
// testCylinderAlgoim();
// test8PlanesAlgoim();
// test8PlanesPowerDirectly();
// test3PlanesAlgoim();
// test3PlanesPowerDirectly();
test4PlanesPowerDirectly();
// test4PlanesAlgoim();
// testQuarterSphere();
}

Write Preview
Loading…
Cancel
Save