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

fix bug of composite tensor

master
gjj 1 year ago
parent
e5d6a869cc
commit
e7d39058d8
5 changed files with 96 additions and 28 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. 6
      algoim/organizer/organizer.hpp
  2. 6
      algoim/organizer/primitive.hpp
  3. 2
      algoim/sparkstack.hpp
  4. 90
      gjj/myDebug.hpp
  5. 2
      gjj/primitiveDebug.hpp

6
algoim/organizer/organizer.hpp
View File

@ -86,9 +86,11 @@ public:
makeMesh(*pt, compositeTensor, planeTensors);
detail::powerTransformation(range, xmin, compositeTensor);
auto planeStdVector1 = xarray2StdVector(planeTensors[0]);
auto planeStdVector2 = xarray2StdVector(planeTensors[1]);
auto compositeTensorStdVector = xarray2StdVector(compositeTensor);
uvector<real, 3> testX(0.8, 0.8, 0.8);
uvector<real, 3> testX(0., 0.75, 0.2);
real textEvalPower = evalPower(compositeTensor, testX);
tensor3 phi(nullptr, 1 + faceCount);
@ -97,6 +99,8 @@ public:
int quadraturePointCount = 0;
ImplicitPolyQuadrature<3> ipquad(phi);
real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
quadraturePointCount++;
auto realX = x * range + xmin;

6
algoim/organizer/primitive.hpp
View File

@ -131,7 +131,7 @@ static void compositePower(const std::vector<xarray<real, 3>>& powers,
auto& power = powers[powerIdx];
for (auto i = power.loop(); ~i; ++i) {
if (power.l(i) == 0) {
factor = 0;
// factor = 0;
continue;
}
compositePower(powers, powerIdx + 1, powerSum + i(), factor * power.l(i), res);
@ -220,7 +220,7 @@ real evalBernstein(const xarray<real, N>& phi, const uvector<real, N>& x)
// auto& tensor = powers[powerIdx].tensor;
// for (auto i = tensor.loop(); ~i; ++i) {
// if (tensor.l(i) == 0) {
// factor = 0;
// //factor = 0;
// continue;
// }
// compositePower(powers, powerIdx + 1, powerSum + i(), factor * tensor.l(i), res);
@ -249,7 +249,7 @@ real evalBernstein(const xarray<real, N>& phi, const uvector<real, N>& x)
// auto& power = powers[powerIdx];
// for (auto i = power.loop(); ~i; ++i) {
// if (power.l(i) == 0) {
// factor = 0;
// // factor = 0;
// continue;
// }
// compositePower(powers, powerIdx + 1, powerSum + i(), factor * power.l(i), res);

2
algoim/sparkstack.hpp
View File

@ -107,7 +107,7 @@ public:
~SparkStack()
{
pos() -= len_;
std::cout << "Here!" << std::endl;
// std::cout << "Here!" << std::endl;
}
};

90
gjj/myDebug.hpp
View File

@ -273,7 +273,7 @@ void qConvBernstein(const xarray<real, N>& phi,
// 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;
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;
@ -708,7 +708,8 @@ void test8PlanesPowerDirectly()
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));
const int planeNum = 8;
std::vector<xarray<real, 3>> phiPowers(planeNum, xarray<real, 3>(nullptr, ext));
algoim_spark_alloc(real,
phiPowers[0],
phiPowers[1],
@ -718,7 +719,7 @@ void test8PlanesPowerDirectly()
phiPowers[5],
phiPowers[6],
phiPowers[7]);
for (int i = 0; i < 8; i++) {
for (int i = 0; i < planeNum; i++) {
xarrayInit(phiPowers[i]);
if (i & 1) {
phiPowers[i].m(uvector<int, 3>(0, 0, 1)) = 1;
@ -737,7 +738,7 @@ void test8PlanesPowerDirectly()
}
phiPowers[i].m(uvector<int, 3>(0, 0, 0)) = -1;
}
uvector<int, 3> resExt = 1 + phiPowers.size();
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);
@ -819,18 +820,19 @@ void test3PlanesPowerDirectly()
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));
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 < 3; i++) { xarrayInit(phiPowers[i]); }
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 + phiPowers.size();
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);
@ -1042,12 +1044,13 @@ void test4PlanesPowerDirectly()
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));
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 < 4; i++) { xarrayInit(phiPowers[i]); }
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;
@ -1057,7 +1060,7 @@ void test4PlanesPowerDirectly()
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();
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);
@ -1106,6 +1109,68 @@ void test4PlanesPowerDirectly()
});
}
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) {
@ -1183,13 +1248,12 @@ void testMain()
// test8PlanesAlgoim();
// test8PlanesPowerDirectly();
// test4PlanesPowerDirectly();
// test4PlanesAlgoim();
test2PlanesPowerDirectly();
// testQuarterSphere();
// testPower2Bernstein();
int a = log(2.01);
std::cout << "log(3) " << a << std::endl;
}

2
gjj/primitiveDebug.hpp
View File

@ -117,4 +117,4 @@ void case2()
std::vector<int> scan = {4, 8, 12, 16, 20, 24};
}
void testPrimitive() { casePolyhedron2(); }
void testPrimitive() { casePolyhedron1(); }
Write Preview
Loading…
Cancel
Save