poly tests

.clang-format

@@ -0,0 +1,268 @@
+---
+BasedOnStyle: Google
+AccessModifierOffset: -4
+AlignAfterOpenBracket: Align
+AlignArrayOfStructures: Left
+AlignConsecutiveAssignments:
+  Enabled: true
+  AcrossEmptyLines: false
+  AcrossComments: true
+  AlignCompound: true
+  PadOperators: true
+AlignConsecutiveBitFields:
+  Enabled: true
+  AcrossEmptyLines: false
+  AcrossComments: true
+  AlignCompound: true
+  PadOperators: true
+AlignConsecutiveDeclarations:
+  Enabled: true
+  AcrossEmptyLines: false
+  AcrossComments: true
+  AlignCompound: true
+  PadOperators: true
+AlignConsecutiveMacros:
+  Enabled: true
+  AcrossEmptyLines: false
+  AcrossComments: true
+  AlignCompound: true
+  PadOperators: true
+AlignConsecutiveShortCaseStatements:
+  Enabled: true
+  AcrossEmptyLines: false
+  AcrossComments: true
+  AlignCaseColons: false
+AlignEscapedNewlines: Left
+AlignOperands: Align
+AlignTrailingComments:
+  Kind: Always
+  OverEmptyLines: 2
+AllowAllArgumentsOnNextLine: false
+AllowAllParametersOfDeclarationOnNextLine: false
+AllowShortBlocksOnASingleLine: Always
+AllowShortCaseLabelsOnASingleLine: true
+AllowShortEnumsOnASingleLine: true
+AllowShortFunctionsOnASingleLine: All
+AllowShortIfStatementsOnASingleLine: WithoutElse
+AllowShortLambdasOnASingleLine: All
+AllowShortLoopsOnASingleLine: true
+AlwaysBreakAfterDefinitionReturnType: None
+AlwaysBreakAfterReturnType: None
+AlwaysBreakBeforeMultilineStrings: false
+AlwaysBreakTemplateDeclarations: Yes
+AttributeMacros:
+  - __capability
+BinPackArguments: false
+BinPackParameters: false
+BitFieldColonSpacing: Both
+BraceWrapping:
+  AfterCaseLabel: false
+  AfterClass: false
+  AfterControlStatement: Never
+  AfterEnum: false
+  AfterFunction: false
+  AfterNamespace: false
+  AfterObjCDeclaration: false
+  AfterStruct: false
+  AfterUnion: false
+  AfterExternBlock: false
+  BeforeCatch: false
+  BeforeElse: false
+  BeforeLambdaBody: false
+  BeforeWhile: false
+  IndentBraces: false
+  SplitEmptyFunction: true
+  SplitEmptyRecord: true
+  SplitEmptyNamespace: true
+BreakAfterAttributes: Never
+BreakAfterJavaFieldAnnotations: false
+BreakArrays: true
+BreakBeforeBinaryOperators: NonAssignment
+BreakBeforeBraces: Linux
+BreakBeforeConceptDeclarations: Always
+BreakBeforeInlineASMColon: OnlyMultiline
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializers: BeforeColon
+BreakInheritanceList: AfterComma
+BreakStringLiterals: true
+ColumnLimit: 128
+CommentPragmas: "^(!<)? IWYU pragma:"
+CompactNamespaces: false
+ConstructorInitializerIndentWidth: 4
+ContinuationIndentWidth: 4
+Cpp11BracedListStyle: true
+DerivePointerAlignment: true
+DisableFormat: false
+EmptyLineAfterAccessModifier: Never
+EmptyLineBeforeAccessModifier: LogicalBlock
+ExperimentalAutoDetectBinPacking: false
+FixNamespaceComments: true
+ForEachMacros:
+  - foreach
+  - Q_FOREACH
+  - BOOST_FOREACH
+IfMacros:
+  - KJ_IF_MAYBE
+IncludeBlocks: Regroup
+IncludeCategories:
+  - Regex: ^<ext/.*\.h>
+    Priority: 2
+    SortPriority: 0
+    CaseSensitive: false
+  - Regex: ^<.*\.h>
+    Priority: 1
+    SortPriority: 0
+    CaseSensitive: false
+  - Regex: ^<.*
+    Priority: 2
+    SortPriority: 0
+    CaseSensitive: false
+  - Regex: .*
+    Priority: 3
+    SortPriority: 0
+    CaseSensitive: false
+IncludeIsMainRegex: ([-_](test|unittest))?$
+IncludeIsMainSourceRegex: ""
+IndentAccessModifiers: false
+IndentCaseBlocks: false
+IndentCaseLabels: true
+IndentExternBlock: AfterExternBlock
+IndentGotoLabels: false
+IndentPPDirectives: None
+IndentRequiresClause: true
+IndentWidth: 4
+IndentWrappedFunctionNames: true
+InsertBraces: false
+InsertNewlineAtEOF: false
+InsertTrailingCommas: None
+IntegerLiteralSeparator:
+  Binary: 0
+  BinaryMinDigits: 0
+  Decimal: 0
+  DecimalMinDigits: 0
+  Hex: 0
+  HexMinDigits: 0
+JavaScriptQuotes: Leave
+JavaScriptWrapImports: true
+KeepEmptyLinesAtEOF: false
+KeepEmptyLinesAtTheStartOfBlocks: false
+LambdaBodyIndentation: Signature
+Language: Cpp
+LineEnding: DeriveLF
+MacroBlockBegin: ""
+MacroBlockEnd: ""
+MaxEmptyLinesToKeep: 2
+NamespaceIndentation: None
+ObjCBinPackProtocolList: Never
+ObjCBlockIndentWidth: 2
+ObjCBreakBeforeNestedBlockParam: true
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: true
+PPIndentWidth: -1
+PackConstructorInitializers: NextLine
+PenaltyBreakAssignment: 2
+PenaltyBreakBeforeFirstCallParameter: 1
+PenaltyBreakComment: 300
+PenaltyBreakFirstLessLess: 120
+PenaltyBreakOpenParenthesis: 0
+PenaltyBreakString: 1000
+PenaltyBreakTemplateDeclaration: 10
+PenaltyExcessCharacter: 1000000
+PenaltyIndentedWhitespace: 0
+PenaltyReturnTypeOnItsOwnLine: 200
+PointerAlignment: Right
+QualifierAlignment: Leave
+RawStringFormats:
+  - Language: Cpp
+    Delimiters:
+      - cc
+      - CC
+      - cpp
+      - Cpp
+      - CPP
+      - c++
+      - C++
+    CanonicalDelimiter: ""
+    BasedOnStyle: google
+  - Language: TextProto
+    Delimiters:
+      - pb
+      - PB
+      - proto
+      - PROTO
+    EnclosingFunctions:
+      - EqualsProto
+      - EquivToProto
+      - PARSE_PARTIAL_TEXT_PROTO
+      - PARSE_TEST_PROTO
+      - PARSE_TEXT_PROTO
+      - ParseTextOrDie
+      - ParseTextProtoOrDie
+      - ParseTestProto
+      - ParsePartialTestProto
+    CanonicalDelimiter: pb
+    BasedOnStyle: google
+ReferenceAlignment: Pointer
+ReflowComments: true
+RemoveBracesLLVM: false
+RemoveParentheses: Leave
+RemoveSemicolon: false
+RequiresClausePosition: OwnLine
+RequiresExpressionIndentation: OuterScope
+SeparateDefinitionBlocks: Always
+ShortNamespaceLines: 0
+SortIncludes: Never
+SortJavaStaticImport: Before
+SortUsingDeclarations: LexicographicNumeric
+SpaceAfterCStyleCast: false
+SpaceAfterLogicalNot: false
+SpaceAfterTemplateKeyword: true
+SpaceAroundPointerQualifiers: Default
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeCaseColon: false
+SpaceBeforeCpp11BracedList: false
+SpaceBeforeCtorInitializerColon: true
+SpaceBeforeInheritanceColon: true
+SpaceBeforeJsonColon: false
+SpaceBeforeParens: Custom
+SpaceBeforeParensOptions:
+  AfterControlStatements: true
+  AfterForeachMacros: true
+  AfterFunctionDeclarationName: false
+  AfterFunctionDefinitionName: false
+  AfterIfMacros: true
+  AfterOverloadedOperator: false
+  AfterRequiresInClause: false
+  AfterRequiresInExpression: true
+  BeforeNonEmptyParentheses: false
+SpaceBeforeRangeBasedForLoopColon: true
+SpaceBeforeSquareBrackets: false
+SpaceInEmptyBlock: false
+SpacesBeforeTrailingComments: 1
+SpacesInAngles: Never
+SpacesInContainerLiterals: true
+SpacesInLineCommentPrefix:
+  Minimum: 1
+  Maximum: -1
+SpacesInParens: Never
+SpacesInParensOptions:
+  InConditionalStatements: false
+  InCStyleCasts: false
+  InEmptyParentheses: false
+  Other: false
+SpacesInSquareBrackets: false
+Standard: Auto
+StatementAttributeLikeMacros:
+  - Q_EMIT
+StatementMacros:
+  - Q_UNUSED
+  - QT_REQUIRE_VERSION
+TabWidth: 8
+UseTab: Never
+VerilogBreakBetweenInstancePorts: true
+WhitespaceSensitiveMacros:
+  - BOOST_PP_STRINGIZE
+  - CF_SWIFT_NAME
+  - NS_SWIFT_NAME
+  - PP_STRINGIZE
+  - STRINGIZE

