diff --git a/examples/examples_mytest.cpp b/examples/examples_mytest.cpp index e69de29..8b13789 100644 --- a/examples/examples_mytest.cpp +++ b/examples/examples_mytest.cpp @@ -0,0 +1 @@ + diff --git a/gjj/myDebug.hpp b/gjj/myDebug.hpp index 817eddd..30863fb 100644 --- a/gjj/myDebug.hpp +++ b/gjj/myDebug.hpp @@ -15,6 +15,7 @@ #include "uvector.hpp" #include "vector" #include "xarray.hpp" +#include using namespace algoim; @@ -87,8 +88,8 @@ void qConvMultiPloy(const F1& fphi1, // bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2); // Build quadrature hierarchy - // ImplicitPolyQuadrature ipquad(phi1, phi2); - ImplicitPolyQuadrature ipquad(phi1); + ImplicitPolyQuadrature ipquad(phi1, phi2); + // ImplicitPolyQuadrature ipquad(phi1); // Functional to evaluate volume and surface integrals of given integrand real volume, surf; @@ -97,8 +98,8 @@ void qConvMultiPloy(const F1& fphi1, surf = 0.0; // compute volume integral over {phi < 0} using AutoMixed strategy ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { - // if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0) - if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin)); + 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& x, real w, const uvector& wn) { @@ -112,7 +113,76 @@ void qConvMultiPloy(const F1& fphi1, std::cout << "q volume: " << volume << std::endl; } -void testMultiPloys() +template +void qConvMultiPloy2( + real xmin, + real xmax, + const uvector& P, + const F& integrand, + int q, + std::string qfile) +{ + // Construct phi by mapping [0,1] onto bounding box [xmin,xmax] + + // bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi1(xmin + x * (xmax - xmin)); }, phi1); + // bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2); + + std::vector> phis; // 大量sphere + 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.1 + rand() % 9 / 10.0; + auto fphi = [&](const uvector& xx) { + real x = xx(0) - centerX; + real y = xx(1) - centerY; + real z = xx(2) - centerZ; + return x * x + y * y + z * z - 1; + }; + xarray phi(nullptr, P); + algoim_spark_alloc(real, phi); + bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi(xmin + x * (xmax - xmin)); }, phi); + } + auto t = std::chrono::system_clock::now(); + + ImplicitPolyQuadrature ipquad(phis); + + + + // Build quadrature hierarchy + // ImplicitPolyQuadrature ipquad(phi1, phi2); + // ImplicitPolyQuadrature 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& 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)); + volume += w * integrand(xmin + x * (xmax - xmin)); + }); + // compute surface integral over {phi == 0} using AutoMixed strategy + // ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& 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 elapsed_seconds = tAfter - t; + std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n"; + std::cout << "q volume: " << volume << std::endl; +} + +void fromCommon2Bernstein() { + +} + +void testMultiPolys() { // Visusalisation of a 2D implicitly-defined domain involving the intersection of two polynomials; this example // corresponds to the top-left example of Figure 15, https://doi.org/10.1016/j.jcp.2021.110720 @@ -124,14 +194,14 @@ void testMultiPloys() return x * x + y * y + z * z - 1; }; auto phi1 = [](const uvector& 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); // real y = xx(1); // real z = xx(2); - // return x * y * z; + // return x * x + y * y + z * z - 1; + real x = xx(0); + real y = xx(1); + real z = xx(2); + return x * y * z; }; // auto phi0 = [](const uvector& xx) { // real x = xx(0) * 2 - 1; @@ -148,7 +218,8 @@ void testMultiPloys() // + x * (-0.6169311810674598 - 0.19449299999999994 * y - 0.005459163675646665 * y * y); // }; auto integrand = [](const uvector& x) { return 1.0; }; - qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC"); + // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC"); + qConvMultiPloy2<3>(-1.0, 1.0, 3, integrand, 20, "exampleC"); std::cout << "\n\nQuadrature visualisation of a 2D implicitly-defined domain involving the\n"; std::cout << "intersection of two polynomials, corresponding to the top-left example of Figure 15,\n"; std::cout << "https://doi.org/10.1016/j.jcp.2021.110720, written to exampleC-xxxx.vtp files\n"; @@ -175,5 +246,5 @@ void testBooluarray() void testMain() { testBooluarray(); - testMultiPloys(); + testMultiPolys(); }