fix calculate area bug by Chebyshev.

Author: BIMCoderLiang

src/LNLib/Algorithm/Integrator.cpp

@@ -196,6 +196,63 @@ namespace LNLib
         }
         return integration;
     }
+    double Integrator::ClenshawCurtisQuadrature2(IntegrationFunction& function, void* customData, double start, double end, std::vector<double> series, double epsilon)
+    {
+        double integration;
+        int j, k, l;
+        double err, esf, eref, erefh, hh, ir, iback, irback, ba, ss, x, y, fx, errir;
+        int lenw = series.size() - 1;
+        esf = 10;
+        ba = 0.5 * (end - start);
+        ss = 2 * series[lenw];
+        x = ba * series[lenw];
+        series[0] = 0.5 * (function)(start, customData);
+        series[3] = 0.5 * (function)(end, customData);
+        series[2] = (function)(start + x, customData);
+        series[4] = (function)(end - x, customData);
+        series[1] = (function)(start + ba, customData);
+        eref = 0.5 * (fabs(series[0]) + fabs(series[1]) + fabs(series[2]) + fabs(series[3]) + fabs(series[4]));
+        series[0] += series[3];
+        series[2] += series[4];
+        ir = series[0] + series[1] + series[2];
+        integration = series[0] * series[lenw - 1] + series[1] * series[lenw - 2] + series[2] * series[lenw - 3];
+        erefh = eref * sqrt(epsilon);
+        eref *= epsilon;
+        hh = 0.25;
+        l = 2;
+        k = lenw - 5;
+        do {
+            iback = integration;
+            irback = ir;
+            x = ba * series[k + 1];
+            y = 0;
+            integration = series[0] * series[k];
+            for (j = 1; j <= l; j++) {
+                x += y;
+                y += ss * (ba - x);
+                fx = (function)(start + x, customData) + (function)(end - x, customData);
+                ir += fx;
+                integration += series[j] * series[k - j] + fx * series[k - j - l];
+                series[j + l] = fx;
+            }
+            ss = 2 * series[k + 1];
+            err = esf * l * fabs(integration - iback);
+            hh *= 0.25;
+            errir = hh * fabs(ir - 2 * irback);
+            l *= 2;
+            k -= l + 2;
+        } while ((err > erefh || errir > eref) && k > 4 * l);
+        integration *= end - start;
+        if (err > erefh || errir > eref)
+        {
+            err *= -fabs(end - start);
+        }
+        else
+        {
+            err = eref * fabs(end - start);
+        }
+        return integration;
+    }
 }
 
 

src/LNLib/Geometry/Surface/NurbsSurface.cpp

@@ -66,7 +66,7 @@ namespace LNLib
 			AreaCoreFunction areaCoreFunction;
 			AreaData* data = (AreaData*)customData;
 			data->ParameterV = parameter;
-			return Integrator::ClenshawCurtisQuadrature(areaCoreFunction, data, data->CurrentKnotU, data->NextKnotU, series);
+			return Integrator::ClenshawCurtisQuadrature2(areaCoreFunction, data, data->CurrentKnotU, data->NextKnotU, data->Series);
 		}
 	};
 
@@ -2913,10 +2913,10 @@ double LNLib::NurbsSurface::ApproximateArea(const LN_NurbsSurface& surface, Inte
 				data.NextKnotU = knotVectorU[i + 1];
 				for (int j = degreeV; j < controlPoints[0].size(); j++) 
 				{
-					data.CurrentKnotV = knotVectorV[i];
-					data.NextKnotV = knotVectorV[i + 1];
+					data.CurrentKnotV = knotVectorV[j];
+					data.NextKnotV = knotVectorV[j + 1];
 
-					area += Integrator::ClenshawCurtisQuadrature(function, (void*)&data, data.CurrentKnotV, data.NextKnotV, series);
+					area += Integrator::ClenshawCurtisQuadrature2(function, (void*)&data, data.CurrentKnotV, data.NextKnotV, series);
 				}
 			}
 			break;

src/LNLib/include/Integrator.h

@@ -39,7 +39,8 @@ namespace LNLib
 		/// According to https://github.com/chrisidefix/nurbs
 		/// </summary>
 		static std::vector<double> ChebyshevSeries(int size = 100);
-		static double ClenshawCurtisQuadrature(IntegrationFunction& function, void* customData, double start, double end, std::vector<double>& series, double epsilon = Constants::DistanceEpsilon);	
+		static double ClenshawCurtisQuadrature(IntegrationFunction& function, void* customData, double start, double end, std::vector<double>& series, double epsilon = Constants::DistanceEpsilon);
+		static double ClenshawCurtisQuadrature2(IntegrationFunction& function, void* customData, double start, double end, std::vector<double> series, double epsilon = Constants::DistanceEpsilon);
 	};
 }
 

tests/T_Additional.cpp

@@ -119,8 +119,8 @@ TEST(Test_Additional, Area)
 
 	double simpson = NurbsSurface::ApproximateArea(surface, IntegratorType::Simpson);
 	double gaussLegendre = NurbsSurface::ApproximateArea(surface, IntegratorType::GaussLegendre);
+	double chebyshev = NurbsSurface::ApproximateArea(surface, IntegratorType::Chebyshev);
 	EXPECT_FALSE(MathUtils::IsAlmostEqualTo(simpson, standardArea)); // not accuracy when use Simpson
 	EXPECT_TRUE(MathUtils::IsAlmostEqualTo(gaussLegendre, standardArea));
-	//Need fix, to be continued...
-	//double chebyshev = NurbsSurface::ApproximateArea(surface, IntegratorType::Chebyshev);
+	EXPECT_TRUE(MathUtils::IsAlmostEqualTo(chebyshev, standardArea));
 }