diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..52ae38d --- /dev/null +++ b/.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: ^ + 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 diff --git a/CMakeLists.txt b/CMakeLists.txt index a775ff1..bf73fd9 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/algoim/bernstein.hpp b/algoim/bernstein.hpp index 91def97..cde10c9 100644 --- a/algoim/bernstein.hpp +++ b/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); } } diff --git a/algoim/booluarray.hpp b/algoim/booluarray.hpp index b02795c..81d6648 100644 --- a/algoim/booluarray.hpp +++ b/algoim/booluarray.hpp @@ -8,56 +8,45 @@ namespace algoim { - namespace booluarray_detail +namespace booluarray_detail +{ +constexpr int pow(int base, int exp) { return exp == 0 ? 1 : base * pow(base, exp - 1); } + +template +constexpr int furl(const uvector& i) // 为什么这里可以是constexpr，uvector不是变量吗？ +{ + int ind = i(0); + 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 +// specialisation of uarray +template +class booluarray +{ + constexpr static int size = booluarray_detail::pow(E, N); + + +public: + std::bitset bits; + + booluarray() {} + + booluarray(bool val) { - constexpr int pow(int base, int exp) - { - return exp == 0 ? 1 : base * pow(base, exp - 1); - } - - template - constexpr int furl(const uvector& i) - { - int ind = i(0); - for (int j = 1; j < N; ++j) - ind = ind * E + i(j); - return ind; - } + if (val) bits.set(); // 全部置为1 } - // booluarray implements a simple N-dimensional array of booleans, with - // compile-time extent E across all dimensions; it is essentially a basic - // specialisation of uarray - template - class booluarray - { - constexpr static int size = booluarray_detail::pow(E, N); - std::bitset bits; - public: - booluarray() {} - - booluarray(bool val) - { - if (val) - bits.set(); - } - - bool operator() (const uvector& i) const - { - return bits[booluarray_detail::furl(i)]; - } - - auto operator() (const uvector& i) - { - return bits[booluarray_detail::furl(i)]; - } - - // returns true iff the entire array is false - bool none() const - { - return bits.none(); - } - }; + bool operator()(const uvector& i) const { return bits[booluarray_detail::furl(i)]; } + + auto operator()(const uvector& i) { return bits[booluarray_detail::furl(i)]; } + + // returns true iff the entire array is false + bool none() const { return bits.none(); } +}; } // namespace algoim #endif diff --git a/algoim/polyset.hpp b/algoim/polyset.hpp index 08c8fd6..f95cebb 100644 --- a/algoim/polyset.hpp +++ b/algoim/polyset.hpp @@ -6,50 +6,48 @@ #include #include #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 - struct PolySet - { - struct Poly - { - uvector ext; // Degree/extent of polynomial - size_t offset; // Offset into buffer, storing the xarray polynomial data - booluarray mask; // Mask - }; - std::vector buff; // Memory buffer containing polynomial data - std::vector items; // Record of contained polynomials - - // Access polynomial by index - xarray poly(size_t ind) - { - assert(0 <= ind && ind < items.size()); - return xarray(&buff[items[ind].offset], items[ind].ext); - } - - // Access mask by index - booluarray& mask(size_t ind) - { - assert(0 <= ind && ind < items.size()); - return items[ind].mask; - } - - // Add a polynomial/mask pair to the container - void push_back(const xarray& p, const booluarray& m) - { - items.push_back({p.ext(), buff.size(), m}); - buff.resize(buff.size() + p.size()); - poly(items.size() - 1) = p; - } - - size_t count() const - { - return items.size(); - } +// PolySet implements a simple container to hold one or more Bernstein polynomials +// and their associated masks +template +struct PolySet { + struct Poly { + uvector ext; // Degree/extent of polynomial + size_t offset; // Offset into buffer, storing the xarray polynomial data + booluarray mask; // Mask }; + + std::vector buff; // Memory buffer containing polynomial data + std::vector items; // Record of contained polynomials + + // Access polynomial by index + xarray poly(size_t ind) + { + assert(0 <= ind && ind < items.size()); + return xarray(&buff[items[ind].offset], items[ind].ext); + } + + // Access mask by index + booluarray& mask(size_t ind) + { + assert(0 <= ind && ind < items.size()); + return items[ind].mask; + } + + // Add a polynomial/mask pair to the container + void push_back(const xarray& p, const booluarray& m) + { + items.push_back({p.ext(), buff.size(), m}); + buff.resize(buff.size() + p.size()); + poly(items.size() - 1) = p; + } + + size_t count() const { return items.size(); } +}; } // namespace algoim #endif diff --git a/algoim/quadrature_multipoly.hpp b/algoim/quadrature_multipoly.hpp index ec23ff0..071dbb7 100644 --- a/algoim/quadrature_multipoly.hpp +++ b/algoim/quadrature_multipoly.hpp @@ -3,7 +3,7 @@ // High-order accurate quadrature algorithms for multi-component domains implicitly // defined by (one or more) multivariate Bernstein polynomials, based on the algorithms -// developed in the paper +// developed in the paper // R. I. Saye, High-order quadrature on multi-component domains implicitly defined // by multivariate polynomials, Journal of Computational Physics, 448, 110720 (2022), // https://doi.org/10.1016/j.jcp.2021.110720 @@ -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 #define ALGOIM_M 8 #include @@ -32,788 +33,730 @@ namespace algoim { - namespace detail - { - template - booluarray mask_driver(const xarray& f, const booluarray& fmask, const xarray* g, const booluarray* gmask) - { - booluarray mask(false); - auto helper = [&](auto&& self, uvector a, uvector b) - { - bool overlap = false; - for (MultiLoop i(a,b); ~i; ++i) - if (fmask(i()) && (!gmask || (*gmask)(i()))) - overlap = true; - if (!overlap) - return; - - real eps = 0.015625 / ALGOIM_M; - uvector xa, xb; - 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) - { - xarray 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 - { - xarray fab(nullptr, f.ext()); - algoim_spark_alloc(real, fab); - bernstein::deCasteljau(f, xa, xb, fab); - if (bernstein::uniformSign(fab) != 0) - return; - } - - if (b(0) - a(0) == 1) - { - assert(all(b - a == 1)); - assert(fmask(a) && (!gmask || (*gmask)(a))); - mask(a) = true; - return; - } - - assert(all(b - a > 1) && all((b - a) % 2 == 0)); - uvector delta = (b - a) / 2; - for (MultiLoop i(0,2); ~i; ++i) - self(self, a + i() * delta, a + (i() + 1) * delta); - }; - helper(helper, 0, ALGOIM_M); - return mask; +namespace detail +{ +template +booluarray mask_driver(const xarray& f, + const booluarray& fmask, + const xarray* g, + const booluarray* gmask) +{ + booluarray mask(false); + auto helper = [&](auto&& self, uvector a, uvector b) { + bool overlap = false; + for (MultiLoop i(a, b); ~i; ++i) + if (fmask(i()) && (!gmask || (*gmask)(i()))) overlap = true; + if (!overlap) return; + + real eps = 0.015625 / ALGOIM_M; + uvector xa, xb; + for (int dim = 0; dim < N; ++dim) { + xa(dim) = real(a(dim)) / ALGOIM_M - eps; + xb(dim) = real(b(dim)) / ALGOIM_M + eps; } - // Using orthant tests, compute a mask for the possible subrectangles for which f and g share - // common zeros, i.e., intersecting zero level sets; if the mask is false somewhere, then it is - // guaranteed that f and/or g were originally masked off at the same place, or that they - // 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 - booluarray intersectionMask(const xarray& f, const booluarray& fmask, const xarray& g, const booluarray& gmask) - { - return mask_driver(f, fmask, &g, &gmask); + if (g) { + xarray 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 { + xarray fab(nullptr, f.ext()); + algoim_spark_alloc(real, fab); + bernstein::deCasteljau(f, xa, xb, fab); + if (bernstein::uniformSign(fab) != 0) return; } - // Using orthant tests, compute a mask for the possible subrectangles for which f has a zero - // level set; if the mask is false somewhere, then it is guaranteed that f was originally - // masked off at the same place, or that it definitively does not have any zeros in that - // subrectangle; if the returned mask is true, then shared zeros may exist (and with high likelihood) - template - booluarray nonzeroMask(const xarray& f, const booluarray& fmask) - { - return mask_driver(f, fmask, nullptr, nullptr); + if (b(0) - a(0) == 1) { + assert(all(b - a == 1)); + assert(fmask(a) && (!gmask || (*gmask)(a))); + mask(a) = true; + return; } - // Collapse a mask along dimension k by bitwise-OR-ing along columns - template - booluarray collapseMask(const booluarray& mask, int k) - { - booluarray r(false); - for (MultiLoop i(0,ALGOIM_M); ~i; ++i) - if (mask(i())) - r(remove_component(i(), k)) = true; - return r; - } + assert(all(b - a > 1) && all((b - a) % 2 == 0)); + uvector delta = (b - a) / 2; + for (MultiLoop i(0, 2); ~i; ++i) self(self, a + i() * delta, a + (i() + 1) * delta); + }; + helper(helper, 0, ALGOIM_M); + return mask; +} + +// Using orthant tests, compute a mask for the possible subrectangles for which f and g share +// common zeros, i.e., intersecting zero level sets; if the mask is false somewhere, then it is +// guaranteed that f and/or g were originally masked off at the same place, or that they +// 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 +booluarray intersectionMask(const xarray& f, + const booluarray& fmask, + const xarray& g, + const booluarray& gmask) +{ + return mask_driver(f, fmask, &g, &gmask); +} + +// Using orthant tests, compute a mask for the possible subrectangles for which f has a zero +// level set; if the mask is false somewhere, then it is guaranteed that f was originally +// masked off at the same place, or that it definitively does not have any zeros in that +// subrectangle; if the returned mask is true, then shared zeros may exist (and with high likelihood) +template +booluarray nonzeroMask(const xarray& f, const booluarray& fmask) +{ + return mask_driver(f, fmask, nullptr, nullptr); +} - // Test if a mask is empty, i.e., all entries are false - template - bool maskEmpty(const booluarray& mask) - { - return mask.none(); - } +// Collapse a mask along dimension k by bitwise-OR-ing along columns +template +booluarray collapseMask(const booluarray& mask, int k) +{ + booluarray r(false); + for (MultiLoop i(0, ALGOIM_M); ~i; ++i) + if (mask(i())) r(remove_component(i(), k)) = true; + return r; +} + +// Test if a mask is empty, i.e., all entries are false +template +bool maskEmpty(const booluarray& mask) +{ + return mask.none(); +} - // Test if a point x \in [0,1]^N is in a true subrectangle of a mask; if x is exactly - // on the border between two subrectangles, the left subrectangle shall be used - template - bool pointWithinMask(const booluarray& mask, const uvector& x) - { - using std::floor; - uvector cell; - for (int dim = 0; dim < N; ++dim) +// Test if a point x \in [0,1]^N is in a true subrectangle of a mask; if x is exactly +// on the border between two subrectangles, the left subrectangle shall be used +template +bool pointWithinMask(const booluarray& mask, const uvector& x) +{ + using std::floor; + uvector cell; + for (int dim = 0; dim < N; ++dim) + cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast(floor(x(dim) * ALGOIM_M)))); + return mask(cell); +} + +// Test if a point {x + alpha e_k} is in a true subrectangle of a mask for some alpha \in [0,1] +template +bool lineIntersectsMask(const booluarray& mask, const uvector& x, int k) +{ + using std::floor; + if constexpr (N > 1) { + uvector cell; + for (int dim = 0; dim < N; ++dim) + if (dim < k) cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast(floor(x(dim) * ALGOIM_M)))); - return mask(cell); + else if (dim > k) + cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast(floor(x(dim - 1) * ALGOIM_M)))); + for (int i = 0; i < ALGOIM_M; ++i) { + cell(k) = i; + if (mask(cell)) return true; } + return false; + } else + return !maskEmpty(mask); +} - // Test if a point {x + alpha e_k} is in a true subrectangle of a mask for some alpha \in [0,1] - template - bool lineIntersectsMask(const booluarray& mask, const uvector& x, int k) - { - using std::floor; - if constexpr (N > 1) - { - uvector cell; - for (int dim = 0; dim < N; ++dim) - if (dim < k) - cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast(floor(x(dim) * ALGOIM_M)))); - else if (dim > k) - cell(dim) = std::max(0, std::min(ALGOIM_M - 1, static_cast(floor(x(dim - 1) * ALGOIM_M)))); - for (int i = 0; i < ALGOIM_M; ++i) - { - cell(k) = i; - if (mask(cell)) - return true; - } - return false; +template +void restrictToFace(const xarray& a, int k, int side, xarray& out) +{ + 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) { + uvector j; + 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); + } +} + +template +booluarray restrictToFace(const booluarray& a, int k, int side) +{ + assert(0 <= k && k < N && (side == 0 || side == 1)); + booluarray r; + for (MultiLoop i(0, ALGOIM_M); ~i; ++i) { + uvector j; + 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; +} + +// Compute determinant of the given matrix using QR + Givens rotations + column +// pivoting, along with approximated rank +// in: square matrix A, which will be overwritten +template +T det_qr(xarray& A, int& rank, T tol = 10.0) +{ + assert(A.ext(0) == A.ext(1) && A.ext(0) > 0); + 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) { + T m = -1; + int k = -1; + 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) { + m = mag; + k = i; } - else - return !maskEmpty(mask); } - - template - void restrictToFace(const xarray& a, int k, int side, xarray& out) - { - 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) - { - uvector j; - 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); - } + assert(j <= k && k < n); + if (k != j) { + for (int a = 0; a < n; ++a) std::swap(A(a, j), A(a, k)); + det *= -1.0; } - - template - booluarray restrictToFace(const booluarray& a, int k, int side) - { - assert(0 <= k && k < N && (side == 0 || side == 1)); - booluarray r; - for (MultiLoop i(0,ALGOIM_M); ~i; ++i) - { - uvector j; - 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; + 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); } + det *= A(j, j); + max_diag_r = max(max_diag_r, abs(A(j, j))); + } - // Compute determinant of the given matrix using QR + Givens rotations + column - // pivoting, along with approximated rank - // in: square matrix A, which will be overwritten - template - T det_qr(xarray& A, int& rank, T tol = 10.0) - { - assert(A.ext(0) == A.ext(1) && A.ext(0) > 0); - using std::max; - using std::abs; - T det = 1.0; - int n = A.ext(0); - T max_diag_r = 0.0; - for (int j = 0; j < n; ++j) - { - T m = -1; - int k = -1; - 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) - { - 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)); - det *= -1.0; - } - 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); - } - det *= A(j,j); - max_diag_r = max(max_diag_r, abs(A(j,j))); - } + tol *= max_diag_r * n * std::numeric_limits::epsilon(); + rank = 0; + for (int i = 0; i < n; ++i) + if (abs(A(i, i)) > tol) ++rank; - tol *= max_diag_r * n * std::numeric_limits::epsilon(); - rank = 0; - for (int i = 0; i < n; ++i) - if (abs(A(i,i)) > tol) - ++rank; + return det; +} - return det; - } +// Determine the largest possible degree of the resultant of two general polynomials +template +uvector resultantExtent(const uvector& p, const uvector& q, int dim) +{ + uvector ext; + 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; + } + return ext; +} + +// Determine the largest possible degree of the discriminant of a polynomial +template +uvector discriminantExtent(const uvector& p, int dim) +{ + uvector ext; + 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; + } + return ext; +} + +// Compute the resultant of p and q and store the result in out +// ========================================= NOTE ========================================= +// This is a heavily simplified method and does not handle rank deficiency caused by, e.g., +// common polynomial factors, nor does it handle ill-conditioning caused by extreme values, +// among various other aspects. If your application requires this kind of special handling, +// consider contacting the author of this code for suggestions. +// ========================================= NOTE ========================================= +template +bool resultant_core(const xarray& p, const xarray* q, int k, xarray& out) +{ + assert(0 <= k && k < N); - // Determine the largest possible degree of the resultant of two general polynomials - template - uvector resultantExtent(const uvector& p, const uvector& q, int dim) - { - uvector ext; - 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; - } - return ext; - } + int P = p.ext(k); + int Q = q ? q->ext(k) : P - 1; + int M = (P == Q) ? P - 1 : P + Q - 2; + xarray mat(nullptr, uvector{M, M}); - // Determine the largest possible degree of the discriminant of a polynomial - template - uvector discriminantExtent(const uvector& p, int dim) - { - uvector ext; - 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; - } - return ext; - } + assert(P >= 2 && Q >= 1 && M >= 1); - // Compute the resultant of p and q and store the result in out - // ========================================= NOTE ========================================= - // This is a heavily simplified method and does not handle rank deficiency caused by, e.g., - // common polynomial factors, nor does it handle ill-conditioning caused by extreme values, - // among various other aspects. If your application requires this kind of special handling, - // consider contacting the author of this code for suggestions. - // ========================================= NOTE ========================================= - template - bool resultant_core(const xarray& p, const xarray* q, int k, xarray& out) - { - assert(0 <= k && k < N); - - int P = p.ext(k); - int Q = q ? q->ext(k) : P - 1; - int M = (P == Q) ? P - 1 : P + Q - 2; - assert(P >= 2 && Q >= 1 && M >= 1); - - xarray f(nullptr, out.ext()); - xarray mat(nullptr, uvector{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) - { - uvector x; - 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) - bernstein::collapseAlongAxis(*q, x, k, qk); - else - bernstein::bernsteinDerivative(pk, P, qk); + xarray f(nullptr, out.ext()); - if (P == Q) - bernstein::bezoutMatrix(pk, qk, P, mat); - else - bernstein::sylvesterMatrix(pk, P, qk, Q, mat); + real *pk, *qk; + algoim_spark_alloc(real, f, mat); + algoim_spark_alloc(real, &pk, P, &qk, Q); - int rank; - f.l(i) = det_qr(mat, rank); - } + for (auto i = f.loop(); ~i; ++i) { + uvector x; + for (int dim = 0; dim < N - 1; ++dim) x(dim) = bernstein::modifiedChebyshevNode(i(dim), f.ext(dim)); - // Interpolate the resultant on the tensor-product grid - bernstein::normalise(f); - bernstein::bernsteinInterpolate(f, std::pow(100.0 * std::numeric_limits::epsilon(), 1.0 / (N - 1)), out); + bernstein::collapseAlongAxis(p, x, k, pk); + if (q) + bernstein::collapseAlongAxis(*q, x, k, qk); + else + bernstein::bernsteinDerivative(pk, P, qk); - // Try for polynomial degree reduction - bool b = bernstein::autoReduction(out, 1e4 * std::numeric_limits::epsilon()); + if (P == Q) + bernstein::bezoutMatrix(pk, qk, P, mat); + else + bernstein::sylvesterMatrix(pk, P, qk, Q, mat); - // 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); - return true; - } + int rank; + f.l(i) = det_qr(mat, rank); + } - // Compute the pseudo-resultant R(p,q) along dimension k - template - bool resultant(const xarray& p, const xarray& q, int k, xarray& out) - { - return resultant_core(p, &q, k, out); - } + // Interpolate the resultant on the tensor-product grid + bernstein::normalise(f); + bernstein::bernsteinInterpolate(f, std::pow(100.0 * std::numeric_limits::epsilon(), 1.0 / (N - 1)), out); - // Compute the (intentially unnormalised) pseudo-discriminant R(p,p') along dimension k - template - bool discriminant(const xarray& p, int k, xarray& out) - { - xarray prime(nullptr, inc_component(p.ext(), k, -1)); - algoim_spark_alloc(real, prime); - bernstein::bernsteinDerivative(p, k, prime); - return resultant_core(p, &prime, k, out); - } - - // Using the polynomials in phi, eliminate the axis k by restricting to faces, computing - // discriminants and resultants, and storing the computed polynomials in psi - template - void eliminate_axis(PolySet& phi, int k, PolySet& psi) - { - static_assert(N >= 2, "N >= 2 required to eliminate axis"); - assert(0 <= k && k < N); - assert(psi.count() == 0); + // Try for polynomial degree reduction + bool b = bernstein::autoReduction(out, 1e4 * std::numeric_limits::epsilon()); - // For every phi(i) ... - for (int i = 0; i < phi.count(); ++i) - { - const auto& p = phi.poly(i); - const auto& mask = phi.mask(i); + // 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); + return true; +} - // Examine bottom and top faces in the k'th dimension - for (int side = 0; side <= 1; ++side) - { - xarray 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)) - { - bernstein::autoReduction(p_face); - bernstein::normalise(p_face); - psi.push_back(p_face, p_face_mask); - } - } +// Compute the pseudo-resultant R(p,q) along dimension k +template +bool resultant(const xarray& p, const xarray& q, int k, xarray& out) +{ + return resultant_core(p, &q, k, out); +} - // Consider discriminant - xarray p_k(nullptr, p.ext()); - 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)) - { - // note: computed disc might have lower degree than the following - uvector R = discriminantExtent(p.ext(), k); - xarray disc(nullptr, R); - algoim_spark_alloc(real, disc); - if (discriminant(p, k, disc)) - { - bernstein::normalise(disc); - psi.push_back(disc, collapseMask(disc_mask, k)); - } - } +// Compute the (intentially unnormalised) pseudo-discriminant R(p,p') along dimension k +template +bool discriminant(const xarray& p, int k, xarray& out) +{ + xarray prime(nullptr, inc_component(p.ext(), k, -1)); + algoim_spark_alloc(real, prime); + bernstein::bernsteinDerivative(p, k, prime); + return resultant_core(p, &prime, k, out); +} + +// Using the polynomials in phi, eliminate the axis k by restricting to faces, computing +// discriminants and resultants, and storing the computed polynomials in psi +template +void eliminate_axis(PolySet& phi, int k, PolySet& psi) +{ + static_assert(N >= 2, "N >= 2 required to eliminate axis"); + assert(0 <= k && k < N); + assert(psi.count() == 0); + + // For every phi(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) { + xarray 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)) { + bernstein::autoReduction(p_face); + bernstein::normalise(p_face); + psi.push_back(p_face, p_face_mask); } - - // Consider every pairwise combination of resultants ... - 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)) - { - // note: computed resultant might have lower degree than the following - uvector R = resultantExtent(p.ext(), q.ext(), k); - xarray res(nullptr, R); - algoim_spark_alloc(real, res); - if (resultant(p, q, k, res)) - { - bernstein::normalise(res); - psi.push_back(res, collapseMask(mask, k)); - } - } - }; } - // Compute the 'score' across all dimensions - template - uvector score_estimate(PolySet& phi, uvector& has_disc) - { - static_assert(N > 1, "score_estimate of practical use only with N > 1"); - using std::abs; - uvector s = 0; - has_disc = false; - // For every phi(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 j(0,ALGOIM_M); ~j; ++j) if (mask(j())) - { - uvector x = (j() + 0.5)/real(ALGOIM_M); - uvector g = bernstein::evalBernsteinPolyGradient(p, x); - real sum = 0; - for (int dim = 0; dim < N; ++dim) - { - g(dim) = abs(g(dim)); - sum += g(dim); - } - if (sum > 0) - s += g / sum; + // Consider discriminant + xarray p_k(nullptr, p.ext()); + 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)) { + // note: computed disc might have lower degree than the following + uvector R = discriminantExtent(p.ext(), k); + xarray disc(nullptr, R); + algoim_spark_alloc(real, disc); + if (discriminant(p, k, disc)) { + bernstein::normalise(disc); + psi.push_back(disc, collapseMask(disc_mask, k)); + } + } + } + + // Consider every pairwise combination of resultants ... + 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)) { + // note: computed resultant might have lower degree than the following + uvector R = resultantExtent(p.ext(), q.ext(), k); + xarray res(nullptr, R); + algoim_spark_alloc(real, res); + if (resultant(p, q, k, res)) { + bernstein::normalise(res); + psi.push_back(res, collapseMask(mask, k)); } + } + }; +} - // Consider discriminant - xarray p_k(nullptr, p.ext()); - algoim_spark_alloc(real, p_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); +// Compute the 'score' across all dimensions +template +uvector score_estimate(PolySet& phi, uvector& has_disc) +{ + static_assert(N > 1, "score_estimate of practical use only with N > 1"); + using std::abs; + uvector s = 0; + has_disc = false; + // For every phi(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 j(0, ALGOIM_M); ~j; ++j) + if (mask(j())) { + uvector x = (j() + 0.5) / real(ALGOIM_M); + uvector g = bernstein::evalBernsteinPolyGradient(p, x); + real sum = 0; + for (int dim = 0; dim < N; ++dim) { + g(dim) = abs(g(dim)); + sum += g(dim); } + if (sum > 0) s += g / sum; } - return s; + + // Consider discriminant + xarray p_k(nullptr, p.ext()); + algoim_spark_alloc(real, p_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); } - } // namespace detail - - // Main engine for generating high-order accurature quadrature schemes on multi-component domains - // implicitly defined by multivariate Bernstein polynomials in the unit hyperrectangle [0,1]^N. - // See examples_quad_multipoly.cpp for examples demonstrating the usage of these methods. - // See also additional examples provided on the GitHub documentation page, - // https://algoim.github.io/ - - // After building the quadrature hierarchy (via dimension reduction), the specific kind - // of quadrature scheme applied to the one-dimensional line integrals is chosen by the - // user via a QuadStrategy parameter: - // AlwaysGL: applies Gauss-Legendre quadrature on every interval, regardless of - // the level; this strategy is generally good for small-to-medium q - // AlwaysTS: applies tanh-sinh quadrature on every interval, regardless of the level; - // this strategy is generally the worst of all three and is provided mainly - // for exploratory purposes - // AutoMixed: apply the automated strategy discussed in detail in the paper - // https://doi.org/10.1016/j.jcp.2021.110720; briefly, the method: (i) applies - // tanh-sinh on base integrals for which vertical tangents are detected in - // the height-function-based representation of their implicitly-defined - // geometry; (ii) applies Gauss-Legendre quadrature on the inner-most - // integral, and on all base integrals whose geometry is the graph of a - // multi-valued height function devoid of vertical tangents/branching. - enum QuadStrategy { AlwaysGL, AlwaysTS, AutoMixed }; - - template - struct ImplicitPolyQuadrature + } + return s; +} +} // namespace detail + +// Main engine for generating high-order accurature quadrature schemes on multi-component domains +// implicitly defined by multivariate Bernstein polynomials in the unit hyperrectangle [0,1]^N. +// See examples_quad_multipoly.cpp for examples demonstrating the usage of these methods. +// See also additional examples provided on the GitHub documentation page, +// https://algoim.github.io/ + +// After building the quadrature hierarchy (via dimension reduction), the specific kind +// of quadrature scheme applied to the one-dimensional line integrals is chosen by the +// user via a QuadStrategy parameter: +// AlwaysGL: applies Gauss-Legendre quadrature on every interval, regardless of +// the level; this strategy is generally good for small-to-medium q +// AlwaysTS: applies tanh-sinh quadrature on every interval, regardless of the level; +// this strategy is generally the worst of all three and is provided mainly +// for exploratory purposes +// AutoMixed: apply the automated strategy discussed in detail in the paper +// https://doi.org/10.1016/j.jcp.2021.110720; briefly, the method: (i) applies +// tanh-sinh on base integrals for which vertical tangents are detected in +// the height-function-based representation of their implicitly-defined +// geometry; (ii) applies Gauss-Legendre quadrature on the inner-most +// integral, and on all base integrals whose geometry is the graph of a +// multi-valued height function devoid of vertical tangents/branching. +enum QuadStrategy { AlwaysGL, AlwaysTS, AutoMixed }; + +template +struct ImplicitPolyQuadrature { + enum IntegralType { Inner, OuterSingle, OuterAggregate }; + + PolySet phi; // Given N-dimensional polynomials + int k; // Elimination axis/height direction; k = N if there are no interfaces + // 这种struct的递归定义一般是不允许的，因为会内存爆炸，但是模板类可以通过特化定义递归的终止条件 + // 所以有template<> struct ImplicitPolyQuadrature<0> {}; + ImplicitPolyQuadrature 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 ? tanh-sinh? + IntegralType type; // Whether an inner integral, or outer of two kinds + std::array>, N - 1> + base_other; // Stores other base cases, besides k, when in aggregate mode + + // Default ctor sets to an uninitialised state + ImplicitPolyQuadrature() : k(-1) {} + + // Build quadrature hierarchy for a domain implicitly defined by a single polynomial + ImplicitPolyQuadrature(const xarray& p) + { + auto mask = detail::nonzeroMask(p, booluarray(true)); + if (!detail::maskEmpty(mask)) phi.push_back(p, mask); + build(true, false); + } + + // Build quadrature hierarchy for a domain implicitly defined by two polynomials + ImplicitPolyQuadrature(const xarray& p, const xarray& q) { - enum IntegralType { Inner, OuterSingle, OuterAggregate }; - - PolySet phi; // Given N-dimensional polynomials - int k; // Elimination axis/height direction; k = N if there are no interfaces - // 这种struct的递归定义一般是不允许的，因为会内存爆炸，但是模板类可以通过特化定义递归的终止条件 - // 所以有template<> struct ImplicitPolyQuadrature<0> {}; - ImplicitPolyQuadrature 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 ? - IntegralType type; // Whether an inner integral, or outer of two kinds - std::array>,N-1> base_other; // Stores other base cases, besides k, when in aggregate mode - - // Default ctor sets to an uninitialised state - ImplicitPolyQuadrature() : k(-1) {} - - // Build quadrature hierarchy for a domain implicitly defined by a single polynomial - ImplicitPolyQuadrature(const xarray& p) { - auto mask = detail::nonzeroMask(p, booluarray(true)); - if (!detail::maskEmpty(mask)) - phi.push_back(p, mask); - build(true, false); + // auto tmp = booluarray(false); + auto mask = detail::nonzeroMask(p, booluarray(true)); + // tmp(0) = true; + // auto res = !detail::maskEmpty(tmp); // 似乎是测试代码？ + if (!detail::maskEmpty(mask)) phi.push_back(p, mask); } - - // Build quadrature hierarchy for a domain implicitly defined by two polynomials - ImplicitPolyQuadrature(const xarray& p, const xarray& q) { - { - auto tmp = booluarray(false); - auto mask = detail::nonzeroMask(p, booluarray(true)); - tmp(0) = true; - auto res = !detail::maskEmpty(tmp); - if (!detail::maskEmpty(mask)) - phi.push_back(p, mask); - } - { - auto mask = detail::nonzeroMask(q, booluarray(true)); - if (!detail::maskEmpty(mask)) - phi.push_back(q, mask); - } - build(true, false); + auto mask = detail::nonzeroMask(q, booluarray(true)); + if (!detail::maskEmpty(mask)) phi.push_back(q, mask); } + build(true, false); + } - // Build quadrature hierarchy for a given domain implicitly defined by two polynomials with user-defined masks - ImplicitPolyQuadrature(const xarray& p, const booluarray& pmask, const xarray& q, const booluarray& qmask) + // Build quadrature hierarchy for a domain implicitly defined by multiple polynomials + ImplicitPolyQuadrature(const std::vector>& ps) + { + for (const auto& p : ps) { + auto mask = detail::nonzeroMask(p, booluarray(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& p, + const booluarray& pmask, + const xarray& q, + const booluarray& qmask) + { { - { - auto mask = detail::nonzeroMask(p, pmask); - if (!maskEmpty(mask)) - phi.push_back(p, mask); - } - { - auto mask = detail::nonzeroMask(q, qmask); - if (!maskEmpty(mask)) - phi.push_back(q, mask); - } - build(true, false); + auto mask = detail::nonzeroMask(p, pmask); + if (!maskEmpty(mask)) phi.push_back(p, mask); } - - // Assuming phi has been instantiated, determine elimination axis and build base - void build(bool outer, bool auto_apply_TS) { - type = outer ? OuterSingle : Inner; - this->auto_apply_TS = auto_apply_TS; - - // If phi is empty, apply a tensor-product Gaussian quadrature - if (phi.count() == 0) - { - k = N; - this->auto_apply_TS = false; - return; - } + auto mask = detail::nonzeroMask(q, qmask); + if (!maskEmpty(mask)) phi.push_back(q, mask); + } + build(true, false); + } - 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 - { - // Compute score; penalise any directions which likely contain vertical tangents - uvector has_disc; - uvector score = detail::score_estimate(phi, has_disc); - //if (max(abs(score)) == 0) - // score(0) = 1.0; - assert(max(abs(score)) > 0); - score /= 2 * max(abs(score)); + // Assuming phi has been instantiated, determine elimination axis and build base + void build(bool outer, bool auto_apply_TS) + { + type = outer ? OuterSingle : Inner; + this->auto_apply_TS = auto_apply_TS; + + // If phi is empty, apply a tensor-product Gaussian quadrature + if (phi.count() == 0) { + k = N; + this->auto_apply_TS = false; + return; + } + + 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 { + // Compute score; penalise any directions which likely contain vertical tangents + uvector has_disc; + uvector score = detail::score_estimate(phi, has_disc); + // if (max(abs(score)) == 0) + // score(0) = 1.0; + 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; + + // 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 + // discriminant mask has been found + k = argmax(score); + detail::eliminate_axis(phi, k, base.phi); + base.build(false, this->auto_apply_TS || has_disc(k)); + + // If this is the outer integral, and surface quadrature schemes are required, apply + // the dimension-aggregated scheme when necessary + if (outer && has_disc(k)) { + type = OuterAggregate; for (int i = 0; i < N; ++i) - 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 - // discriminant mask has been found - k = argmax(score); - detail::eliminate_axis(phi, k, base.phi); - base.build(false, this->auto_apply_TS || has_disc(k)); - - // If this is the outer integral, and surface quadrature schemes are required, apply - // the dimension-aggregated scheme when necessary - if (outer && has_disc(k)) - { - type = OuterAggregate; - for (int i = 0; i < N; ++i) if (i != k) - { + if (i != k) { auto& [kother, base] = base_other[i < k ? i : i - 1]; - kother = i; + kother = i; detail::eliminate_axis(phi, kother, base.phi); // In aggregate mode, triggered by non-empty discriminant mask, // base integrals always have T-S suggested base.build(false, this->auto_apply_TS || true); } - } } } + } - // Integrate a functional via quadrature of the base integral, adding the dimension k - template - void integrate(QuadStrategy strategy, int q, const F& f) - { - assert(0 <= k && k <= N); - - // If there are no interfaces, apply rectangular tensor-product Gauss-Legendre quadrature - // 所以这个是积分域包裹了整个box时的做法？ - if (k == N) - { - assert(!auto_apply_TS); - for (MultiLoop i(0, q); ~i; ++i) - { - uvector x; - real w = 1.0; - for (int dim = 0; dim < N; ++dim) - { - x(dim) = GaussQuad::x(q, i(dim)); // 0 <= i(dim) < q - w *= GaussQuad::w(q, i(dim)); - } - f(x, w); + // Integrate a functional via quadrature of the base integral, adding the dimension k + template + void integrate(QuadStrategy strategy, int q, const F& f) + { + assert(0 <= k && k <= N); + + // If there are no interfaces, apply rectangular tensor-product Gauss-Legendre quadrature + // 所以这个是积分域包裹了整个box时的做法？ + if (k == N) { + assert(!auto_apply_TS); + for (MultiLoop i(0, q); ~i; ++i) { + uvector x; + real w = 1.0; + for (int dim = 0; dim < N; ++dim) { + x(dim) = GaussQuad::x(q, i(dim)); // 0 <= i(dim) < q + w *= GaussQuad::w(q, i(dim)); } - return; + f(x, w); } + return; + } - // 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; - - // Base integral invokes the following integrand - auto integrand = [&](const uvector& xbase, real w) - { - // Allocate node buffer of sufficient size and initialise with {0, 1} - real *nodes; - algoim_spark_alloc(real, &nodes, max_count); - nodes[0] = 0.0; - nodes[1] = 1.0; - int count = 2; - - // For every phi(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; - - // Restrict polynomial to axis-aligned line and compute its roots - real *pline, *roots; - algoim_spark_alloc(real, &pline, P, &roots, P - 1); - bernstein::collapseAlongAxis(p, xbase, k, pline); - 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) - { - auto x = add_component(xbase, k, roots[j]); - if (detail::pointWithinMask(mask, x)) - nodes[count++] = roots[j]; - } - }; - - // Sort the nodes - std::sort(nodes, nodes + count); - assert(nodes[0] == real(0) && nodes[count-1] == real(1)); - - // Force nearly-degenerate sub-intervals to be exactly degenerate - real tol = 10.0 * std::numeric_limits::epsilon(); - using std::abs; - for (int i = 1; i < count - 1; ++i) - if (abs(nodes[i]) < tol) - nodes[i] = 0.0; - else if (abs(nodes[i] - 1) < tol) - nodes[i] = 1.0; - else if (abs(nodes[i] - nodes[i+1]) < tol) - nodes[i+1] = nodes[i]; - 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) - { - real x0 = nodes[i]; - real x1 = nodes[i + 1]; - 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 (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)); + // 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; + + // Base integral invokes the following integrand + auto integrand = [&](const uvector& xbase, real w) { + // Allocate node buffer of sufficient size and initialise with {0, 1} + real* nodes; + algoim_spark_alloc(real, &nodes, max_count); + nodes[0] = 0.0; + nodes[1] = 1.0; + int count = 2; + + // For every phi(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; + + // Restrict polynomial to axis-aligned line and compute its roots + real *pline, *roots; + algoim_spark_alloc(real, &pline, P, &roots, P - 1); + bernstein::collapseAlongAxis(p, xbase, k, pline); + 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) { + auto x = add_component(xbase, k, roots[j]); + if (detail::pointWithinMask(mask, x)) nodes[count++] = roots[j]; } }; - // When N = 1, the base case of recursion is invoked on a zero-dimensional point with unit weight - if constexpr (N > 1) - base.integrate(strategy, q, integrand); - else - integrand(uvector(), real(1)); - } + // Sort the nodes + std::sort(nodes, nodes + count); + assert(nodes[0] == real(0) && nodes[count - 1] == real(1)); - // Surface-integrate a functional via quadrature of the base integral, adding the dimension k - template - void integrate_surf(QuadStrategy strategy, int q, const F& f) - { - static_assert(N > 1, "surface integral only implemented in N > 1 dimensions"); - assert(type == OuterSingle || type == OuterAggregate); - - // If there is no interface, there is no surface integral - 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& xbase, real w) - { - assert(0 <= k_active && k_active < N); - // For every phi(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; - - // Compute roots of { x \mapsto phi(xbase + x e_k) } - real *pline, *roots; - algoim_spark_alloc(real, &pline, P, &roots, P - 1); - bernstein::collapseAlongAxis(p, xbase, k_active, pline); - int rcount = bernstein::bernsteinUnitIntervalRealRoots(pline, P, roots); - - // 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) - { - auto x = add_component(xbase, k_active, roots[j]); - if (detail::pointWithinMask(mask, x)) - { - using std::abs; - uvector g = bernstein::evalBernsteinPolyGradient(p, x); - 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) - { - g /= alpha; - alpha = abs(g(k_active)) / norm(g); - } - // Simplistic method to compute sign(n_k). NOTE: This method relies on a reasonable consistency - // between the gradient calculation of original polynomial, and that of the roots computed from - // 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(0.0), k_active, w * util::sign(g(k_active)))); - } - 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 n = g; - if (norm(n) > 0) - n *= real(1.0) / norm(n); - real alpha = w * norm(g) / abs(g(k_active)); - f(x, alpha, alpha * n); + // Force nearly-degenerate sub-intervals to be exactly degenerate + real tol = 10.0 * std::numeric_limits::epsilon(); + using std::abs; + for (int i = 1; i < count - 1; ++i) + if (abs(nodes[i]) < tol) + nodes[i] = 0.0; + else if (abs(nodes[i] - 1) < tol) + nodes[i] = 1.0; + else if (abs(nodes[i] - nodes[i + 1]) < tol) + nodes[i + 1] = nodes[i]; + 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) { + real x0 = nodes[i]; + real x1 = nodes[i + 1]; + 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 (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)); + } + }; + + // When N = 1, the base case of recursion is invoked on a zero-dimensional point with unit weight + if constexpr (N > 1) + base.integrate(strategy, q, integrand); + else + integrand(uvector(), real(1)); + } + + // Surface-integrate a functional via quadrature of the base integral, adding the dimension k + template + void integrate_surf(QuadStrategy strategy, int q, const F& f) + { + static_assert(N > 1, "surface integral only implemented in N > 1 dimensions"); + assert(type == OuterSingle || type == OuterAggregate); + + // If there is no interface, there is no surface integral + 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& xbase, real w) { + assert(0 <= k_active && k_active < N); + // For every phi(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; + + // Compute roots of { x \mapsto phi(xbase + x e_k) } + real *pline, *roots; + algoim_spark_alloc(real, &pline, P, &roots, P - 1); + bernstein::collapseAlongAxis(p, xbase, k_active, pline); + int rcount = bernstein::bernsteinUnitIntervalRealRoots(pline, P, roots); + + // 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) { + auto x = add_component(xbase, k_active, roots[j]); + if (detail::pointWithinMask(mask, x)) { + using std::abs; + uvector g = bernstein::evalBernsteinPolyGradient(p, x); + 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) { + g /= alpha; + alpha = abs(g(k_active)) / norm(g); } + // Simplistic method to compute sign(n_k). NOTE: This method relies on a reasonable consistency + // between the gradient calculation of original polynomial, and that of the roots computed from + // 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(0.0), k_active, w * util::sign(g(k_active)))); + } 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 n = g; + if (norm(n) > 0) n *= real(1.0) / norm(n); + real alpha = w * norm(g) / abs(g(k_active)); + f(x, alpha, alpha * n); } } - }; + } }; + }; - // Apply primary base integral - k_active = k; - base.integrate(strategy, q, integrand); + // Apply primary base integral + k_active = k; + 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) - { - auto& [k, base] = base_other[i]; - k_active = k; - 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) { + auto& [k, base] = base_other[i]; + k_active = k; + base.integrate(strategy, q, integrand); } } - }; + } +}; - template<> struct ImplicitPolyQuadrature<0> {}; +template <> +struct ImplicitPolyQuadrature<0> { +}; } // namespace algoim #endif diff --git a/gjj/myDebugTool.cpp b/examples/examples_mytest.cpp similarity index 100% rename from gjj/myDebugTool.cpp rename to examples/examples_mytest.cpp diff --git a/examples/examples_quad_multipoly.cpp b/examples/examples_quad_multipoly.cpp index 5ce4124..3cbcbe2 100644 --- a/examples/examples_quad_multipoly.cpp +++ b/examples/examples_quad_multipoly.cpp @@ -1,6 +1,6 @@ // Examples to demonstrate Algoim's methods for computing high-order accurate quadrature schemes // on multi-component domains implicitly-defined by (one or more) multivariate Bernstein -// polynomials. Additional examples are provided on the GitHub documentation page, +// polynomials. Additional examples are provided on the GitHub documentation page, // https://algoim.github.io/ #include @@ -15,137 +15,137 @@ #include "vector" #include "xarray.hpp" +#include "myDebug.hpp" + using namespace algoim; const static std::vector> binomial_table = { - {1, 0, 0, 0, 0, 0, 0}, - {1, 1, 0, 0, 0, 0, 0}, - {1, 2, 1, 0, 0, 0, 0}, - {1, 3, 3, 1, 0, 0, 0}, - {1, 4, 6, 4, 1, 0, 0}, - {1, 5, 10, 10, 5, 1, 0}, + {1, 0, 0, 0, 0, 0, 0}, + {1, 1, 0, 0, 0, 0, 0}, + {1, 2, 1, 0, 0, 0, 0}, + {1, 3, 3, 1, 0, 0, 0}, + {1, 4, 6, 4, 1, 0, 0}, + {1, 5, 10, 10, 5, 1, 0}, {1, 6, 15, 20, 15, 6, 1} }; // function模板不允许直接<部分>特化，所以用类模板 -template +template struct DebugXArray { - template - void operator()(const xarray& iData, Phi&& phi) { - } + template + void operator()(const xarray& iData, Phi&& phi) + { + } }; -template<> -struct DebugXArray<2>{ - template - void operator()(const xarray& iData, Phi&& phi) { - std::vector> data(iData.ext(0), std::vector(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); - } - } - real inputX1 = 0.2, inputX2 = 0.3; - std::vector b1(iData.ext(0)), b2(iData.ext(1)); - int n1 = iData.ext(0) - 1, n2 = iData.ext(1) - 1; - real res = 0; - for (int i = 0; i < iData.ext(0); ++i) { - b1[i] = binomial_table[n1][i] * std::pow(inputX1, i) * std::pow(1 - inputX1, n1-i); - } - for (int i = 0; i < iData.ext(1); ++i) { - b2[i] = binomial_table[n2][i] * std::pow(inputX2, i) * std::pow(1 - inputX2, n2-i); - } - 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]; - } - res += tmp * b1[i]; - } +template <> +struct DebugXArray<2> { + template + void operator()(const xarray& iData, Phi&& phi) + { + std::vector> data(iData.ext(0), std::vector(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); } + } + real inputX1 = 0.2, inputX2 = 0.3; + std::vector b1(iData.ext(0)), b2(iData.ext(1)); + int n1 = iData.ext(0) - 1, n2 = iData.ext(1) - 1; + real res = 0; + for (int i = 0; i < iData.ext(0); ++i) { + b1[i] = binomial_table[n1][i] * std::pow(inputX1, i) * std::pow(1 - inputX1, n1 - i); + } + for (int i = 0; i < iData.ext(1); ++i) { + b2[i] = binomial_table[n2][i] * std::pow(inputX2, i) * std::pow(1 - inputX2, n2 - i); + } + 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]; } + res += tmp * b1[i]; + } - // original phi function evaluation - const uvector x(inputX1, inputX2); - real phiEval = phi(x); - } + // original phi function evaluation + const uvector x(inputX1, inputX2); + real phiEval = phi(x); + } }; // 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 -void qConv(const Phi& phi, real xmin, real xmax, uvector P, const F& integrand, int qMax, real volume_exact, real surf_exact) +template +void qConv(const Phi& phi, + real xmin, + real xmax, + uvector 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 phipoly(nullptr, P); + xarray phipoly(nullptr, P); algoim_spark_alloc(real, phipoly); - bernstein::bernsteinInterpolate([&](const uvector& x) { return phi(xmin + x * (xmax - xmin)); }, phipoly); - DebugXArray()(phipoly, [&](const uvector& x) { return phi(xmin + x * (xmax - xmin)); }); + bernstein::bernsteinInterpolate([&](const uvector& x) { return phi(xmin + x * (xmax - xmin)); }, phipoly); + DebugXArray()(phipoly, [&](const uvector& x) { return phi(xmin + x * (xmax - xmin)); }); // Build quadrature hierarchy ImplicitPolyQuadrature ipquad(phipoly); // 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; + surf = 0.0; // compute volume integral over {phi < 0} using AutoMixed strategy - ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) - { - if (bernstein::evalBernsteinPoly(phipoly, x) < 0) - volume += w * integrand(xmin + x * (xmax - xmin)); + ipquad.integrate(AutoMixed, q, [&](const uvector& 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& x, real w, const uvector& wn) - { + ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& wn) { surf += w * integrand(xmin + x * (xmax - xmin)); }); // scale appropriately volume *= pow(xmax - xmin, N); - surf *= pow(xmax - xmin, N - 1); + surf *= pow(xmax - xmin, N - 1); }; // 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; } } // Given a set of quadrature points and weights, output them to an VTP XML file for visualisation // purposes, e.g., using ParaView -template -void outputQuadratureRuleAsVtpXML(const std::vector>& q, std::string fn) +template +void outputQuadratureRuleAsVtpXML(const std::vector>& q, std::string fn) { static_assert(N == 2 || N == 3, "outputQuadratureRuleAsVtpXML only supports 2D and 3D quadrature schemes"); std::ofstream stream(fn); stream << "\n"; stream << "\n"; stream << "\n"; - stream << "\n"; + stream << "\n"; stream << "\n"; stream << " "; - 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 << "\n"; stream << "\n"; stream << "\n"; stream << " "; - for (size_t i = 0; i < q.size(); ++i) - stream << i << ' '; - stream << "\n"; + for (size_t i = 0; i < q.size(); ++i) stream << i << ' '; + stream << "\n"; stream << " "; - for (size_t i = 1; i <= q.size(); ++i) - stream << i << ' '; + for (size_t i = 1; i <= q.size(); ++i) stream << i << ' '; stream << "\n"; stream << "\n"; stream << "\n"; stream << " "; - for (const auto& pt : q) - stream << pt(N) << ' '; + for (const auto& pt : q) stream << pt(N) << ' '; stream << "\n"; stream << "\n"; stream << "\n"; @@ -156,29 +156,27 @@ void outputQuadratureRuleAsVtpXML(const std::vector>& q, std:: // Driver method which takes a functor phi defining a single polynomial in the reference // rectangle [xmin, xmax]^N, 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 -void outputQuadScheme(const F& fphi, real xmin, real xmax, const uvector& P, int q, std::string qfile) +template +void outputQuadScheme(const F& fphi, real xmin, real xmax, const uvector& P, int q, std::string qfile) { // Construct phi by mapping [0,1] onto bounding box [xmin,xmax] - xarray phi(nullptr, P); + xarray phi(nullptr, P); algoim_spark_alloc(real, phi); - bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi(xmin + x * (xmax - xmin)); }, phi); + bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi(xmin + x * (xmax - xmin)); }, phi); // Build quadrature hierarchy ImplicitPolyQuadrature ipquad(phi); // 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> phase0, phase1, surf; - ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) - { + std::vector> phase0, phase1, surf; + ipquad.integrate(AutoMixed, q, [&](const uvector& 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& x, real w, const uvector& wn) - { + ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& wn) { surf.push_back(add_component(x, N, w)); }); @@ -191,14 +189,20 @@ void outputQuadScheme(const F& fphi, real xmin, real xmax, const uvector& // 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 -void outputQuadScheme(const F1& fphi1, const F2& fphi2, real xmin, real xmax, const uvector& P, int q, std::string qfile) +template +void outputQuadScheme(const F1& fphi1, + const F2& fphi2, + real xmin, + real xmax, + const uvector& P, + int q, + std::string qfile) { // Construct phi by mapping [0,1] onto bounding box [xmin,xmax] - xarray phi1(nullptr, P), phi2(nullptr, P); + xarray phi1(nullptr, P), phi2(nullptr, P); algoim_spark_alloc(real, phi1, phi2); - bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi1(xmin + x * (xmax - xmin)); }, phi1); - bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2); + bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi1(xmin + x * (xmax - xmin)); }, phi1); + bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2); // Build quadrature hierarchy ImplicitPolyQuadrature ipquad(phi1, phi2); @@ -206,41 +210,36 @@ void outputQuadScheme(const F1& fphi1, const F2& fphi2, real xmin, real xmax, co // Compute quadrature scheme and record the nodes & weights; one could examine the signs // of phi1 and phi2 in order to separate the nodes into different components, but for // simplicity they are agglomerated - std::vector> vol, surf; - ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) - { - vol.push_back(add_component(x, N, w)); - }); - ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& wn) - { + std::vector> vol, surf; + ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { vol.push_back(add_component(x, N, w)); }); + ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& wn) { surf.push_back(add_component(x, N, w)); }); // output to a file outputQuadratureRuleAsVtpXML(vol, qfile + "-vol.vtp"); - outputQuadratureRuleAsVtpXML(surf, qfile + "-surf.vtp"); + outputQuadratureRuleAsVtpXML(surf, qfile + "-surf.vtp"); } -void module_test() { - if (false){ +void module_test() +{ + if (false) { const uvector P(5); - xarray beta(nullptr, P); + xarray beta(nullptr, P); algoim_spark_alloc(real, beta); - for(int i = 0; i < 5; i++) beta[i]=i+1; + for (int i = 0; i < 5; i++) beta[i] = i + 1; const uvector x(0.5); - auto res = bernstein::evalBernsteinPoly(beta, x); + auto res = bernstein::evalBernsteinPoly(beta, x); std::cout << "res: " << res << std::endl; } if (true) { const uvector P(4, 4); - xarray beta(nullptr, P); + xarray 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 i = 0; i < 4; i++) + for (int j = 0; j < 4; j++) { beta(i, j) = i + j; } const uvector x(0.5, 0.3); - auto res = bernstein::evalBernsteinPoly(beta, x); + auto res = bernstein::evalBernsteinPoly(beta, x); std::cout << "res: " << res << std::endl; } } @@ -253,34 +252,22 @@ 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& x) - { - return x(0)*x(0) + x(1)*x(1)*4 - 1; - }; - auto integrand = [](const uvector& x) - { - return 1.0; - }; + if (true) { + auto ellipse = [](const uvector& x) { return x(0) * x(0) + x(1) * x(1) * 4 - 1; }; + auto integrand = [](const uvector& x) { return 1.0; }; real volume_exact = algoim::util::pi / 2; - real surf_exact = 4.844224110273838099214251598195914705976959198943300412541558176231060; + real surf_exact = 4.844224110273838099214251598195914705976959198943300412541558176231060; std::cout << "\n\nEllipse q-convergence test\n"; - std::cout << "q area(q) perim(q) area error perim error\n"; // perimeter: 周长 + std::cout << "q area(q) perim(q) area error perim error\n"; // perimeter: 周长 qConv<2>(ellipse, -1.1, 1.1, 3, integrand, 50, volume_exact, surf_exact); } // q-convergence study for a 3D ellipsoid if (false) { - auto ellipsoid = [](const uvector& x) - { - return x(0)*x(0) + x(1)*x(1)*4 + x(2)*x(2)*9 - 1; - }; - auto integrand = [](const uvector& x) - { - return 1.0; - }; + auto ellipsoid = [](const uvector& x) { return x(0) * x(0) + x(1) * x(1) * 4 + x(2) * x(2) * 9 - 1; }; + auto integrand = [](const uvector& x) { return 1.0; }; real volume_exact = (algoim::util::pi * 2) / 9; - real surf_exact = 4.400809564664970341600200389229705943483674323377145800356686868037845; + real surf_exact = 4.400809564664970341600200389229705943483674323377145800356686868037845; std::cout << "\n\nEllipsoid q-convergence test\n"; std::cout << "q volume(q) surf(q) vol error surf error\n"; qConv<3>(ellipsoid, -1.1, 1.1, 3, integrand, 50, volume_exact, surf_exact); @@ -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& 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); + auto phi = [](const uvector& 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); }; 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"; @@ -305,21 +291,22 @@ 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& 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)); + if (false) { + auto phi = [](const uvector& 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)); }; 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"; @@ -329,22 +316,20 @@ 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& 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); + if (false) { + auto phi0 = [](const uvector& 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& 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 phi1 = [](const uvector& 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); }; 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; } diff --git a/gjj/myDebug.hpp b/gjj/myDebug.hpp new file mode 100644 index 0000000..817eddd --- /dev/null +++ b/gjj/myDebug.hpp @@ -0,0 +1,179 @@ +#include +#include +#include + + +#include +#include +#include +#include +#include +#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 +void qConv1(const Phi& phi, + real xmin, + real xmax, + uvector 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 phipoly(nullptr, P); + algoim_spark_alloc(real, phipoly); + bernstein::bernsteinInterpolate([&](const uvector& x) { return phi(xmin + x * (xmax - xmin)); }, phipoly); + + // Build quadrature hierarchy + ImplicitPolyQuadrature 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& 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& x, real w, const uvector& wn) { + surf += w * integrand(xmin + x * (xmax - xmin)); + }); + // scale appropriately + volume *= pow(xmax - xmin, N); + surf *= pow(xmax - xmin, N - 1); + }; + + // Compute 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 +void qConvMultiPloy(const F1& fphi1, + const F2& fphi2, + real xmin, + real xmax, + const uvector& P, + const F& integrand, + int q, + std::string qfile) +{ + // Construct phi by mapping [0,1] onto bounding box [xmin,xmax] + xarray phi1(nullptr, P), phi2(nullptr, P); + algoim_spark_alloc(real, phi1, phi2); + bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi1(xmin + x * (xmax - xmin)); }, phi1); + // bernstein::bernsteinInterpolate([&](const uvector& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2); + + // Build quadrature hierarchy + // ImplicitPolyQuadrature ipquad(phi1, phi2); + ImplicitPolyQuadrature ipquad(phi1); + + // Functional to evaluate volume and surface integrals of given integrand + real volume, surf; + auto compute = [&](int q) { + volume = 0.0; + surf = 0.0; + // compute volume integral over {phi < 0} using AutoMixed strategy + ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { + // if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0) + if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin)); + }); + // compute surface integral over {phi == 0} using AutoMixed strategy + ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& wn) { + 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& 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& 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& 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& 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& 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(); +} diff --git a/gjj/myDebugTool.hpp b/gjj/myDebugTool.hpp new file mode 100644 index 0000000..e69de29