CMakeLists.txt

@@ -5,7 +5,7 @@ set(CMAKE_STANDARD 17)
 
 set(CMAKE_EXPORT_COMPILECOMMANDS ON)
 
-include_directories(algoim)
+include_directories(algoim gjj)
 
 # find_package(LAPACK REQUIRED)
 find_path(LAPACKE_INCLUDE_DIR NAMES lapacke.h PATH_SUFFIXES include)

algoim/bernstein.hpp

@@ -297,7 +297,7 @@ namespace algoim::bernstein
         for (int i = 1; i < P; ++i)
             for (int j = P - 1; j >= i; --j)
             {
-                alpha.a(j) *= tau;
+                alpha.a(j) *= tau;  // here the ptr to the array is included in the return
                 alpha.a(j) += alpha.a(j - 1) * (1.0 - tau);
             }
     }

algoim/booluarray.hpp

@@ -10,20 +10,16 @@ namespace algoim
 {
 namespace booluarray_detail
 {
-        constexpr int pow(int base, int exp)
-        {
-            return exp == 0 ? 1 : base * pow(base, exp - 1);
-        }
+constexpr int pow(int base, int exp) { return exp == 0 ? 1 : base * pow(base, exp - 1); }
 
 template <int E, int N>
-        constexpr int furl(const uvector<int,N>& i)
+constexpr int furl(const uvector<int, N>& i) // 为什么这里可以是constexpr，uvector不是变量吗？
 {
     int ind = i(0);
-            for (int j = 1; j < N; ++j)
-                ind = ind * E + i(j);
+    for (int j = 1; j < N; ++j) ind = ind * E + i(j);
     return ind;
 }
-    }
+} // namespace booluarray_detail
 
 // booluarray implements a simple N-dimensional array of booleans, with
 // compile-time extent E across all dimensions; it is essentially a basic
@@ -32,31 +28,24 @@ namespace algoim
 class booluarray
 {
     constexpr static int size = booluarray_detail::pow(E, N);
-        std::bitset<size> bits;
+
+
 public:
+    std::bitset<size> bits;
+
     booluarray() {}
 
     booluarray(bool val)
     {
-            if (val)
-                bits.set();
+        if (val) bits.set(); // 全部置为1
     }
 
-        bool operator() (const uvector<int,N>& i) const
-        {
-            return bits[booluarray_detail::furl<E>(i)];
-        }
+    bool operator()(const uvector<int, N>& i) const { return bits[booluarray_detail::furl<E>(i)]; }
 
-        auto operator() (const uvector<int,N>& i)
-        {
-            return bits[booluarray_detail::furl<E>(i)];
-        }
+    auto operator()(const uvector<int, N>& i) { return bits[booluarray_detail::furl<E>(i)]; }
 
     // returns true iff the entire array is false
-        bool none() const
-        {
-            return bits.none();
-        }
+    bool none() const { return bits.none(); }
 };
 } // namespace algoim
 

algoim/polyset.hpp

@@ -6,20 +6,21 @@
 #include <cassert>
 #include <vector>
 #include "booluarray.hpp"
