test: add comments and unittest for integral.h

range: source/module_base

Author: pxlxingliang

source/module_base/math_integral.h

@@ -13,7 +13,14 @@ class Integral
     Integral();
     ~Integral();
 
-	// Peize Lin accelerate 2017-10-02
+    /**
+     * @brief simpson integral. Peize Lin accelerate 2017-10-02
+     * 
+     * @param mesh [in] mhe number of grid points (should be odd)
+     * @param func [in] function to be integrated
+     * @param rab [in] a list of interval
+     * @param asum [out] the final integral value
+     */
     static void Simpson_Integral
     (
         const int mesh,
@@ -22,7 +29,14 @@ class Integral
         double &asum
     );
 
-	// Peize Lin accelerate 2017-10-02
+    /**
+     * @brief simpson integral. Peize Lin accelerate 2017-10-02
+     * 
+     * @param mesh [in] mhe number of grid points (should be odd)
+     * @param func [in] function to be integrated
+     * @param dr [in] interval
+     * @param asum [out] the final integral value
+     */
 	static void Simpson_Integral
 	(
 		const int mesh,
@@ -31,7 +45,15 @@ class Integral
 		double &asum
 	);
 
-    // Peize Lin add 2016-02-14
+    /**
+     * @brief simpson integral. Peize Lin add 2016-02-14.
+     *  
+     * 
+     * @param mesh [in] mhe number of grid points 
+     * @param func [in] function to be integrated
+     * @param rab [in] a list of interval
+     * @param asum [out] a list of integral value. asum[i] = integral from 0 to i. The max index of asum is an even (mesh-1 or mesh).
+     */
     static void Simpson_Integral_0toall
     (
         const int mesh,
@@ -40,7 +62,14 @@ class Integral
         double * const asum
     );
 
-    // Peize Lin add 2016-02-14
+    /**
+     * @brief simpson integral. Peize Lin add 2016-02-14.
+     * 
+     * @param mesh [in] mhe number of grid points 
+     * @param func [in] function to be integrated 
+     * @param rab [in] r(i) * dr(i)/di * di or (b-a)/2n 
+     * @param asum [out] a list of integral value. sum[i] = integral from i to mesh-1. 
+     */
     static void Simpson_Integral_alltoinf
     (
         const int mesh,

source/module_base/test/CMakeLists.txt

@@ -18,3 +18,7 @@ AddTest(
   TARGET tool_quit
   SOURCES tool_quit_test.cpp
 )
+AddTest(
+  TARGET base_integral
+  SOURCES math_integral_test.cpp
+)

source/module_base/test/math_integral_test.cpp

@@ -0,0 +1,93 @@
+#include"../math_integral.h"
+#include"gtest/gtest.h"
+
+#include<math.h>
+
+#define PI 3.14159265358979
+#define doublethreshold 1e-12
+
+/************************************************
+*  unit test of class Integral
+***********************************************/
+
+/**
+ * Tested functions:
+ *  - function Simpson_Integral with rab as input
+ *  - function Simpson_Integral with dr as input
+ *  - function Simpson_Integral_0toall
+ *  - function Simpson_Integral_alltoinf
+ * 
+ */
+
+class Simpson_Integral_sinx : public testing::Test
+{
+    /**
+     * test the integral of sinx between [0,PI],
+     * devide to mesh-1 parts
+     */
+    
+    protected:
+
+    double*     func;
+    double*     rab;
+    int         mesh = 10001;
+    double      dr   = PI/(mesh-1);
+    double      asum;
+    double*     asumlist;
+    double      expectvalue = 2.0;
+
+    void SetUp()
+    {
+        func        = new double[mesh];
+        rab         = new double[mesh];
+        asumlist    = new double[mesh];
+
+        for (int i=0;i<mesh;i++)
+        {
+            func[i] = sin(PI*i/(mesh-1));
+            rab[i] = dr;
+        }
+    }
+
+    void TearDown()
+    {
+        delete func;
+        delete rab;
+        delete asumlist;
+    }
+};
+
+TEST_F(Simpson_Integral_sinx,Simpson_Integral_rab)
+{
+    ModuleBase::Integral::Simpson_Integral(mesh,func,rab,asum);
+    EXPECT_NEAR(asum,expectvalue,doublethreshold);
+    EXPECT_DEATH(ModuleBase::Integral::Simpson_Integral(mesh-1,func,rab,asum),"");
+}
+
+TEST_F(Simpson_Integral_sinx,Simpson_Integral_dr)
+{
+    ModuleBase::Integral::Simpson_Integral(mesh,func,dr,asum);
+    EXPECT_NEAR(asum,expectvalue,doublethreshold);
+    EXPECT_DEATH(ModuleBase::Integral::Simpson_Integral(mesh-1,func,dr,asum),"");
+}
+
+TEST_F(Simpson_Integral_sinx,Simpson_Integral_0toall)
+{
+    int halfmesh = round(mesh/2);
+    ModuleBase::Integral::Simpson_Integral_0toall(mesh,func,rab,asumlist);
+    EXPECT_NEAR(asumlist[mesh-1],expectvalue,doublethreshold);
+    EXPECT_NEAR(asumlist[mesh-2],1.9999999506519738901,doublethreshold);
+    EXPECT_NEAR(asumlist[halfmesh],expectvalue/2.0,doublethreshold);
+    EXPECT_NEAR(asumlist[halfmesh-1],0.99968584073722277505,doublethreshold);
+}
+
+
+TEST_F(Simpson_Integral_sinx,Simpson_Integral_alltoinf)
+{
+    int halfmesh = round(mesh/2);
+    ModuleBase::Integral::Simpson_Integral_alltoinf(mesh,func,rab,asumlist);
+    EXPECT_NEAR(asumlist[0],expectvalue,doublethreshold);
+    EXPECT_NEAR(asumlist[1],1.9999999506519754444,doublethreshold);
+    EXPECT_NEAR(asumlist[halfmesh],expectvalue/2.0,doublethreshold);
+    EXPECT_NEAR(asumlist[halfmesh-1],1.0003141592627740053,doublethreshold);
+}