From ac0df77087e867bfe899b751e9d32a91301dcf03 Mon Sep 17 00:00:00 2001 From: gjj Date: Thu, 22 Aug 2024 20:40:04 +0800 Subject: [PATCH] unified interfaces --- algoim/organizer/organizer.hpp | 48 ++- algoim/organizer/primitive.hpp | 573 ++++++++++++++++----------- algoim/sparkstack.hpp | 7 + algoim/xarray.hpp | 2 + examples/examples_quad_general.cpp | 53 ++- examples/examples_quad_multipoly.cpp | 6 +- gjj/myDebug.hpp | 2 + gjj/primitiveDebug.hpp | 47 +-- 8 files changed, 431 insertions(+), 307 deletions(-) diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index 76dcb21..c44060c 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -29,11 +29,11 @@ namespace algoim::Organizer class BasicTask { public: - std::vector primitives; + // std::vector> primitives; - BasicTask(std::vector> ps) {}; + // BasicTask(std::vector> ps) {}; - BasicTask(std::shared_ptr p) + BasicTask(std::shared_ptr p) { int q = 20; real volume = 0; @@ -43,25 +43,45 @@ public: uvector3 range = xmax - xmin; - if (auto pt = std::dynamic_pointer_cast(p)) { - detail::powerTransformation(range, xmin, pt->tensor); - auto primitive = BernsteinTensor(*pt); + if (auto pt = std::dynamic_pointer_cast(p)) { + tensor3 tensor(nullptr, 3); + algoim_spark_alloc(real, tensor); + makeSphere(*pt, tensor); + detail::powerTransformation(range, xmin, tensor); + + tensor3 phi(nullptr, 3); + algoim_spark_alloc(real, phi); + detail::power2BernsteinTensor(tensor, phi); uvector testX(0., 0., 0.25); - real testEvalBernstein = bernstein::evalBernsteinPoly(primitive.tensor, testX); + real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX); // auto vec1 = xarray2StdVector(phi); std::cout << "eval bernstein without interpolation:" << testEvalBernstein << std::endl; - ImplicitPolyQuadrature<3> ipquad(primitive.tensor); + ImplicitPolyQuadrature<3> ipquad(phi); ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { - if (primitive.isInside(x)) volume += w * integrand(xmin + x * (xmax - xmin)); + if (isInsideBernstein(phi, x)) volume += w * integrand(xmin + x * (xmax - xmin)); }); - } else if (auto pc = std::dynamic_pointer_cast(p)) { - detail::powerTransformation(range, xmin, pc->compositeTensor); - auto primitive = BernsteinTensorComplex(*pc); - ImplicitPolyQuadrature<3> ipquad(primitive.compositeTensor); + } else if (auto pt = std::dynamic_pointer_cast(p)) { + const int faceCount = pt->indexInclusiveScan.size(); + assert(faceCount > 1); + std::vector planeTensors(faceCount, tensor3(nullptr, 2)); + algoim_spark_alloc(real, planeTensors); + + tensor3 compositeTensor(nullptr, 1 + faceCount); + algoim_spark_alloc(real, compositeTensor); + + makeMesh(*pt, compositeTensor, planeTensors); + detail::powerTransformation(range, xmin, compositeTensor); + + tensor3 phi(nullptr, 1 + faceCount); + algoim_spark_alloc(real, phi); + detail::power2BernsteinTensor(compositeTensor, phi); + + ImplicitPolyQuadrature<3> ipquad(phi); ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { - if (primitive.isInside(x)) volume += w * integrand(xmin + x * (xmax - xmin)); + auto realX = x * range + xmin; + if (isInsidePowers(planeTensors, realX)) volume += w * integrand(realX); }); } volume *= pow(xmax - xmin, 3); diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 8c36c3f..41496ae 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -1,5 +1,6 @@ #pragma once #include +#include #include #include #include "iostream" @@ -17,36 +18,6 @@ namespace algoim::Organizer namespace detail { -void compositePower(const std::vector>& powers, - int powerIdx, - uvector powerSum, - real factor, - xarray& 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 compositePower(const std::vector>& powers) {} template @@ -113,7 +84,7 @@ class Primitive public: virtual void print() = 0; - virtual real eval(uvector3) { return 0; } + virtual real eval(const uvector3&) { return 0; } }; template @@ -135,216 +106,243 @@ real evalBernstein(const xarray& phi, const uvector& x) return bernstein::evalBernsteinPoly(phi, x); } -class PowerTensor : public Primitive +// class PowerTensor : public Primitive +// { +// public: +// xarray tensor; + +// // SparkStack* sparkStackPtr; + +// void print() override { std::cout << "Power" << std::endl; } + +// real eval(const uvector3& p) override { return evalPower(tensor, p); } + +// // PowerTensor() {} + +// PowerTensor(const xarray& t_) : tensor(t_) {} + +// // const xarray& getTensor() { return tensor; } + +// ~PowerTensor() = default; +// }; + +// class PowerTensorComplex : public Primitive +// { +// public: +// xarray compositeTensor; // 复合后的张量 +// SparkStack* sparkStackPtr; +// std::vector> tensors; // 原始张量 + +// // std::vector powerTensors; // 原始张量 + +// void print() override { std::cout << "PowerTensorComplex" << std::endl; } + +// static void compositePower(const std::vector& powers, +// int powerIdx, +// uvector powerSum, +// real factor, +// xarray& res) + +// { +// // const xarray& tensor +// if (powerIdx == 0) { +// { +// uvector3 ext(1, 1, 1); +// for (auto& t : powers) { ext += t.tensor.ext() - 1; } +// assert(all(ext == res.ext())); +// } +// xarrayInit(res); +// } +// if (powerIdx == powers.size()) { +// res.m(powerSum) += factor; +// return; +// } +// auto& tensor = powers[powerIdx].tensor; +// for (auto i = tensor.loop(); ~i; ++i) { +// if (tensor.l(i) == 0) { +// factor = 0; +// continue; +// } +// compositePower(powers, powerIdx + 1, powerSum + i(), factor * tensor.l(i), res); +// } +// } + +// static void compositePower(const std::vector>& powers, +// int powerIdx, +// uvector powerSum, +// real factor, +// xarray& 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); +// } +// } + +// // PowerTensorComplex() {} + +// PowerTensorComplex(const std::vector>& ts_, xarray& ct_) : tensors(ts_), compositeTensor(ct_) +// { +// uvector3 ext(1); +// for (auto& t : ts_) { ext += t.ext() - 1; } +// assert(all(ext == compositeTensor.ext())); +// compositePower(tensors, 0, uvector3(0), 1, compositeTensor); +// } + +// PowerTensorComplex(const std::vector& pts_, xarray& ct_) : compositeTensor(ct_) +// { +// uvector3 ext(1); +// tensors.resize(pts_.size()); +// for (int i = 0; i < pts_.size(); ++i) { +// tensors[i] = pts_[i].tensor; +// ext += pts_[i].tensor.ext() - 1; +// } +// assert(all(ext == compositeTensor.ext())); +// compositePower(tensors, 0, uvector3(0), 1, compositeTensor); +// } + +// PowerTensorComplex(const std::vector& ptcs_, xarray& ct_) : compositeTensor(ct_) +// { +// std::vector> originCompositeTensors; +// uvector3 ext(1); + +// for (auto& ptc : ptcs_) { +// for (auto& t : ptc.tensors) { tensors.emplace_back(t); } +// originCompositeTensors.push_back(ptc.compositeTensor); +// ext += ptc.compositeTensor.ext() - 1; +// } +// compositePower(originCompositeTensors, 0, uvector3(0, 0, 0), 1, compositeTensor); +// } + +// real eval(const uvector3& p) override { return evalPower(compositeTensor, p); } + +// bool isInside(const uvector3& p) +// { +// for (auto& t : tensors) { +// if (evalPower(t, p) >= 0) { return false; } +// } +// return true; +// } +// }; + +// class BernsteinPrimitive : public Primitive +// { +// }; + +// class BernsteinTensor : public BernsteinPrimitive +// { +// public: +// xarray tensor; + +// void print() override { std::cout << "Bernstein" << std::endl; } + +// real eval(const uvector3& p) override { return evalBernstein(tensor, p); } + +// BernsteinTensor(const PowerTensor& pt_) +// { +// auto v0 = xarray2StdVector(pt_.tensor); + +// tensor.ext_ = pt_.tensor.ext(); +// algoim_spark_alloc(real, tensor); +// auto v1 = xarray2StdVector(pt_.tensor); +// auto v2 = xarray2StdVector(tensor); +// detail::power2BernsteinTensor(pt_.tensor, tensor); + +// uvector3 x(0.2, 0.5, 0.6); +// real BernsteinValue = bernstein::evalBernsteinPoly(tensor, x); +// real PowerValue = evalPower(pt_.tensor, x); +// int a = 1; +// } + +// bool isInside(uvector3 p) { return eval(p) < 0; } +// }; + +// class BernsteinTensorComplex : public BernsteinPrimitive +// { +// public: +// bool fromPower; +// xarray compositeTensor; // 复合后的张量 +// std::vector> tensors; // 原始张量 + +// void print() override { std::cout << "Bernstein Complex" << std::endl; } + +// real eval(const 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); +// // }; + +// bool isInside(uvector3 p) +// { +// if (fromPower) { +// for (auto& t : tensors) { +// if (evalPower(t, p) >= 0) { return false; } +// } +// return true; +// } else { +// for (auto& t : tensors) { +// if (eval(p) >= 0) { return false; } +// } +// return true; +// } +// return true; +// }; +// }; + +bool isInsidePowers(const std::vector& tensors, const uvector3& p) { -public: - xarray tensor; - - // SparkStack* sparkStackPtr; - - void print() override { std::cout << "Power" << std::endl; } - - real eval(uvector3 p) override { return evalPower(tensor, p); } - - // PowerTensor() {} - - PowerTensor(const xarray& t_) : tensor(t_) {} - - // const xarray& getTensor() { return tensor; } - - ~PowerTensor() = default; -}; - -class PowerTensorComplex : public Primitive -{ -public: - xarray compositeTensor; // 复合后的张量 - SparkStack* sparkStackPtr; - // std::vector> tensors; // 原始张量 - std::vector tensors; // 原始张量 - - void print() override { std::cout << "PowerTensorComplex" << std::endl; } - - // PowerTensorComplex() {} - - // PowerTensorComplex(const std::vector>& 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& pts_, xarray& ptc_) : tensors(pts_), compositeTensor(ptc_) - { - uvector3 ext(1); - for (auto& pt : pts_) { ext += pt.tensor.ext() - 1; } - assert(all(ext == compositeTensor.ext())); - detail::compositePower(tensors, 0, uvector3(0, 0, 0), 1, compositeTensor); - } - - PowerTensorComplex(const std::vector& pcs_) - { - for (auto& pc : pcs_) { - for (auto& t : pc.tensors) { tensors.push_back(t); } - } - std::vector> originCompositeTensors; - uvector3 ext(1, 1, 1); - for (auto& pc : pcs_) { - originCompositeTensors.push_back(pc.compositeTensor); - ext += pc.compositeTensor.ext() - 1; - } - compositeTensor.ext_ = ext; - sparkStackPtr = algoim_spark_alloc_heap(real, compositeTensor); - detail::compositePower(originCompositeTensors, 0, uvector3(0, 0, 0), 1, compositeTensor); - } - - PowerTensorComplex add(const PowerTensorComplex& pc) - { - std::vector pcs; - pcs.emplace_back(*this); - pcs.emplace_back(pc); - return PowerTensorComplex(pcs); - } - - PowerTensorComplex add(const PowerTensor& pt) - { - std::vector pcs; - pcs.emplace_back(*this); - // pcs.emplace_back(PowerTensorComplex({pt.tensor})); - return PowerTensorComplex(pcs); - } - - real eval(uvector3 p) override { return evalPower(compositeTensor, p); } - - bool isInside(uvector3 p) - { - for (auto& t : tensors) { - // if (evalPower(t, p) >= 0) { return false; } - } - return true; + for (auto& t : tensors) { + if (evalPower(t, p) >= 0) { return false; } } + return true; }; -PowerTensorComplex makeMesh(const std::vector& vertices, const std::vector>& indices) -{ - uvector3 ext(1 + indices.size()); - // PowerTensorComplex pc(ext); - std::vector pts; - for (const auto& index : indices) { - // xarray tensor(nullptr, ); // 最高1次 - // algoim_spark_alloc(real, tensor); - PowerTensor pt(uvector3(2)); - - uvector3 V01 = vertices[index(1)] - vertices[index(0)]; - uvector3 V02 = vertices[index(2)] - vertices[index(0)]; - uvector3 N = cross(V01, V02); - N /= norm(N); - real d = -dot(N, vertices[index(0)]); - // 法线所指方向为>0区域 - pt.tensor.m(uvector3(0, 0, 0)) = d; - pt.tensor.m(uvector3(1, 0, 0)) = N(0); - pt.tensor.m(uvector3(0, 1, 0)) = N(1); - pt.tensor.m(uvector3(0, 0, 1)) = N(2); - pts.push_back(pt); - } - PowerTensorComplex pc(pts); - pc.compositeTensor.ext_ = ext; - algoim_spark_alloc(real, pc.compositeTensor); - // detail::compositePower(pc.tensors, 0, 0, 1, pc.compositeTensor); - return pc; -} +bool isInsideBernstein(const tensor3& t, const uvector3& p) { return evalBernstein(t, p) < 0; } -PowerTensor makeSphere(const real r, const uvector3& c = 0, const uvector3& a = 1) -{ - uvector ext = 3; - PowerTensor pt(ext); - - for (int dim = 0; dim < 3; ++dim) { - uvector idx = 0; - pt.tensor.m(idx) += c(dim) * c(dim); - idx(dim) = 1; - pt.tensor.m(idx) = 2 * a(dim) * (-c(dim)); - idx(dim) = 2; - pt.tensor.m(idx) = a(dim) * a(dim); - } - pt.tensor.m(0) -= r * r; - return pt; -} +bool isInsidePower(const tensor3& t, const uvector3& p) { return evalPower(t, p) < 0; } -PowerTensor makeCylinder(uvector3 startPt, uvector3 endPt, real r) -{ - // PowerTensor pt; - // TODO: - return PowerTensor({}); -} - -class BernsteinPrimitive : public Primitive +class FRep : public Primitive { -}; +private: + std::function f; -class BernsteinTensor : public BernsteinPrimitive -{ public: - xarray tensor; - - void print() override { std::cout << "Bernstein" << std::endl; } - - real eval(uvector3 p) override { return evalBernstein(tensor, p); } + void print() override { std::cout << "FRep" << std::endl; } - BernsteinTensor(const PowerTensor& pt_) - { - auto v0 = xarray2StdVector(pt_.tensor); - - tensor.ext_ = pt_.tensor.ext(); - algoim_spark_alloc(real, tensor); - auto v1 = xarray2StdVector(pt_.tensor); - auto v2 = xarray2StdVector(tensor); - detail::power2BernsteinTensor(pt_.tensor, tensor); - - uvector3 x(0.2, 0.5, 0.6); - real BernsteinValue = bernstein::evalBernsteinPoly(tensor, x); - real PowerValue = evalPower(pt_.tensor, x); - int a = 1; - } + real eval(const uvector3& p) override { return f(p); } - bool isInside(uvector3 p) { return eval(p) < 0; } + FRep(std::function f_) : f(f_) {} }; -class BernsteinTensorComplex : public BernsteinPrimitive +enum PrimitiveType { Sphere, Cylinder, Cone, Mesh, BRep }; + +class PrimitiveDesc { public: - bool fromPower; - xarray compositeTensor; // 复合后的张量 - std::vector> tensors; // 原始张量 - - void print() override { std::cout << "Bernstein Complex" << std::endl; } - - real eval(uvector3 p) override { return evalBernstein(compositeTensor, p); } + // const static PrimitiveType type; + PrimitiveDesc() = default; - // 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) - { - if (fromPower) { - for (auto& t : tensors) { - if (evalPower(t, p) >= 0) { return false; } - } - return true; - } else { - for (auto& t : tensors) { - if (eval(p) >= 0) { return false; } - } - return true; - } - return true; - }; + virtual void print() {} // 空定义也可以，但是一定要有定义 }; class ParametricSurface @@ -388,27 +386,30 @@ private: int degree; }; -class BRep : public Primitive +class BRepDesc : virtual public PrimitiveDesc { + const static PrimitiveType type = BRep; std::vector vertices; std::vector curves; std::vector surfaces; public: - void print() override { std::cout << "FRep" << std::endl; } + void print() override { std::cout << "BRep Description" << std::endl; } - real eval(uvector3 p) override + real eval(const uvector3& p) { // DOTO: the implicit conversion of Parametric BRep return 0; } - BRep(const std::vector& vs_, const std::vector& cs_, const std::vector& ss_) + BRepDesc(const std::vector& vs_, + const std::vector& cs_, + const std::vector& ss_) : vertices(vs_), curves(cs_), surfaces(ss_) { } - BRep(const BRep& brep) + BRepDesc(const BRepDesc& brep) { vertices = brep.vertices; curves = brep.curves; @@ -416,16 +417,116 @@ public: } }; -class FRep : public Primitive +class SphereDesc : virtual public PrimitiveDesc { -private: - std::function f; +public: + const static PrimitiveType type = Sphere; + real radius; + uvector3 center; + uvector3 amplitude; + + SphereDesc(real r_, const uvector3& c_, const uvector3& a_) : PrimitiveDesc(), radius(r_), center(c_), amplitude(a_) {} + + void print() override { std::cout << "Sphere Description" << std::endl; } +}; + +class CylinderDesc : virtual public PrimitiveDesc +{ + const static PrimitiveType type = Cylinder; + uvector3 node1; + uvector3 node2; + real radius; + + CylinderDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {} + + void print() override { std::cout << "Cylinder Description" << std::endl; } +}; + +class ConeDesc : virtual public PrimitiveDesc +{ + const static PrimitiveType type = Cone; + uvector3 node1; + uvector3 node2; + real radius; + + ConeDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {} + + void print() override { std::cout << "Cone Description" << std::endl; } +}; +class MeshDesc : virtual public PrimitiveDesc +{ public: - void print() override { std::cout << "FRep" << std::endl; } + const static PrimitiveType type = Mesh; + std::vector vertices; + std::vector indices; + std::vector indexInclusiveScan; + + MeshDesc(const std::vector& vertices_, + const std::vector& indices_, + const std::vector& indexInclusiveScan_) + : PrimitiveDesc(), vertices(vertices_), indices(indices_), indexInclusiveScan(indexInclusiveScan_) + { + } - real eval(uvector3 p) override { return f(p); } + void print() override { std::cout << "Mesh Description" << std::endl; } +}; - FRep(std::function f_) : f(f_) {} +void makeMesh(const MeshDesc& mesh, xarray& tensor, std::vector>& planeTensors) +{ + uvector3 ext(1 + mesh.indexInclusiveScan.size()); + assert(all(ext == tensor.ext())); + assert(planeTensors.size() == mesh.indexInclusiveScan.size()); + // for (const auto& index : indices) { + for (int i = 0; i < mesh.indexInclusiveScan.size(); ++i) { + const int indexBeg = i == 0 ? 0 : mesh.indexInclusiveScan[i - 1]; + const int indexSize = mesh.indexInclusiveScan[i] - indexBeg; + assert(indexSize >= 3); + auto& planeTensor = planeTensors[i]; + xarrayInit(planeTensor); + auto& vertices = mesh.vertices; + auto& indices = mesh.indices; + uvector3 V01 = vertices[indices[indexBeg + 1]] - mesh.vertices[indices[indexBeg]]; + uvector3 V02 = vertices[indices[indexBeg + 2]] - mesh.vertices[indices[indexBeg]]; + uvector3 N = cross(V01, V02); + N /= norm(N); + real d = -dot(N, vertices[indices[indexBeg]]); + // 法线所指方向为>0区域 + planeTensor.m(uvector3(0, 0, 0)) = d; + planeTensor.m(uvector3(1, 0, 0)) = N(0); + planeTensor.m(uvector3(0, 1, 0)) = N(1); + planeTensor.m(uvector3(0, 0, 1)) = N(2); + // test other vertices + for (int j = indexBeg + 3; j < mesh.indexInclusiveScan[i]; ++j) { + assert(dot(N, vertices[indices[j]]) + d < std::numeric_limits::epsilon()); + } + } + // compositePower(planeTensors, 0, 0, 1, tensor); + // return PowerTensorComplex(planeTensors, tensor); }; + +void makeSphere(const SphereDesc& sphereDesc, xarray& tensor) +{ + uvector ext = 3; + assert(all(ext == tensor.ext())); + + for (int dim = 0; dim < 3; ++dim) { + uvector idx = 0; + tensor.m(idx) += sphereDesc.center(dim) * sphereDesc.center(dim); + idx(dim) = 1; + tensor.m(idx) = 2 * sphereDesc.amplitude(dim) * (-sphereDesc.center(dim)); + idx(dim) = 2; + tensor.m(idx) = sphereDesc.amplitude(dim) * sphereDesc.amplitude(dim); + } + tensor.m(0) -= sphereDesc.radius * sphereDesc.radius; + // return PowerTensor(tensor); +}; + +void makeCylinder(xarray& tensor, uvector3 startPt, uvector3 endPt, real r) +{ + // PowerTensor pt; + // TODO: +} + + } // namespace algoim::Organizer \ No newline at end of file diff --git a/algoim/sparkstack.hpp b/algoim/sparkstack.hpp index 5c49db4..1c8239c 100644 --- a/algoim/sparkstack.hpp +++ b/algoim/sparkstack.hpp @@ -96,6 +96,13 @@ public: len_ = (alloc(&a.data_, a.size()) + ...); } + template + explicit SparkStack(std::vector>& as) + { + len_ = 0; + for (auto& a : as) len_ += alloc(&a.data_, a.size()); + } + // Release memory when the SparkStack object goes out of scope ~SparkStack() { diff --git a/algoim/xarray.hpp b/algoim/xarray.hpp index e569a18..c586bab 100644 --- a/algoim/xarray.hpp +++ b/algoim/xarray.hpp @@ -375,6 +375,8 @@ public: xarray& ref() { return *this; } }; +using tensor3 = xarray; + template void xarrayInit(xarray& arr) { diff --git a/examples/examples_quad_general.cpp b/examples/examples_quad_general.cpp index bd0e67f..fafd61a 100644 --- a/examples/examples_quad_general.cpp +++ b/examples/examples_quad_general.cpp @@ -7,25 +7,24 @@ #include #include "quadrature_general.hpp" -template -struct Ellipsoid -{ - template - T operator() (const algoim::uvector& x) const +template +struct Ellipsoid { + template + T operator()(const algoim::uvector& x) const { if constexpr (N == 2) - return x(0)*x(0) + 4.0*x(1)*x(1) - 1.0; + return x(0) * x(0) + 4.0 * x(1) * x(1) - 1.0; else - return x(0)*x(0) + 4.0*x(1)*x(1) + 9.0*x(2)*x(2) - 1.0; + return x(0) * x(0) + 4.0 * x(1) * x(1) + 9.0 * x(2) * x(2) - 1.0; } - template - algoim::uvector grad(const algoim::uvector& x) const + template + algoim::uvector grad(const algoim::uvector& x) const { if constexpr (N == 2) - return algoim::uvector(2.0*x(0), 8.0*x(1)); + return algoim::uvector(2.0 * x(0), 8.0 * x(1)); else - return algoim::uvector(2.0*x(0), 8.0*x(1), 18.0*x(2)); + return algoim::uvector(2.0 * x(0), 8.0 * x(1), 18.0 * x(2)); } }; @@ -39,8 +38,8 @@ int main(int argc, char* argv[]) { std::cout << "Area of a 2D ellipse using automatic subdivision:\n"; Ellipsoid<2> phi; - auto q = algoim::quadGen<2>(phi, algoim::HyperRectangle(-1.1, 1.1), -1, -1, 4); - double area = q([](const auto& x) { return 1.0; }); + auto q = algoim::quadGen<2>(phi, algoim::HyperRectangle(-1.1, 1.1), -1, -1, 4); + double area = q([](const auto& x) { return 1.0; }); std::cout << " computed area = " << area << "\n"; std::cout << " (exact area = 1.5707963267948966)\n\n"; } @@ -49,8 +48,8 @@ int main(int argc, char* argv[]) { std::cout << "Volume of a 3D ellipsoid using automatic subdivision:\n"; Ellipsoid<3> phi; - auto q = algoim::quadGen<3>(phi, algoim::HyperRectangle(-1.1, 1.1), -1, -1, 4); - double volume = q([](const auto& x) { return 1.0; }); + auto q = algoim::quadGen<3>(phi, algoim::HyperRectangle(-1.1, 1.1), -1, -1, 4); + double volume = q([](const auto& x) { return 1.0; }); std::cout << " computed volume = " << volume << "\n"; std::cout << " (exact volume = 0.6981317007977318)\n\n"; } @@ -59,15 +58,15 @@ int main(int argc, char* argv[]) { int n = 16; std::cout << "Area of a 2D ellipse, computed via the cells of a " << n << " by " << n << " Cartesian grid:\n"; - double dx = 2.2 / n; + double dx = 2.2 / n; Ellipsoid<2> phi; - double area = 0.0; - for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) - { - algoim::uvector xmin{-1.1 + i*dx, -1.1 + j*dx}; - algoim::uvector xmax{-1.1 + i*dx + dx, -1.1 + j*dx + dx}; - area += algoim::quadGen<2>(phi, algoim::HyperRectangle(xmin, xmax), -1, -1, 4).sumWeights(); - } + double area = 0.0; + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) { + algoim::uvector xmin{-1.1 + i * dx, -1.1 + j * dx}; + algoim::uvector xmax{-1.1 + i * dx + dx, -1.1 + j * dx + dx}; + area += algoim::quadGen<2>(phi, algoim::HyperRectangle(xmin, xmax), -1, -1, 4).sumWeights(); + } std::cout << " computed area = " << area << "\n"; std::cout << " (exact area = 1.5707963267948966)\n\n"; } @@ -76,8 +75,8 @@ int main(int argc, char* argv[]) { std::cout << "Surface area of a 3D ellipsoid using automatic subdivision:\n"; Ellipsoid<3> phi; - auto q = algoim::quadGen<3>(phi, algoim::HyperRectangle(-1.1, 1.1), 3, -1, 4); - double surface_area = q.sumWeights(); + auto q = algoim::quadGen<3>(phi, algoim::HyperRectangle(-1.1, 1.1), 3, -1, 4); + double surface_area = q.sumWeights(); std::cout << " computed surface area = " << surface_area << "\n"; std::cout << " (exact surface area = 4.4008095646649703)\n\n"; } @@ -85,8 +84,8 @@ int main(int argc, char* argv[]) // Visualisation of a quadrature scheme in ParaView via XML VTP file { std::cout << "Visualisation of a quadrature scheme in ParaView via XML VTP file:\n"; - Ellipsoid<3> phi; - auto q = algoim::quadGen<3>(phi, algoim::HyperRectangle(-1.1, 1.1), -1, -1, 2); + Ellipsoid<3> phi; + auto q = algoim::quadGen<3>(phi, algoim::HyperRectangle(-1.1, 1.1), -1, -1, 2); std::ofstream f("scheme.vtp"); algoim::outputQuadratureRuleAsVtpXML(q, f); std::cout << " scheme.vtp file written, containing " << q.nodes.size() << " quadrature points\n"; diff --git a/examples/examples_quad_multipoly.cpp b/examples/examples_quad_multipoly.cpp index 3430595..d23a1ac 100644 --- a/examples/examples_quad_multipoly.cpp +++ b/examples/examples_quad_multipoly.cpp @@ -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 diff --git a/gjj/myDebug.hpp b/gjj/myDebug.hpp index d9a3203..8f02e53 100644 --- a/gjj/myDebug.hpp +++ b/gjj/myDebug.hpp @@ -1190,4 +1190,6 @@ void testMain() // testQuarterSphere(); // testPower2Bernstein(); + int a = log(2.01); + std::cout << "log(3) " << a << std::endl; } diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index 432a049..de0f0c6 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -30,44 +30,37 @@ using namespace algoim; void case0() { - std::vector> ps; + // std::vector> ps; // mesh - std::vector vertices = {uvector3(-0.8, -0.8, -0.8), - uvector3(-0.8, -0.8, 0.8), - uvector3(-0.8, 0.8, -0.8), - uvector3(-0.8, 0.8, 0.8), - uvector3(0.8, -0.8, -0.8), - uvector3(0.8, -0.8, 0.8), - uvector3(0.8, 0.8, -0.8), - uvector3(0.8, 0.8, 0.8) + std::vector vertices = {uvector3(-0.8, -0.8, -0.8), + uvector3(-0.8, -0.8, 0.8), + uvector3(-0.8, 0.8, -0.8), + uvector3(-0.8, 0.8, 0.8), + uvector3(0.8, -0.8, -0.8), + uvector3(0.8, -0.8, 0.8), + uvector3(0.8, 0.8, -0.8), + uvector3(0.8, 0.8, 0.8) }; - std::vector> indixes = { - uvector(2, 1, 0), - uvector(1, 2, 3), - uvector(0, 1, 4), - uvector(4, 1, 5), - uvector(6, 2, 0), - uvector(6, 0, 4), - uvector(0, 3, 2), - uvector(0, 2, 6), - uvector(1, 3, 7), - uvector(1, 7, 5), - uvector(5, 7, 6), - uvector(4, 5, 6), - }; + std::vector indices = {2, 3, 1, 0, // force break here + 4, 0, 1, 5, // + 4, 6, 2, 0, // + 6, 7, 3, 2, // + 7, 6, 4, 5, // + 3, 7, 5, 1}; + std::vector scan = {4, 8, 12, 16, 20, 24}; // ps.emplace_back(std::make_shared(makeMesh(vertices, indixes))); // ps.emplace_back(std::make_shared(makeSphere(0.2))); - auto phi0 = std::make_shared(makeMesh(vertices, indixes)); - phi0->add(*std::make_shared(makeSphere(0.2))); - auto basicTask = BasicTask(phi0); + // ps.emplace_back(std::make_shared(MeshDesc(vertices, indices, scan))); + // auto basicTask = BasicTask(ps); } void case1() { - auto phi0 = std::make_shared(makeSphere(0.2)); + auto phi0 = std::make_shared(SphereDesc(0.8, 0., 1.)); + // SphereDesc sphere(0.8, 0., 1.); auto basicTask = BasicTask(phi0); }