| 
						
						
							
								
							
						
						
					 | 
					@ -36,6 +36,15 @@ | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					#include "tanhsinh.hpp" | 
					 | 
					 | 
					#include "tanhsinh.hpp" | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					#include "bernstein.hpp" | 
					 | 
					 | 
					#include "bernstein.hpp" | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					
 | 
					 | 
					 | 
					
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#define STOP_WHEN_BLOCKED true | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#if STOP_WHEN_BLOCKED | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#include <chrono> | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					auto         timerStart            = std::chrono::high_resolution_clock::now(); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					bool         stopCurrentQuadrature = false; | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					const double MAX_DURATION          = 10; // seconds
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#endif | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					// #define
 | 
					 | 
					 | 
					// #define
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					
 | 
					 | 
					 | 
					
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					namespace algoim | 
					 | 
					 | 
					namespace algoim | 
				
			
			
		
	
	
		
		
			
				
					| 
						
							
								
							
						
						
							
								
							
						
						
					 | 
					@ -404,6 +413,14 @@ void eliminate_axis(PolySet<N, ALGOIM_M>& phi, int k, PolySet<N - 1, ALGOIM_M>& | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    // R
 | 
					 | 
					 | 
					    // R
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    for (int i = 0; i < phi.count(); ++i) | 
					 | 
					 | 
					    for (int i = 0; i < phi.count(); ++i) | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        for (int j = i + 1; j < phi.count(); ++j) { | 
					 | 
					 | 
					        for (int j = i + 1; j < phi.count(); ++j) { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#if STOP_WHEN_BLOCKED | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					            auto                          timerEnd = std::chrono::high_resolution_clock::now(); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					            std::chrono::duration<double> duration = timerEnd - timerStart; | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					            if (duration.count() > MAX_DURATION) { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					                stopCurrentQuadrature = true; | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					                return; | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					            } // 超过20s
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#endif | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            const auto& p     = phi.poly(i); | 
					 | 
					 | 
					            const auto& p     = phi.poly(i); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            const auto& pmask = phi.mask(i); | 
					 | 
					 | 
					            const auto& pmask = phi.mask(i); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            const auto& q     = phi.poly(j); | 
					 | 
					 | 
					            const auto& q     = phi.poly(j); | 
				
			
			
		
	
	
		
		
			
				
					| 
						
							
								
							
						
						
							
								
							
						
						
					 | 
					@ -551,6 +568,9 @@ struct ImplicitPolyQuadrature { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    // Build quadrature hierarchy for a domain implicitly defined by a single polynomial
 | 
					 | 
					 | 
					    // Build quadrature hierarchy for a domain implicitly defined by a single polynomial
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    ImplicitPolyQuadrature(const xarray<real, N>& p) | 
					 | 
					 | 
					    ImplicitPolyQuadrature(const xarray<real, N>& p) | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    { | 
					 | 
					 | 
					    { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#if STOP_WHEN_BLOCKED | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					        timerStart = std::chrono::high_resolution_clock::now(); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#endif | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true)); | 
					 | 
					 | 
					        auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true)); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					
 | 
					 | 
					 | 
					
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        if (!detail::maskEmpty(mask)) { | 
					 | 
					 | 
					        if (!detail::maskEmpty(mask)) { | 
				
			
			
		
	
	
		
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
					@ -563,6 +583,9 @@ struct ImplicitPolyQuadrature { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    // Build quadrature hierarchy for a domain implicitly defined by two polynomials
 | 
					 | 
					 | 
					    // Build quadrature hierarchy for a domain implicitly defined by two polynomials
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    ImplicitPolyQuadrature(const xarray<real, N>& p, const xarray<real, N>& q) | 
					 | 
					 | 
					    ImplicitPolyQuadrature(const xarray<real, N>& p, const xarray<real, N>& q) | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    { | 
					 | 
					 | 
					    { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#if STOP_WHEN_BLOCKED | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					        timerStart = std::chrono::high_resolution_clock::now(); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#endif | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        { | 
					 | 
					 | 
					        { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            // auto tmp  = booluarray<N, ALGOIM_M>(false);
 | 
					 | 
					 | 
					            // auto tmp  = booluarray<N, ALGOIM_M>(false);
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true)); | 
					 | 
					 | 
					            auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true)); | 
				
			
			
		
	
	
		
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
					@ -580,6 +603,9 @@ struct ImplicitPolyQuadrature { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    // Build quadrature hierarchy for a domain implicitly defined by multiple polynomials
 | 
					 | 
					 | 
					    // Build quadrature hierarchy for a domain implicitly defined by multiple polynomials
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    ImplicitPolyQuadrature(const std::vector<xarray<real, N>>& ps) | 
					 | 
					 | 
					    ImplicitPolyQuadrature(const std::vector<xarray<real, N>>& ps) | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    { | 
					 | 
					 | 
					    { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#if STOP_WHEN_BLOCKED | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					        timerStart = std::chrono::high_resolution_clock::now(); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#endif | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        for (const auto& p : ps) { | 
					 | 
					 | 
					        for (const auto& p : ps) { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true)); | 
					 | 
					 | 
					            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); | 
				
			
			
		
	
	
		
		
			
				
					| 
						
							
								
							
						
						
							
								
							
						
						
					 | 
					@ -608,6 +634,9 @@ struct ImplicitPolyQuadrature { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					
 | 
					 | 
					 | 
					
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    ImplicitPolyQuadrature(const std::vector<tensor3>& tensors, const std::vector<int>& polyIndices) | 
					 | 
					 | 
					    ImplicitPolyQuadrature(const std::vector<tensor3>& tensors, const std::vector<int>& polyIndices) | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    { | 
					 | 
					 | 
					    { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#if STOP_WHEN_BLOCKED | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					        timerStart = std::chrono::high_resolution_clock::now(); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#endif | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        for (auto i : polyIndices) { | 
					 | 
					 | 
					        for (auto i : polyIndices) { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            auto mask = detail::nonzeroMask(tensors[i], booluarray<N, ALGOIM_M>(true)); | 
					 | 
					 | 
					            auto mask = detail::nonzeroMask(tensors[i], booluarray<N, ALGOIM_M>(true)); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            if (!detail::maskEmpty(mask)) phi.push_back(tensors[i], mask); | 
					 | 
					 | 
					            if (!detail::maskEmpty(mask)) phi.push_back(tensors[i], mask); | 
				
			
			
		
	
	
		
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
					@ -617,6 +646,9 @@ struct ImplicitPolyQuadrature { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					
 | 
					 | 
					 | 
					
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    ImplicitPolyQuadrature(const std::vector<organizer::MinimalPrimitiveRep>& minimalReps, const std::vector<int>& polyIndices) | 
					 | 
					 | 
					    ImplicitPolyQuadrature(const std::vector<organizer::MinimalPrimitiveRep>& minimalReps, const std::vector<int>& polyIndices) | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    { | 
					 | 
					 | 
					    { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#if STOP_WHEN_BLOCKED | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					        timerStart = std::chrono::high_resolution_clock::now(); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#endif | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        for (auto i : polyIndices) { | 
					 | 
					 | 
					        for (auto i : polyIndices) { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            auto mask = detail::nonzeroMask(minimalReps[i].tensor, booluarray<N, ALGOIM_M>(true)); | 
					 | 
					 | 
					            auto mask = detail::nonzeroMask(minimalReps[i].tensor, booluarray<N, ALGOIM_M>(true)); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            if (!detail::maskEmpty(mask)) phi.push_back(minimalReps[i].tensor, mask); | 
					 | 
					 | 
					            if (!detail::maskEmpty(mask)) phi.push_back(minimalReps[i].tensor, mask); | 
				
			
			
		
	
	
		
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
					@ -630,6 +662,9 @@ struct ImplicitPolyQuadrature { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					                           const xarray<real, N>&         q, | 
					 | 
					 | 
					                           const xarray<real, N>&         q, | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					                           const booluarray<N, ALGOIM_M>& qmask) | 
					 | 
					 | 
					                           const booluarray<N, ALGOIM_M>& qmask) | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    { | 
					 | 
					 | 
					    { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#if STOP_WHEN_BLOCKED | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					        timerStart = std::chrono::high_resolution_clock::now(); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					#endif | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        { | 
					 | 
					 | 
					        { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            auto mask = detail::nonzeroMask(p, pmask); | 
					 | 
					 | 
					            auto mask = detail::nonzeroMask(p, pmask); | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					            if (!maskEmpty(mask)) phi.push_back(p, mask); | 
					 | 
					 | 
					            if (!maskEmpty(mask)) phi.push_back(p, mask); | 
				
			
			
		
	
	
		
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
					@ -645,6 +680,7 @@ struct ImplicitPolyQuadrature { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    // Assuming phi has been instantiated, determine elimination axis and build base
 | 
					 | 
					 | 
					    // Assuming phi has been instantiated, determine elimination axis and build base
 | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    void build(bool outer, bool auto_apply_TS) | 
					 | 
					 | 
					    void build(bool outer, bool auto_apply_TS) | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					    { | 
					 | 
					 | 
					    { | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					 | 
					 | 
					 | 
					        if (stopCurrentQuadrature) { return; } | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        type                = outer ? OuterSingle : Inner; | 
					 | 
					 | 
					        type                = outer ? OuterSingle : Inner; | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					        this->auto_apply_TS = auto_apply_TS; | 
					 | 
					 | 
					        this->auto_apply_TS = auto_apply_TS; | 
				
			
			
		
	
		
		
			
				
					 | 
					 | 
					
 | 
					 | 
					 | 
					
 | 
				
			
			
		
	
	
		
		
			
				
					| 
						
							
								
							
						
						
						
					 | 
					
  |