+#include "real.hpp"
+#include "xarray.hpp"
 
 namespace algoim
 {
 // PolySet implements a simple container to hold one or more Bernstein polynomials
 // and their associated masks
 template <int N, int E>
-    struct PolySet
-    {
-        struct Poly
-        {
+struct PolySet {
+    struct Poly {
         uvector<int, N>  ext;    // Degree/extent of polynomial
         size_t           offset; // Offset into buffer, storing the xarray<Real,N> polynomial data
         booluarray<N, E> mask;   // Mask
     };
+
     std::vector<real> buff;      // Memory buffer containing polynomial data
     std::vector<Poly> items;     // Record of contained polynomials
 
@@ -45,10 +46,7 @@ namespace algoim
         poly(items.size() - 1) = p;
     }
 
-        size_t count() const
-        {
-            return items.size();
-        }
+    size_t count() const { return items.size(); }
 };
 } // namespace algoim
 

algoim/quadrature_multipoly.hpp

@@ -16,6 +16,7 @@
 // binary value 0 or 1 and indicates whether its accompanying polynomial is provably
 // nonzero on that subcell, or if its roots can be ignored. Typically, M = 4 or 8 is
 // a good choice, and this is specified by the following macro def.
+#include <vector>
 #define ALGOIM_M 8
 
 #include <algorithm>
@@ -35,46 +36,39 @@ namespace algoim
 namespace detail
 {
 template <int N>
-        booluarray<N,ALGOIM_M> mask_driver(const xarray<real,N>& f, const booluarray<N,ALGOIM_M>& fmask, const xarray<real,N>* g, const booluarray<N,ALGOIM_M>* gmask)
+booluarray<N, ALGOIM_M> mask_driver(const xarray<real, N>&         f,
+                                    const booluarray<N, ALGOIM_M>& fmask,
+                                    const xarray<real, N>*         g,
+                                    const booluarray<N, ALGOIM_M>* gmask)
 {
     booluarray<N, ALGOIM_M> mask(false);
-            auto helper = [&](auto&& self, uvector<int,N> a, uvector<int,N> b)
-            {
+    auto                    helper = [&](auto&& self, uvector<int, N> a, uvector<int, N> b) {
         bool overlap = false;
         for (MultiLoop<N> i(a, b); ~i; ++i)
-                    if (fmask(i()) && (!gmask || (*gmask)(i())))
-                        overlap = true;
-                if (!overlap)
-                    return;
+            if (fmask(i()) && (!gmask || (*gmask)(i()))) overlap = true;
+        if (!overlap) return;
 
         real             eps = 0.015625 / ALGOIM_M;
         uvector<real, N> xa, xb;
-                for (int dim = 0; dim < N; ++dim)
-                {
+        for (int dim = 0; dim < N; ++dim) {
             xa(dim) = real(a(dim)) / ALGOIM_M - eps;
             xb(dim) = real(b(dim)) / ALGOIM_M + eps;
         }
 
-                if (g)
-                {
+        if (g) {
             xarray<real, N> fab(nullptr, f.ext()), gab(nullptr, g->ext());
             algoim_spark_alloc(real, fab, gab);
             bernstein::deCasteljau(f, xa, xb, fab);
             bernstein::deCasteljau(*g, xa, xb, gab);
-                    if (bernstein::orthantTest(fab, gab))
-                        return;
-                }
-                else
-                {
+            if (bernstein::orthantTest(fab, gab)) return;
+        } else {
             xarray<real, N> fab(nullptr, f.ext());
             algoim_spark_alloc(real, fab);
             bernstein::deCasteljau(f, xa, xb, fab);
-                    if (bernstein::uniformSign(fab) != 0)
-                        return;
+            if (bernstein::uniformSign(fab) != 0) return;
         }
 
-                if (b(0) - a(0) == 1)
-                {
+        if (b(0) - a(0) == 1) {
             assert(all(b - a == 1));
             assert(fmask(a) && (!gmask || (*gmask)(a)));
             mask(a) = true;
@@ -83,8 +77,7 @@ namespace algoim
 
         assert(all(b - a > 1) && all((b - a) % 2 == 0));
         uvector<int, N> delta = (b - a) / 2;
-                for (MultiLoop<N> i(0,2); ~i; ++i)
-                    self(self, a + i() * delta, a + (i() + 1) * delta);
+        for (MultiLoop<N> i(0, 2); ~i; ++i) self(self, a + i() * delta, a + (i() + 1) * delta);
     };
     helper(helper, 0, ALGOIM_M);
     return mask;
@@ -96,7 +89,10 @@ namespace algoim
 // definitively do not share common zeros in that subrectangle; if the returned mask is true,
 // then shared zeros may exist (and with high likelihood)
 template <int N>
-        booluarray<N,ALGOIM_M> intersectionMask(const xarray<real,N>& f, const booluarray<N,ALGOIM_M>& fmask, const xarray<real,N>& g, const booluarray<N,ALGOIM_M>& gmask)
+booluarray<N, ALGOIM_M> intersectionMask(const xarray<real, N>&         f,
+                                         const booluarray<N, ALGOIM_M>& fmask,
+                                         const xarray<real, N>&         g,
+                                         const booluarray<N, ALGOIM_M>& gmask)
 {
     return mask_driver(f, fmask, &g, &gmask);
 }
@@ -117,8 +113,7 @@ namespace algoim
 {
     booluarray<N - 1, ALGOIM_M> r(false);
     for (MultiLoop<N> i(0, ALGOIM_M); ~i; ++i)
-                if (mask(i()))
-                    r(remove_component(i(), k)) = true;
+        if (mask(i())) r(remove_component(i(), k)) = true;
     return r;
 }
 
@@ -146,23 +141,19 @@ namespace algoim
 bool lineIntersectsMask(const booluarray<N, ALGOIM_M>& mask, const uvector<real, N - 1>& x, int k)
 {
     using std::floor;
-            if constexpr (N > 1)
-            {
+    if constexpr (N > 1) {
         uvector<int, N> cell;
         for (int dim = 0; dim < N; ++dim)
             if (dim < k)
                 cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast<int>(floor(x(dim) * ALGOIM_M))));
             else if (dim > k)
                 cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast<int>(floor(x(dim - 1) * ALGOIM_M))));
-                for (int i = 0; i < ALGOIM_M; ++i)
-                {
+        for (int i = 0; i < ALGOIM_M; ++i) {
             cell(k) = i;
-                    if (mask(cell))
-                        return true;
+            if (mask(cell)) return true;
         }
         return false;
-            }
-            else
+    } else
         return !maskEmpty(mask);
 }
 
@@ -172,11 +163,9 @@ namespace algoim
     assert(0 <= k && k < N && (side == 0 || side == 1));
     assert(all(out.ext() == remove_component(a.ext(), k)));
     int P = a.ext(k);
-            for (auto i = out.loop(); ~i; ++i)
-            {
+    for (auto i = out.loop(); ~i; ++i) {
         uvector<int, N> j;
-                for (int dim = 0; dim < N; ++dim)
-                    j(dim) = (dim < k) ? i(dim) : ( (dim == k) ? side*(P-1) : i(dim - 1) );
+        for (int dim = 0; dim < N; ++dim) j(dim) = (dim < k) ? i(dim) : ((dim == k) ? side * (P - 1) : i(dim - 1));
         out.l(i) = a.m(j);
     }
 }
@@ -186,11 +175,9 @@ namespace algoim
 {
     assert(0 <= k && k < N && (side == 0 || side == 1));
     booluarray<N - 1, ALGOIM_M> r;
-            for (MultiLoop<N-1> i(0,ALGOIM_M); ~i; ++i)
-            {
+    for (MultiLoop<N - 1> i(0, ALGOIM_M); ~i; ++i) {
         uvector<int, N> j;
-                for (int dim = 0; dim < N; ++dim)
-                    j(dim) = (dim < k)? i(dim) : ( (dim == k) ? side*(ALGOIM_M - 1) : i(dim - 1) );
+        for (int dim = 0; dim < N; ++dim) j(dim) = (dim < k) ? i(dim) : ((dim == k) ? side * (ALGOIM_M - 1) : i(dim - 1));
         r(i()) = a(j);
     }
     return r;
@@ -203,39 +190,31 @@ namespace algoim
 T det_qr(xarray<T, 2>& A, int& rank, T tol = 10.0)
 {
     assert(A.ext(0) == A.ext(1) && A.ext(0) > 0);
-            using std::max;
     using std::abs;
+    using std::max;
     T   det        = 1.0;
     int n          = A.ext(0);
     T   max_diag_r = 0.0;
-            for (int j = 0; j < n; ++j)
-            {
+    for (int j = 0; j < n; ++j) {
         T   m = -1;
         int k = -1;
-                for (int i = j; i < n; ++i)
-                {
+        for (int i = j; i < n; ++i) {
             T mag = 0;
-                    for (int a = 0; a < n; ++a)
-                        mag += util::sqr(A(a,i));
-                    if (k == -1 || mag >= m)
-                    {
+            for (int a = 0; a < n; ++a) mag += util::sqr(A(a, i));
+            if (k == -1 || mag >= m) {
                 m = mag;
                 k = i;
             }
         }
         assert(j <= k && k < n);
-                if (k != j)
-                {
-                    for (int a = 0; a < n; ++a)
-                        std::swap(A(a,j), A(a,k));
+        if (k != j) {
+            for (int a = 0; a < n; ++a) std::swap(A(a, j), A(a, k));
             det *= -1.0;
         }
-                for (int i = n - 1; i >= j + 1; --i)
-                {
+        for (int i = n - 1; i >= j + 1; --i) {
             T c, s;
             util::givens_get(A(i - 1, j), A(i, j), c, s);
-                    for (int k = j; k < n; ++k)
-                        util::givens_rotate(A(i-1,k), A(i,k), c, s);
+            for (int k = j; k < n; ++k) util::givens_rotate(A(i - 1, k), A(i, k), c, s);
         }
         det        *= A(j, j);
         max_diag_r  = max(max_diag_r, abs(A(j, j)));
@@ -244,8 +223,7 @@ namespace algoim
     tol  *= max_diag_r * n * std::numeric_limits<T>::epsilon();
     rank  = 0;
     for (int i = 0; i < n; ++i)
-                if (abs(A(i,i)) > tol)
-                    ++rank;
+        if (abs(A(i, i)) > tol) ++rank;
 
     return det;
 }
@@ -255,8 +233,7 @@ namespace algoim
 uvector<int, N - 1> resultantExtent(const uvector<int, N>& p, const uvector<int, N>& q, int dim)
 {
     uvector<int, N - 1> ext;
-            for (int i = 0; i < N - 1; ++i)
-            {
+    for (int i = 0; i < N - 1; ++i) {
         int ii = (i < dim) ? i : i + 1;
         ext(i) = (p(dim) - 1) * (q(ii) - 1) + (q(dim) - 1) * (p(ii) - 1) + 1;
     }
@@ -268,8 +245,7 @@ namespace algoim
 uvector<int, N - 1> discriminantExtent(const uvector<int, N>& p, int dim)
 {
     uvector<int, N - 1> ext;
-            for (int i = 0; i < N - 1; ++i)
-            {
+    for (int i = 0; i < N - 1; ++i) {
         int ii = (i < dim) ? i : i + 1;
         ext(i) = (2 * p(dim) - 3) * (p(ii) - 1) + 1;
     }
@@ -291,19 +267,19 @@ namespace algoim
     int             P = p.ext(k);
     int             Q = q ? q->ext(k) : P - 1;
     int             M = (P == Q) ? P - 1 : P + Q - 2;
+    xarray<real, 2> mat(nullptr, uvector<int, 2>{M, M});
+
     assert(P >= 2 && Q >= 1 && M >= 1);
 
     xarray<real, N - 1> f(nullptr, out.ext());
-            xarray<real,2> mat(nullptr, uvector<int,2>{M, M});
+
     real *pk, *qk;
     algoim_spark_alloc(real, f, mat);
     algoim_spark_alloc(real, &pk, P, &qk, Q);
 
-            for (auto i = f.loop(); ~i; ++i)
-            {
+    for (auto i = f.loop(); ~i; ++i) {
         uvector<real, N - 1> x;
-                for (int dim = 0; dim < N - 1; ++dim)
-                    x(dim) = bernstein::modifiedChebyshevNode(i(dim), f.ext(dim));
+        for (int dim = 0; dim < N - 1; ++dim) x(dim) = bernstein::modifiedChebyshevNode(i(dim), f.ext(dim));
 
         bernstein::collapseAlongAxis(p, x, k, pk);
         if (q)
@@ -329,8 +305,7 @@ namespace algoim
 
     // If able to reduce the degree, recompute the resultant on the lower-degree poly space
     // which is expected to have better conditioning
-            if (b)
-                resultant_core(p, q, k, out);
+    if (b) resultant_core(p, q, k, out);
     return true;
 }
 
@@ -361,20 +336,17 @@ namespace algoim
     assert(psi.count() == 0);
 
     // For every phi(i) ...
-            for (int i = 0; i < phi.count(); ++i)
-            {
+    for (int i = 0; i < phi.count(); ++i) {
         const auto& p    = phi.poly(i);
         const auto& mask = phi.mask(i);
 
         // Examine bottom and top faces in the k'th dimension
-                for (int side = 0; side <= 1; ++side)
-                {
+        for (int side = 0; side <= 1; ++side) {
             xarray<real, N - 1> p_face(nullptr, remove_component(p.ext(), k));
             algoim_spark_alloc(real, p_face);
             restrictToFace(p, k, side, p_face);
             auto p_face_mask = nonzeroMask(p_face, restrictToFace(mask, k, side));
-                    if (!maskEmpty(p_face_mask))
-                    {
+            if (!maskEmpty(p_face_mask)) {
                 bernstein::autoReduction(p_face);
                 bernstein::normalise(p_face);
                 psi.push_back(p_face, p_face_mask);
@@ -386,14 +358,12 @@ namespace algoim
         algoim_spark_alloc(real, p_k);
         bernstein::elevatedDerivative(p, k, p_k);
         auto disc_mask = intersectionMask(p, mask, p_k, mask);
-                if (!maskEmpty(disc_mask))
-                {
+        if (!maskEmpty(disc_mask)) {
             // note: computed disc might have lower degree than the following
             uvector<int, N - 1> R = discriminantExtent(p.ext(), k);
             xarray<real, N - 1> disc(nullptr, R);
             algoim_spark_alloc(real, disc);
-                    if (discriminant(p, k, disc))
-                    {
+            if (discriminant(p, k, disc)) {
                 bernstein::normalise(disc);
                 psi.push_back(disc, collapseMask(disc_mask, k));
             }
@@ -401,21 +371,19 @@ namespace algoim
     }
 
     // Consider every pairwise combination of resultants ...
-            for (int i = 0; i < phi.count(); ++i) for (int j = i + 1; j < phi.count(); ++j)
-            {
+    for (int i = 0; i < phi.count(); ++i)
+        for (int j = i + 1; j < phi.count(); ++j) {
             const auto& p     = phi.poly(i);
             const auto& pmask = phi.mask(i);
             const auto& q     = phi.poly(j);
             const auto& qmask = phi.mask(j);
             auto        mask  = intersectionMask(p, pmask, q, qmask);
-                if (!maskEmpty(mask))
-                {
+            if (!maskEmpty(mask)) {
                 // note: computed resultant might have lower degree than the following
                 uvector<int, N - 1> R = resultantExtent(p.ext(), q.ext(), k);
                 xarray<real, N - 1> res(nullptr, R);
                 algoim_spark_alloc(real, res);
-                    if (resultant(p, q, k, res))
-                    {
+                if (resultant(p, q, k, res)) {
                     bernstein::normalise(res);
                     psi.push_back(res, collapseMask(mask, k));
                 }
@@ -432,31 +400,27 @@ namespace algoim
     uvector<real, N> s = 0;
     has_disc           = false;
     // For every phi(i) ...
-            for (int i = 0; i < phi.count(); ++i)
-            {
+    for (int i = 0; i < phi.count(); ++i) {
         const auto& p    = phi.poly(i);
         const auto& mask = phi.mask(i);
 
         // Accumulate onto score by sampling at midpoint of every subcell of mask
-                for (MultiLoop<N> j(0,ALGOIM_M); ~j; ++j) if (mask(j()))
-                {
+        for (MultiLoop<N> j(0, ALGOIM_M); ~j; ++j)
+            if (mask(j())) {
                 uvector<real, N> x   = (j() + 0.5) / real(ALGOIM_M);
                 uvector<real, N> g   = bernstein::evalBernsteinPolyGradient(p, x);
                 real             sum = 0;
-                    for (int dim = 0; dim < N; ++dim)
-                    {
+                for (int dim = 0; dim < N; ++dim) {
                     g(dim)  = abs(g(dim));
                     sum    += g(dim);
                 }
-                    if (sum > 0)
-                        s += g / sum;
+                if (sum > 0) s += g / sum;
             }
 
         // Consider discriminant
         xarray<real, N> p_k(nullptr, p.ext());
         algoim_spark_alloc(real, p_k);
-                for (int k = 0; k < N; ++k)
-                {
+        for (int k = 0; k < N; ++k) {
             bernstein::elevatedDerivative(p, k, p_k);
             auto disc_mask  = intersectionMask(p, mask, p_k, mask);
             has_disc(k)    |= !maskEmpty(disc_mask);
@@ -490,8 +454,7 @@ namespace algoim
 enum QuadStrategy { AlwaysGL, AlwaysTS, AutoMixed };
 
 template <int N>
-    struct ImplicitPolyQuadrature
-    {
+struct ImplicitPolyQuadrature {
     enum IntegralType { Inner, OuterSingle, OuterAggregate };
 
     PolySet<N, ALGOIM_M>          phi; // Given N-dimensional polynomials
@@ -499,9 +462,10 @@ namespace algoim
     // 这种struct的递归定义一般是不允许的，因为会内存爆炸，但是模板类可以通过特化定义递归的终止条件
     // 所以有template<> struct ImplicitPolyQuadrature<0> {};
     ImplicitPolyQuadrature<N - 1> base; // Base polynomials corresponding to removal of axis k
-        bool auto_apply_TS;                                                     // If quad method is auto chosen, indicates whether TS is applied  // what is ts ? 
+    bool         auto_apply_TS; // If quad method is auto chosen, indicates whether TS is applied  // what is ts ?  tanh-sinh?
     IntegralType type;          // Whether an inner integral, or outer of two kinds
-        std::array<std::tuple<int,ImplicitPolyQuadrature<N-1>>,N-1> base_other; // Stores other base cases, besides k, when in aggregate mode
+    std::array<std::tuple<int, ImplicitPolyQuadrature<N - 1>>, N - 1>
+        base_other;             // Stores other base cases, besides k, when in aggregate mode
 
     // Default ctor sets to an uninitialised state
     ImplicitPolyQuadrature() : k(-1) {}
@@ -510,8 +474,7 @@ namespace algoim
     ImplicitPolyQuadrature(const xarray<real, N>& p)
     {
         auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true));
-            if (!detail::maskEmpty(mask))
-                phi.push_back(p, mask);
+        if (!detail::maskEmpty(mask)) phi.push_back(p, mask);
         build(true, false);
     }
 
@@ -519,33 +482,42 @@ namespace algoim
     ImplicitPolyQuadrature(const xarray<real, N>& p, const xarray<real, N>& q)
     {
         {
-								auto tmp = booluarray<N,ALGOIM_M>(false);
+            // auto tmp  = booluarray<N, ALGOIM_M>(false);
             auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true));
-								tmp(0) = true;
-								auto res = !detail::maskEmpty(tmp);
-                if (!detail::maskEmpty(mask))
-                    phi.push_back(p, mask);
+            // tmp(0)    = true;
+            // auto res  = !detail::maskEmpty(tmp);  // 似乎是测试代码？
+            if (!detail::maskEmpty(mask)) phi.push_back(p, mask);
         }
         {
             auto mask = detail::nonzeroMask(q, booluarray<N, ALGOIM_M>(true));
-                if (!detail::maskEmpty(mask))
-                    phi.push_back(q, mask);
+            if (!detail::maskEmpty(mask)) phi.push_back(q, mask);
+        }
+        build(true, false);
+    }
+
+    // Build quadrature hierarchy for a domain implicitly defined by multiple polynomials
+    ImplicitPolyQuadrature(const std::vector<xarray<real, N>>& ps)
+    {
+        for (const auto& p : ps) {
+            auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true));
+            if (!detail::maskEmpty(mask)) phi.push_back(p, mask);
         }
         build(true, false);
     }
 
     // Build quadrature hierarchy for a given domain implicitly defined by two polynomials with user-defined masks
-        ImplicitPolyQuadrature(const xarray<real,N>& p, const booluarray<N,ALGOIM_M>& pmask, const xarray<real,N>& q, const booluarray<N,ALGOIM_M>& qmask)
+    ImplicitPolyQuadrature(const xarray<real, N>&         p,
+                           const booluarray<N, ALGOIM_M>& pmask,
+                           const xarray<real, N>&         q,
+                           const booluarray<N, ALGOIM_M>& qmask)
     {
         {
             auto mask = detail::nonzeroMask(p, pmask);
-                if (!maskEmpty(mask))
-                    phi.push_back(p, mask);
+            if (!maskEmpty(mask)) phi.push_back(p, mask);
         }
         {
             auto mask = detail::nonzeroMask(q, qmask);
-                if (!maskEmpty(mask))
-                    phi.push_back(q, mask);
+            if (!maskEmpty(mask)) phi.push_back(q, mask);
         }
         build(true, false);
     }
@@ -557,22 +529,18 @@ namespace algoim
         this->auto_apply_TS = auto_apply_TS;
 
         // If phi is empty, apply a tensor-product Gaussian quadrature
-            if (phi.count() == 0)
-            {
+        if (phi.count() == 0) {
             k                   = N;
             this->auto_apply_TS = false;
             return;
         }
 
-            if constexpr (N == 1)
-            {
+        if constexpr (N == 1) {
             // If in one dimension, there is only one choice of height direction and
             // the recursive process halts
             k = 0;
             return;
-            }
-            else
-            {
+        } else {
             // Compute score; penalise any directions which likely contain vertical tangents
             uvector<bool, N> has_disc;
             uvector<real, N> score = detail::score_estimate(phi, has_disc);
@@ -581,8 +549,7 @@ namespace algoim
             assert(max(abs(score)) > 0);
             score /= 2 * max(abs(score));
             for (int i = 0; i < N; ++i)
-                    if (!has_disc(i))
-                        score(i) += 1.0;
+                if (!has_disc(i)) score(i) += 1.0;
 
             // Choose height direction and form base polynomials; if tanh-sinh is being used at this
             // level, suggest the same all the way down; moreover, suggest tanh-sinh if a non-empty
@@ -593,11 +560,10 @@ namespace algoim
 
             // If this is the outer integral, and surface quadrature schemes are required, apply
             // the dimension-aggregated scheme when necessary
-                if (outer && has_disc(k))
-                {
+            if (outer && has_disc(k)) {
                 type = OuterAggregate;
-                    for (int i = 0; i < N; ++i) if (i != k)
-                    {
+                for (int i = 0; i < N; ++i)
+                    if (i != k) {
                         auto& [kother, base] = base_other[i < k ? i : i - 1];
                         kother               = i;
                         detail::eliminate_axis(phi, kother, base.phi);
@@ -617,15 +583,12 @@ namespace algoim
 
         // If there are no interfaces, apply rectangular tensor-product Gauss-Legendre quadrature
         // 所以这个是积分域包裹了整个box时的做法？
-            if (k == N)
-            {
+        if (k == N) {
             assert(!auto_apply_TS);
-                for (MultiLoop<N> i(0, q); ~i; ++i)
-                {
+            for (MultiLoop<N> i(0, q); ~i; ++i) {
                 uvector<real, N> x;
                 real             w = 1.0;
-                    for (int dim = 0; dim < N; ++dim)
-                    {
+                for (int dim = 0; dim < N; ++dim) {
                     x(dim)  = GaussQuad::x(q, i(dim)); // 0 <= i(dim) < q
                     w      *= GaussQuad::w(q, i(dim));
                 }
@@ -637,12 +600,10 @@ namespace algoim
         // Determine maximum possible number of roots; used to allocate a small buffer
         // ???
         int max_count = 2;
-            for (size_t i = 0; i < phi.count(); ++i)
-                max_count += phi.poly(i).ext(k) - 1;
+        for (size_t i = 0; i < phi.count(); ++i) max_count += phi.poly(i).ext(k) - 1;
 
         // Base integral invokes the following integrand
-            auto integrand = [&](const uvector<real,N-1>& xbase, real w)
-            {
+        auto integrand = [&](const uvector<real, N - 1>& xbase, real w) {
             // Allocate node buffer of sufficient size and initialise with {0, 1}
             real* nodes;
             algoim_spark_alloc(real, &nodes, max_count);
@@ -651,15 +612,13 @@ namespace algoim
             int count = 2;
 
             // For every phi(i) ...
-                for (size_t i = 0; i < phi.count(); ++i)
-                {
+            for (size_t i = 0; i < phi.count(); ++i) {
                 const auto& p    = phi.poly(i);
                 const auto& mask = phi.mask(i);
                 int         P    = p.ext(k);
 
                 // Ignore phi if its mask is void everywhere above the base point
-                    if (!detail::lineIntersectsMask(mask, xbase, k))
-                        continue;
+                if (!detail::lineIntersectsMask(mask, xbase, k)) continue;
 
                 // Restrict polynomial to axis-aligned line and compute its roots
                 real *pline, *roots;
@@ -668,11 +627,9 @@ namespace algoim
                 int rcount = bernstein::bernsteinUnitIntervalRealRoots(pline, P, roots);
 
                 // Add all real roots in [0,1] which are also within masked region of phi
-                    for (int j = 0; j < rcount; ++j)
-                    {
+                for (int j = 0; j < rcount; ++j) {
                     auto x = add_component(xbase, k, roots[j]);
-                        if (detail::pointWithinMask(mask, x))
-                            nodes[count++] = roots[j];
+                    if (detail::pointWithinMask(mask, x)) nodes[count++] = roots[j];
                 }
             };
 
@@ -693,26 +650,23 @@ namespace algoim
             assert(nodes[0] == real(0) && nodes[count - 1] == real(1));
 
             // Apply quadrature to non-degenerate sub-intervals
-                for (int i = 0; i < count - 1; ++i)
-                {
+            for (int i = 0; i < count - 1; ++i) {
                 real x0 = nodes[i];
                 real x1 = nodes[i + 1];
-                    if (x0 == x1)
-                        continue;
+                if (x0 == x1) continue;
 
                 // Choose between Gauss-Legendre and tanh-sinh according to the user-defined strategy
                 bool GL = true;
-                    if (strategy == AlwaysTS)
-                        GL = false;
-                    if (strategy == AutoMixed)
-                        GL = !this->auto_apply_TS;
+                if (strategy == AlwaysTS) GL = false;
+                if (strategy == AutoMixed) GL = !this->auto_apply_TS;
 
                 if (GL)
                     for (int j = 0; j < q; ++j)
                         f(add_component(xbase, k, x0 + (x1 - x0) * GaussQuad::x(q, j)), w * (x1 - x0) * GaussQuad::w(q, j));
                 else
                     for (int j = 0; j < q; ++j)
-                            f(add_component(xbase, k, TanhSinhQuadrature::x(q, j, x0, x1)), w * TanhSinhQuadrature::w(q, j, x0, x1));
+                        f(add_component(xbase, k, TanhSinhQuadrature::x(q, j, x0, x1)),
+                          w * TanhSinhQuadrature::w(q, j, x0, x1));
             }
         };
 
@@ -731,24 +685,20 @@ namespace algoim
         assert(type == OuterSingle || type == OuterAggregate);
 
         // If there is no interface, there is no surface integral
-            if (k == N)
-                return;
+        if (k == N) return;
 
         // Base integral invokes the following integrand which operates in the height direction k_active
         int  k_active  = -1;
-            auto integrand = [&](const uvector<real,N-1>& xbase, real w)
-            {
+        auto integrand = [&](const uvector<real, N - 1>& xbase, real w) {
             assert(0 <= k_active && k_active < N);
             // For every phi(i) ...
-                for (size_t i = 0; i < phi.count(); ++i) 
-                {
+            for (size_t i = 0; i < phi.count(); ++i) {
                 const auto& p    = phi.poly(i);
                 const auto& mask = phi.mask(i);
                 int         P    = p.ext(k_active);
 
                 // Ignore phi if its mask is void everywhere above the base point
-                    if (!detail::lineIntersectsMask(mask, xbase, k_active))
-                        continue;
+                if (!detail::lineIntersectsMask(mask, xbase, k_active)) continue;
 
                 // Compute roots of { x \mapsto phi(xbase + x e_k) }
                 real *pline, *roots;
@@ -758,20 +708,16 @@ namespace algoim
 
                 // Consider all real roots in (0,1) which are also within masked region of phi; evaluate
                 // integrand at interfacial points, multiplying weights by the effective surface Jacobian
-                    for (int j = 0; j < rcount; ++j)
-                    {
+                for (int j = 0; j < rcount; ++j) {
                     auto x = add_component(xbase, k_active, roots[j]);
-                        if (detail::pointWithinMask(mask, x))
-                        {
+                    if (detail::pointWithinMask(mask, x)) {
                         using std::abs;
                         uvector<real, N> g = bernstein::evalBernsteinPolyGradient(p, x);
-                            if (type == OuterAggregate)
-                            {
+                        if (type == OuterAggregate) {
                             // When in aggregate mode, the scalar surf integral multiplies f by |n_k|^2, whose net effect
                             // is multiply weight by |n_k|; the flux surf integral multiplies f by sign(n_k) = sign(g(k))
                             real alpha = max(abs(g));
-                                if (alpha > 0)
-                                {
+                            if (alpha > 0) {
                                 g     /= alpha;
                                 alpha  = abs(g(k_active)) / norm(g);
                             }
@@ -780,14 +726,11 @@ namespace algoim
                             // pline; when near high-multiplicity roots, this simple method can break down; other, more
                             // sophisticated methods should be used in such cases, but these are not implemented here
                             f(x, w * alpha, set_component<real, N>(real(0.0), k_active, w * util::sign(g(k_active))));
-                            }
-                            else
-                            {
+                        } else {
                             // When in non-aggregate mode, the scalar surf integral multiples f by 1, whose net effect
                             // is multiply weight by 1/|n_k|; the flux surf integral multiplies f by n
                             uvector<real, N> n = g;
-                                if (norm(n) > 0)
-                                    n *= real(1.0) / norm(n);
+                            if (norm(n) > 0) n *= real(1.0) / norm(n);
                             real alpha = w * norm(g) / abs(g(k_active));
                             f(x, alpha, alpha * n);
                         }
@@ -801,10 +744,8 @@ namespace algoim
         base.integrate(strategy, q, integrand);
 
         // If in aggregate mode, apply to other dimensions as well
-            if (type == OuterAggregate)
-            {
-                for (int i = 0; i < N - 1; ++i)
-                {
+        if (type == OuterAggregate) {
+            for (int i = 0; i < N - 1; ++i) {
                 auto& [k, base] = base_other[i];
                 k_active        = k;
                 base.integrate(strategy, q, integrand);
@@ -813,7 +754,9 @@ namespace algoim
     }
 };
 
-    template<> struct ImplicitPolyQuadrature<0> {};
+template <>
+struct ImplicitPolyQuadrature<0> {
+};
 } // namespace algoim
 
 #endif

examples/examples_quad_multipoly.cpp

@@ -15,6 +15,8 @@
 #include "vector"
 #include "xarray.hpp"
 
+#include "myDebug.hpp"
+
 using namespace algoim;
 
 const static std::vector<std::vector<real>> binomial_table = {
@@ -31,19 +33,19 @@ const static std::vector<std::vector<real>> binomial_table = {
 template <int N>
 struct DebugXArray {
     template <typename Phi>
-	void operator()(const xarray<real, N>& iData, Phi&& phi) {	
+    void operator()(const xarray<real, N>& iData, Phi&& phi)
+    {
     }
 };
 
 template <>
 struct DebugXArray<2> {
     template <typename Phi>
-	void operator()(const xarray<real, 2>& iData, Phi&& phi) {	
+    void operator()(const xarray<real, 2>& iData, Phi&& phi)
+    {
         std::vector<std::vector<real>> data(iData.ext(0), std::vector<real>(iData.ext(1)));
         for (int i = 0; i < iData.ext(0); ++i) {
-				for (int j = 0; j < iData.ext(1); ++j) {
-					data[i][j] = iData(i, j);
-				}
+            for (int j = 0; j < iData.ext(1); ++j) { data[i][j] = iData(i, j); }
         }
         real              inputX1 = 0.2, inputX2 = 0.3;
         std::vector<real> b1(iData.ext(0)), b2(iData.ext(1));
@@ -57,9 +59,7 @@ struct DebugXArray<2>{
         }
         for (int i = 0; i < iData.ext(0); ++i) {
             real tmp = 0;
-				for (int j = 0; j < iData.ext(1); j++) {
-					tmp += data[i][j] * b2[j];
-				}
+            for (int j = 0; j < iData.ext(1); j++) { tmp += data[i][j] * b2[j]; }
             res += tmp * b1[i];
         }
 
@@ -74,7 +74,14 @@ struct DebugXArray<2>{
 // 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 qConv(const Phi& phi, real xmin, real xmax, uvector<int,N> P, const F& integrand, int qMax, real volume_exact, real surf_exact)
+void qConv(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);
@@ -87,19 +94,15 @@ void qConv(const Phi& phi, real xmin, real xmax, uvector<int,N> P, const F& inte
 
     // Functional to evaluate volume and surface integrals of given integrand
     real volume, surf;
-    auto compute = [&](int q)
-    {
+    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));
+        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)
-        {
+        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
@@ -108,10 +111,10 @@ void qConv(const Phi& phi, real xmin, real xmax, uvector<int,N> P, const F& inte
     };
 
     // Compute results for all q and output in a convergence table
-    for (int q = 1; q <= qMax; ++q)
-    {
+    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;
+        std::cout << q << ' ' << volume << ' ' << surf << ' ' << std::abs(volume - volume_exact) / volume_exact << ' '
+                  << std::abs(surf - surf_exact) / surf_exact << std::endl;
     }
 }
 
@@ -125,27 +128,24 @@ void outputQuadratureRuleAsVtpXML(const std::vector<uvector<real,N+1>>& q, std::
     stream << "<?xml version=\"1.0\"?>\n";
     stream << "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
     stream << "<PolyData>\n";
-    stream << "<Piece NumberOfPoints=\"" << q.size() << "\" NumberOfVerts=\"" << q.size() << "\" NumberOfLines=\"0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">\n";
+    stream << "<Piece NumberOfPoints=\"" << q.size() << "\" NumberOfVerts=\"" << q.size()
+           << "\" NumberOfLines=\"0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">\n";
     stream << "<Points>\n";
     stream << "  <DataArray type=\"Float32\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">";
-    for (const auto& pt : q)
-        stream << pt(0) << ' ' << pt(1) << ' ' << (N == 3 ? pt(2) : 0.0) << ' ';
+    for (const auto& pt : q) stream << pt(0) << ' ' << pt(1) << ' ' << (N == 3 ? pt(2) : 0.0) << ' ';
     stream << "</DataArray>\n";
     stream << "</Points>\n";
     stream << "<Verts>\n";
     stream << "  <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">";
-    for (size_t i = 0; i < q.size(); ++i)
-        stream << i << ' ';
+    for (size_t i = 0; i < q.size(); ++i) stream << i << ' ';
     stream << "</DataArray>\n";
     stream << "  <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">";
-    for (size_t i = 1; i <= q.size(); ++i)
-        stream << i << ' ';
+    for (size_t i = 1; i <= q.size(); ++i) stream << i << ' ';
     stream << "</DataArray>\n";
     stream << "</Verts>\n";
     stream << "<PointData Scalars=\"w\">\n";
     stream << "  <DataArray type=\"Float32\" Name=\"w\" NumberOfComponents=\"1\" format=\"ascii\">";
-    for (const auto& pt : q)
-        stream << pt(N) << ' ';
+    for (const auto& pt : q) stream << pt(N) << ' ';
     stream << "</DataArray>\n";
     stream << "</PointData>\n";
     stream << "</Piece>\n";
@@ -170,15 +170,13 @@ void outputQuadScheme(const F& fphi, real xmin, real xmax, const uvector<int,N>&
     // Compute quadrature scheme and record the nodes & weights; phase0 corresponds to
     // {phi < 0}, phase1 corresponds to {phi > 0}, and surf corresponds to {phi == 0}.
     std::vector<uvector<real, N + 1>> phase0, phase1, surf;
-    ipquad.integrate(AutoMixed, q, [&](const uvector<real,N>& x, real w)
-    {
+    ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
         if (bernstein::evalBernsteinPoly(phi, x) < 0)
             phase0.push_back(add_component(x, N, w));
         else
             phase1.push_back(add_component(x, N, w));
     });
-    ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real,N>& x, real w, const uvector<real,N>& wn)
-    {
+    ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, N>& x, real w, const uvector<real, N>& wn) {
         surf.push_back(add_component(x, N, w));
     });
 
@@ -192,7 +190,13 @@ void outputQuadScheme(const F& fphi, real xmin, real xmax, const uvector<int,N>&
 // 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>
-void outputQuadScheme(const F1& fphi1, const F2& fphi2, real xmin, real xmax, const uvector<int,N>& P, int q, std::string qfile)
+void outputQuadScheme(const F1&              fphi1,
+                      const F2&              fphi2,
+                      real                   xmin,
+                      real                   xmax,
+                      const uvector<int, N>& P,
+                      int                    q,
+                      std::string            qfile)
 {
     // Construct phi by mapping [0,1] onto bounding box [xmin,xmax]
     xarray<real, N> phi1(nullptr, P), phi2(nullptr, P);
@@ -207,12 +211,8 @@ void outputQuadScheme(const F1& fphi1, const F2& fphi2, real xmin, real xmax, co
     // of phi1 and phi2 in order to separate the nodes into different components, but for
     // simplicity they are agglomerated
     std::vector<uvector<real, N + 1>> vol, surf;
-    ipquad.integrate(AutoMixed, q, [&](const uvector<real,N>& x, real w)
-    {
-        vol.push_back(add_component(x, N, w));
-    });
-    ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real,N>& x, real w, const uvector<real,N>& wn)
-    {
+    ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) { vol.push_back(add_component(x, N, w)); });
+    ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, N>& x, real w, const uvector<real, N>& wn) {
         surf.push_back(add_component(x, N, w));
     });
 
@@ -221,7 +221,8 @@ void outputQuadScheme(const F1& fphi1, const F2& fphi2, real xmin, real xmax, co
     outputQuadratureRuleAsVtpXML<N>(surf, qfile + "-surf.vtp");
 }
 
-void module_test() {
+void module_test()
+{
     if (false) {
         const uvector<int, 1> P(5);
         xarray<real, 1>       beta(nullptr, P);
@@ -236,9 +237,7 @@ void module_test() {
         xarray<real, 2>       beta(nullptr, P);
         algoim_spark_alloc(real, beta);
         for (int i = 0; i < 4; i++)
-            for (int j = 0; j < 4; j++) {
-                beta(i,j)= i + j;
-            }
+            for (int j = 0; j < 4; j++) { beta(i, j) = i + j; }
         const uvector<real, 2> x(0.5, 0.3);
         auto                   res = bernstein::evalBernsteinPoly(beta, x);
         std::cout << "res: " << res << std::endl;
@@ -253,15 +252,9 @@ int main(int argc, char* argv[])
     std::cout << std::scientific << std::setprecision(10);
 
     // q-convergence study for a 2D ellipse
-    if (false) {
-        auto ellipse = [](const uvector<real,2>& x)
-        {
-            return x(0)*x(0) + x(1)*x(1)*4 - 1;
-        };
-        auto integrand = [](const uvector<real,2>& x)
-        {
-            return 1.0;
-        };
+    if (true) {
+        auto ellipse      = [](const uvector<real, 2>& x) { return x(0) * x(0) + x(1) * x(1) * 4 - 1; };
+        auto integrand    = [](const uvector<real, 2>& x) { return 1.0; };
         real volume_exact = algoim::util::pi / 2;
         real surf_exact   = 4.844224110273838099214251598195914705976959198943300412541558176231060;
         std::cout << "\n\nEllipse q-convergence test\n";
@@ -271,14 +264,8 @@ int main(int argc, char* argv[])
 
     // q-convergence study for a 3D ellipsoid
     if (false) {
-        auto ellipsoid = [](const uvector<real,3>& x)
-        {
-            return x(0)*x(0) + x(1)*x(1)*4 + x(2)*x(2)*9 - 1;
-        };
-        auto integrand = [](const uvector<real,3>& x)
-        {
-            return 1.0;
-        };
+        auto ellipsoid    = [](const uvector<real, 3>& x) { return x(0) * x(0) + x(1) * x(1) * 4 + x(2) * x(2) * 9 - 1; };
+        auto integrand    = [](const uvector<real, 3>& x) { return 1.0; };
         real volume_exact = (algoim::util::pi * 2) / 9;
         real surf_exact   = 4.400809564664970341600200389229705943483674323377145800356686868037845;
         std::cout << "\n\nEllipsoid q-convergence test\n";
@@ -289,13 +276,12 @@ int main(int argc, char* argv[])
     // Visusalisation of a 2D case involving a single polynomial; this example corresponds to
     // Figure 3, row 3, left column, https://doi.org/10.1016/j.jcp.2021.110720
     if (false) {
-        auto phi = [](const uvector<real,2>& xx)
-        {
+        auto phi = [](const uvector<real, 2>& xx) {
             real x = xx(0) * 2 - 1;
             real y = xx(1) * 2 - 1;
-            return -0.06225100787918392 + 0.1586472897571363*y + 0.5487135634635731*y*y + 
-                x*(0.3478849533965025 - 0.3321074999999999*y - 0.5595163485848738*y*y) + 
-                x*x*(0.7031095851739786 + 0.29459557349175747*y + 0.030425624999999998*y*y);
+            return -0.06225100787918392 + 0.1586472897571363 * y + 0.5487135634635731 * y * y
+                   + x * (0.3478849533965025 - 0.3321074999999999 * y - 0.5595163485848738 * y * y)
+                   + x * x * (0.7031095851739786 + 0.29459557349175747 * y + 0.030425624999999998 * y * y);
         };
         outputQuadScheme<2>(phi, 0.0, 1.0, 3, 3, "exampleA");
         std::cout << "\n\nQuadrature visualisation of a 2D case involving a single polynomial, corresponding\n";
@@ -306,20 +292,21 @@ int main(int argc, char* argv[])
     // Visusalisation of a 3D case involving a single polynomial; this example corresponds to
     // Figure 3, row 3, right column, https://doi.org/10.1016/j.jcp.2021.110720
     if (false) {
-        auto phi = [](const uvector<real,3>& xx)
-        {
+        auto phi = [](const uvector<real, 3>& xx) {
             real x = xx(0) * 2 - 1;
             real y = xx(1) * 2 - 1;
             real z = xx(2) * 2 - 1;
-            return -0.3003521613375472 - 0.22416584292513722*z + 0.07904600284034838*z*z +
-                y*(-0.022501556528537706 - 0.16299445153615613*z - 0.10968042065096766*z*z) + 
-                y*y*(0.09321375574517882 - 0.07409794846221623*z + 0.09940785133211516*z*z) + 
-                x*(0.094131400740032 - 0.11906280402685224*z - 0.010060302873268541*z*z + 
-                y*y*(0.01448948481714108 - 0.0262370580373332*z - 0.08632912757566019*z*z) + 
-                y*(0.08171132326327647 - 0.09286444275596013*z - 0.07651000354823911*z*z)) + 
-                x*x*(-0.0914370528387867 + 0.09778971384044874*z - 0.1086777644685091*z*z + 
-                y*y*(-0.04283439400630859 + 0.0750156999192893*z + 0.051754527934553866*z*z) + 
-                y*(-0.052642188754328405 - 0.03538476045586772*z + 0.11117016852276898*z*z));
+            return -0.3003521613375472 - 0.22416584292513722 * z + 0.07904600284034838 * z * z
+                   + y * (-0.022501556528537706 - 0.16299445153615613 * z - 0.10968042065096766 * z * z)
+                   + y * y * (0.09321375574517882 - 0.07409794846221623 * z + 0.09940785133211516 * z * z)
+                   + x
+                         * (0.094131400740032 - 0.11906280402685224 * z - 0.010060302873268541 * z * z
+                            + y * y * (0.01448948481714108 - 0.0262370580373332 * z - 0.08632912757566019 * z * z)
+                            + y * (0.08171132326327647 - 0.09286444275596013 * z - 0.07651000354823911 * z * z))
+                   + x * x
+                         * (-0.0914370528387867 + 0.09778971384044874 * z - 0.1086777644685091 * z * z
+                            + y * y * (-0.04283439400630859 + 0.0750156999192893 * z + 0.051754527934553866 * z * z)
+                            + y * (-0.052642188754328405 - 0.03538476045586772 * z + 0.11117016852276898 * z * z));
         };
         outputQuadScheme<3>(phi, 0.0, 1.0, 3, 3, "exampleB");
         std::cout << "\n\nQuadrature visualisation of a 3D case involving a single polynomial, corresponding\n";
@@ -330,21 +317,19 @@ int main(int argc, char* argv[])
     // 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
     if (false) {
-        auto phi0 = [](const uvector<real,2>& xx)
-        {
+        auto phi0 = [](const uvector<real, 2>& xx) {
             real x = xx(0) * 2 - 1;
             real y = xx(1) * 2 - 1;
-            return 0.014836540349115947 + 0.7022484024095262*y + 0.09974561176434385*y*y +
-                x*(0.6863910464417281 + 0.03805619999999999*y - 0.09440658332756446*y*y) + 
-                x*x*(0.19266932968830816 - 0.2325190091204104*y + 0.2957473125000001*y*y);
+            return 0.014836540349115947 + 0.7022484024095262 * y + 0.09974561176434385 * y * y
+                   + x * (0.6863910464417281 + 0.03805619999999999 * y - 0.09440658332756446 * y * y)
+                   + x * x * (0.19266932968830816 - 0.2325190091204104 * y + 0.2957473125000001 * y * y);
         };
-        auto phi1 = [](const uvector<real,2>& xx)
-        {
+        auto phi1 = [](const uvector<real, 2>& xx) {
             real x = xx(0) * 2 - 1;
             real y = xx(1) * 2 - 1;
-            return -0.18792528379702625 + 0.6713882473904913*y + 0.3778666084723582*y*y +
-                x*x*(-0.14480813208127946 + 0.0897755603159206*y - 0.141199875*y*y) + 
-                x*(-0.6169311810674598 - 0.19449299999999994*y - 0.005459163675646665*y*y);
+            return -0.18792528379702625 + 0.6713882473904913 * y + 0.3778666084723582 * y * y
+                   + x * x * (-0.14480813208127946 + 0.0897755603159206 * y - 0.141199875 * y * y)
+                   + x * (-0.6169311810674598 - 0.19449299999999994 * y - 0.005459163675646665 * y * y);
         };
         outputQuadScheme<2>(phi0, phi1, 0.0, 1.0, 3, 3, "exampleC");
         std::cout << "\n\nQuadrature visualisation of a 2D implicitly-defined domain involving the\n";
@@ -353,7 +338,8 @@ int main(int argc, char* argv[])
         std::cout << "(XML VTP file format).\n";
     }
 
-    module_test();
+    // module_test();
+    testMain();
 
     return 0;
 }

gjj/myDebug.hpp

@@ -0,0 +1,179 @@
+#include <bitset>
+#include <iostream>
+#include <booluarray.hpp>
+
+
+#include <cstddef>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <vector>
+#include "bernstein.hpp"
+#include "quadrature_multipoly.hpp"
+
+#include "real.hpp"
+#include "uvector.hpp"
+#include "vector"
+#include "xarray.hpp"
+
+
+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;
+    }
+}
+
+// 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,
+                    std::string            qfile)
+{
+    // 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
+    // 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 (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0)
+            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);
+    std::cout << "q volume: " << volume << std::endl;
+}
+
+void testMultiPloys()
+{
+    // 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
+    if (true) {
+        auto phi0 = [](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;
+        };
+        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);
+            // real y = xx(1);
+            // real z = xx(2);
+            // return x * y * z;
+        };
+        // auto phi0 = [](const uvector<real, 3>& xx) {
+        //     real x = xx(0) * 2 - 1;
+        //     real y = xx(1) * 2 - 1;
+        //     return 0.014836540349115947 + 0.7022484024095262 * y + 0.09974561176434385 * y * y
+        //            + x * (0.6863910464417281 + 0.03805619999999999 * y - 0.09440658332756446 * y * y)
+        //            + x * x * (0.19266932968830816 - 0.2325190091204104 * y + 0.2957473125000001 * y * y);
+        // };
+        // auto phi1 = [](const uvector<real, 3>& xx) {
+        //     real x = xx(0) * 2 - 1;
+        //     real y = xx(1) * 2 - 1;
+        //     return -0.18792528379702625 + 0.6713882473904913 * y + 0.3778666084723582 * y * y
+        //            + x * x * (-0.14480813208127946 + 0.0897755603159206 * y - 0.141199875 * y * y)
+        //            + x * (-0.6169311810674598 - 0.19449299999999994 * y - 0.005459163675646665 * y * y);
+        // };
+        auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
+        qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "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";
+        std::cout << "(XML VTP file format).\n";
+    }
+}
+
+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 testMain()
+{
+    testBooluarray();
+    testMultiPloys();
+}