diff --git a/source/module_base/math_bspline.cpp b/source/module_base/math_bspline.cpp index 378fcd1ed5..63d95e3f97 100644 --- a/source/module_base/math_bspline.cpp +++ b/source/module_base/math_bspline.cpp @@ -32,18 +32,12 @@ namespace ModuleBase } } - void Bspline::cleanp() - { - delete[] bezier; - bezier = NULL; - } - double Bspline::bezier_ele(int n) { return this->bezier[n]; } - void Bspline::getbslpine(double x) + void Bspline::getbspline(double x) { bezier[0] = 1.0; for(int k = 1 ; k <= norder ; ++k) diff --git a/source/module_base/math_bspline.h b/source/module_base/math_bspline.h index 2eb72d8880..86155bb79d 100644 --- a/source/module_base/math_bspline.h +++ b/source/module_base/math_bspline.h @@ -4,57 +4,54 @@ namespace ModuleBase { - -// -//DESCRIPTION: -// A class to treat Cardinal B-spline interpolation. -// qianrui created 2021-09-14 -//MATH: -// Only uniform nodes are considered: xm-x[m-1]=Dx(>= 0) for control node: X={x0,x1,...,xm}; -// Any function p(x) can be written by -// p(x)=\sum_i{ci*M_ik(x)} (k->infinity), -// where M_ik is the i-th k-order Cardinal B-spline base function -// and ci is undetermined coefficient. -// M_i0 = H(x-xi)-H(x-x[i+1]), H(x): step function -// x-xi x[i+k+1]-x -// M_ik(x)= ---------*M_i(k-1)(x)+ ----------------*M_[i+1][k-1](x) ( xi <= x <= x[i+1] ) -// x[i+k]-xi x[i+k+1]-x[i+1] -// For uniform nodes: M_[i+1]k(x+Dx)=M_ik(x) -// If we define Bk[n] stores M_ik(x+n*Dx) for x in (xi,xi+Dx): -// x+n*Dx-xi xi+(k-n+1)*Dx-x -// Bk[n] = -----------*B(k-1)[n] + -----------------*B(k-1)[n-1] -// k*Dx k*Dx -//USAGE： -// ModuleBase::Bspline bp; -// bp.init(10,0.7,2); //Dx = 0.7, xi = 2 -// bp.getbslpine(0.5); //x = 0.5 -// cout<= 0) for control node: X={x0,x1,...,xm}; + * Any function p(x) can be written by + * p(x)=\sum_i{ci*M_ik(x)} (k->infinity), + * where M_ik is the i-th k-order Cardinal B-spline base function + * and ci is undetermined coefficient. + * M_i0 = H(x-xi)-H(x-x[i+1]), H(x): step function + * x-xi x[i+k+1]-x + * M_ik(x)= ---------*M_i(k-1)(x)+ ----------------*M_[i+1][k-1](x) ( xi <= x <= x[i+1] ) + * x[i+k]-xi x[i+k+1]-x[i+1] + * For uniform nodes: M_[i+1]k(x+Dx)=M_ik(x) + * If we define Bk[n] stores M_ik(x+n*Dx) for x in (xi,xi+Dx): + * x+n*Dx-xi xi+(k-n+1)*Dx-x + * Bk[n] = -----------*B(k-1)[n] + -----------------*B(k-1)[n-1] + * k*Dx k*Dx + * USAGE： + * ModuleBase::Bspline bp; + * bp.init(10,0.7,2); //Dx = 0.7, xi = 2 + * bp.getbslpine(0.5); //x = 0.5 + * cout<= 0 - double Dx; //Dx: the interval of control node - double xi; // xi: the starting point - double *bezier; //bezier[n] = Bk[n] + private: + int norder; // the order of bezier base; norder >= 0 + double Dx; // Dx: the interval of control node + double xi; // xi: the starting point + double *bezier; // bezier[n] = Bk[n] - public: - Bspline(); - ~Bspline(); + public: + Bspline(); + ~Bspline(); - //Init norder, Dx, xi - void init(int norderin, double Dxin, double xiin); + void init(int norderin, double Dxin, double xiin); - //delete[] bezier - void cleanp(); + // Get the result of i-th bezier base functions for different input x+xi+n*Dx. + // x should be in [0,Dx] + // n-th result is stored in bezier[n]; + void getbspline(double x); - //Get the result of i-th bezier base functions for different input x+xi+n*Dx. - //x should be in [0,Dx] - //n-th result is stored in bezier[n]; - void getbslpine(double x); - - //get the element of bezier - double bezier_ele(int n); + // get the element of bezier + double bezier_ele(int n); }; -} +} // namespace ModuleBase #endif diff --git a/source/module_base/test/CMakeLists.txt b/source/module_base/test/CMakeLists.txt index 7fa8a6c5e2..46359f2b5e 100644 --- a/source/module_base/test/CMakeLists.txt +++ b/source/module_base/test/CMakeLists.txt @@ -82,3 +82,7 @@ AddTest( LIBS ${math_libs} SOURCES math_ylmreal_test.cpp ../math_ylmreal.cpp ../realarray.cpp ../timer.cpp ../matrix.cpp ) +AddTest( + TARGET base_math_bspline + SOURCES math_bspline_test.cpp ../math_bspline.cpp +) diff --git a/source/module_base/test/math_bspline_test.cpp b/source/module_base/test/math_bspline_test.cpp new file mode 100644 index 0000000000..2de14b4282 --- /dev/null +++ b/source/module_base/test/math_bspline_test.cpp @@ -0,0 +1,62 @@ +#include "../math_bspline.h" +#include "gtest/gtest.h" + +/************************************************ + * unit test of class Bspline + ***********************************************/ + +/** + * - Tested Functions: + * - Init + * - norder must be even + * - norder mush be positive + * - Properties + * - \sum_i M_n(u+i) = 1 (i=0,1,2,...n) + * + */ + +class MathBsplineTest : public testing::Test +{ +protected: + ModuleBase::Bspline bp; + int norder; +}; + +TEST_F(MathBsplineTest,Init) +{ + EXPECT_DEATH( + { + norder = 3; // norder must be even + bp.init(norder,0.05,0); + },"" + ); + EXPECT_DEATH( + { + norder = 0; // norder must be positive + bp.init(norder,0.05,0); + },"" + ); +} + +// summation over n is unity +TEST_F(MathBsplineTest,Properties) +{ + int by = 2; + for (norder=2;norder<=20;norder=norder+by) + { + bp.init(norder,1.0,0); + bp.getbspline(0.2); + double sum=0.0; + //std::cout << "\n" << "norder : "<< norder< *b1, complex *b2, complex ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI; ModuleBase::Bspline bsp; bsp.init(norder, 1, 0); - bsp.getbslpine(1.0); + bsp.getbspline(1.0); for(int ix = 0 ; ix < this->nx ; ++ix) { complex fracx=0; @@ -137,4 +137,4 @@ void PW_Basis:: bsplinecoef(complex *b1, complex *b2, complex