增加了稀疏矩阵的迭代求解功能，并对稀疏矩阵和稠密矩阵进行了封装

Author: yinweijie

MatrixCoeff.cpp

@@ -1,54 +1,67 @@
 #include "MatrixCoeff.h"
+#include "MatrixWrapper.h"
 
-// a_{p} T_{p}=a_{L} T_{L}+a_{R} T_{R}+S_{u}
-void MatrixCoeff::initMatrix(const Mesh* mesh, const Boundary& boundary, const Source& source)
+// a_p T_p = a_L T_L + a_R T_R + ... + S_{u}
+void MatrixCoeff::initDenseMatrix()
 {
-    int N = mesh->get_N();
+    init(&dense_matrix_wrapper);
+
+    b_m = Su;
+}
+
+void MatrixCoeff::initSparseMatrix()
+{
+    SparseMatrix<double> A_m = sparse_matrix_wrapper.getMatrix();
+
+    // 每列预留5个元素的空间，用于插入
+    // ref. http://eigen.tuxfamily.org/dox/group__TutorialSparse.html - Filling a sparse matrix
+    A_m.reserve(VectorXi::Constant(N, 5));
+
+    init(&sparse_matrix_wrapper);
+
+    A_m.makeCompressed(); // optional
+
+    b_m = Su;
+}
 
+void MatrixCoeff::init(MatrixInterface* matrix)
+{
     for(int i = 0; i < N; i++)
     {
         int i_l = mesh->left_of(i);
         int i_r = mesh->right_of(i);
         int i_t = mesh->top_of(i);
         int i_b = mesh->bottom_of(i);
 
-        A_m(i, i) = aP[i];
+        matrix->setNum(i, i, aP[i]);
 
         if(!mesh->is_at_left_boundary(i))
         {
-            A_m(i, i_l) = -aL[i];
+            matrix->setNum(i, i_l, -aL[i]);
         }
         if(!mesh->is_at_right_boundary(i))
         {
-            A_m(i, i_r) = -aR[i];
+            matrix->setNum(i, i_r, -aR[i]);
         }
         if(!mesh->is_at_bottom_boundary(i))
         {
-            A_m(i, i_b) = -aB[i];
+            matrix->setNum(i, i_b, -aB[i]);
         }
         if(!mesh->is_at_top_boundary(i))
         {
-            A_m(i, i_t) = -aT[i];
+            matrix->setNum(i, i_t, -aT[i]);
         }
     }
-
-    b_m = Su;
-
-    cout << "A_m: " << endl << A_m << endl;
-    cout << endl;
-
-    cout << "b_m: " << endl << b_m << endl;
-    cout << endl;
 }
 
-MatrixCoeff& MatrixCoeff::addConvectionTerm(const Mesh* mesh, const Boundary& boundary, const Source& source)
+MatrixCoeff& MatrixCoeff::addConvectionTerm()
 {
-    double T_l = boundary.T_left;
-    double T_r = boundary.T_right;
-    double T_b = boundary.T_bottom;
-    double T_t = boundary.T_top;
+    double T_l = boundary->T_left;
+    double T_r = boundary->T_right;
+    double T_b = boundary->T_bottom;
+    double T_t = boundary->T_top;
 
-    double q_w = boundary.q_w;
+    double q_w = boundary->q_w;
 
     double Ax = mesh->get_Ax();
     double Ay = mesh->get_Ay();
@@ -59,8 +72,6 @@ MatrixCoeff& MatrixCoeff::addConvectionTerm(const Mesh* mesh, const Boundary& bo
     const VectorXd& F_b = mesh->get_F_b();
     const VectorXd& F_t = mesh->get_F_t();
 
-    int N = mesh->get_N();
-
     aL += F_l.cwiseMax(VectorXd::Zero(N));
     aR += (-F_r).cwiseMax(VectorXd::Zero(N));
     aB += F_b.cwiseMax(VectorXd::Zero(N));
@@ -99,13 +110,13 @@ MatrixCoeff& MatrixCoeff::addConvectionTerm(const Mesh* mesh, const Boundary& bo
     return *this;
 }
 
-MatrixCoeff& MatrixCoeff::addDiffusionTerm(const Mesh* mesh, const Boundary& boundary, const Source& source)
+MatrixCoeff& MatrixCoeff::addDiffusionTerm()
 {
-    double T_l = boundary.T_left;
-    double T_r = boundary.T_right;
-    double T_b = boundary.T_bottom;
-    double T_t = boundary.T_top;
-    double q_w = boundary.q_w;
+    double T_l = boundary->T_left;
+    double T_r = boundary->T_right;
+    double T_b = boundary->T_bottom;
+    double T_t = boundary->T_top;
+    double q_w = boundary->q_w;
 
     double Ax = mesh->get_Ax();
     double Ay = mesh->get_Ay();
@@ -116,8 +127,6 @@ MatrixCoeff& MatrixCoeff::addDiffusionTerm(const Mesh* mesh, const Boundary& bou
     const VectorXd& DA_B = mesh->get_DA_B();
     const VectorXd& DA_T = mesh->get_DA_T();
 
-    int N = mesh->get_N();
-
     aL += DA_L;
     aR += DA_R;
     aB += DA_B;
@@ -157,9 +166,9 @@ MatrixCoeff& MatrixCoeff::addDiffusionTerm(const Mesh* mesh, const Boundary& bou
     return *this;
 }
 
-MatrixCoeff& MatrixCoeff::addSourceTerm(const Mesh* mesh, const Boundary& boundary, const Source& source)
+MatrixCoeff& MatrixCoeff::addSourceTerm()
 {
-    double S_bar = source.S_bar;
+    double S_bar = source->S_bar;
 
     const VectorXd& V = mesh->get_V();
 
@@ -168,4 +177,44 @@ MatrixCoeff& MatrixCoeff::addSourceTerm(const Mesh* mesh, const Boundary& bounda
     Su += S_bar * V;
 
     return *this;
+}
+
+void MatrixCoeff::DebugSolve()
+{
+    initDenseMatrix();
+
+    MatrixXd& A_m = dense_matrix_wrapper.getMatrix();
+
+    // 求解矩阵
+    x = A_m.fullPivLu().solve(b_m);
+
+    // 输出结果
+    cout << "A_m: " << endl << A_m << endl;
+    cout << endl;
+
+    cout << "b_m: " << endl << b_m << endl;
+    cout << endl;
+
+    cout << "Solution: " << endl << x << endl;
+}
+
+void MatrixCoeff::Solve()
+{
+    initSparseMatrix();
+
+    SparseMatrix<double>& A_m = sparse_matrix_wrapper.getMatrix();
+
+    // ref. http://eigen.tuxfamily.org/dox/classEigen_1_1BiCGSTAB.html
+    // ref. http://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html#TutorialSparseSolverConcept
+    // ref. http://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html
+
+    // BiCGSTAB<SparseMatrix<double>, Eigen::IdentityPreconditioner> solver;
+    BiCGSTAB<SparseMatrix<double>, Eigen::IncompleteLUT<double>> solver;
+    solver.compute(A_m);
+    x = solver.solve(b_m);
+
+    std::cout << "#iterations:     " << solver.iterations() << std::endl;
+    std::cout << "estimated error: " << solver.error()      << std::endl;
+
+    std::cout << "x: " << std::endl << x << std::endl;
 }

MatrixCoeff.h

@@ -2,26 +2,43 @@
 #define __MATRIX_COEFF__
 
 #include <iostream>
+#include <Eigen/Sparse>
+#include <Eigen/IterativeLinearSolvers>
 #include "Inputs.h"
 #include "Mesh.h"
+#include "MatrixWrapper.h"
 
+using Eigen::SparseMatrix;
+using Eigen::BiCGSTAB;
 using namespace std;
 
 // N: 网格数
 class MatrixCoeff
 {
 private:
+    const Mesh* mesh;
+    const Boundary* boundary;
+    const Source* source;
+    // 网格数量
+    int N;
+
     VectorXd aL, aR, aB, aT, aP, Sp, Su;
 
-    MatrixXd A_m;
     VectorXd b_m;
     VectorXd x;
 
+    /**
+     * MatrixXd和SparseMatrix的封装类，封装后的Wrapper类共同继承接口MatrixInterface，
+     * 并实现接口函数MatrixInterface::setNum，这个函数用来往矩阵中插入元素
+     */
+    DenseMatrixWrapper dense_matrix_wrapper;
+    SparseMatrixWrapper sparse_matrix_wrapper;
+
 public:
-    MatrixCoeff(const Mesh* mesh, const Boundary& boundary, const Source& source)
+    MatrixCoeff(const Mesh* mesh, const Boundary* boundary, const Source* source) 
+        : mesh(mesh), boundary(boundary), source(source), N(mesh->get_N()),
+          dense_matrix_wrapper(N, N), sparse_matrix_wrapper(N, N)
     {
-        int N = mesh->get_N();
-
         aL = VectorXd::Zero(N);
         aR = VectorXd::Zero(N);
         aB = VectorXd::Zero(N);
@@ -31,21 +48,30 @@ class MatrixCoeff
         Sp = VectorXd::Zero(N);
         Su = VectorXd::Zero(N);
 
-        A_m = MatrixXd::Zero(N, N);
         b_m = VectorXd::Zero(N);
         x = VectorXd::Zero(N);
     }
 
-    MatrixXd& get_A_m() { return A_m; }
     VectorXd& get_b_m() { return b_m; }
     VectorXd& get_x() { return x; }
 
-    MatrixCoeff& addConvectionTerm(const Mesh* mesh, const Boundary& boundary, const Source& source);
-    MatrixCoeff& addDiffusionTerm(const Mesh* mesh, const Boundary& boundary, const Source& source);
-    MatrixCoeff& addSourceTerm(const Mesh* mesh, const Boundary& boundary, const Source& source);
-    
-    // Matrix系数初始化
-    void initMatrix(const Mesh* mesh, const Boundary& boundary, const Source& source);
+    // MatrixXd系数初始化
+    void initDenseMatrix();
+
+    // SparseMatrix系数初始化
+    void initSparseMatrix();
+
+    void init(MatrixInterface* matrix);
+
+    MatrixCoeff& addConvectionTerm();
+    MatrixCoeff& addDiffusionTerm();
+    MatrixCoeff& addSourceTerm();
+
+    // 系数矩阵存在非稀疏矩阵MatrixXd中，方便打印输出
+    void DebugSolve();
+
+    // 系数矩阵存在稀疏矩阵SparseMatrix中，用于迭代求解
+    void Solve();
 };
 
 #endif

MatrixInterface.cpp

@@ -0,0 +1,7 @@
+#include "MatrixInterface.h"
+
+// 纯虚析构函数必须保留一个类外实现，因为子类在析构的时候会调用父类的析构函数，如果不定义会编译报错
+MatrixInterface::~MatrixInterface()
+{
+
+}

MatrixInterface.h

@@ -0,0 +1,11 @@
+#ifndef __MATRIX_INTERFACE__
+#define __MATRIX_INTERFACE__
+
+class MatrixInterface
+{
+public:
+    virtual void setNum(int row, int col, double value) = 0;
+    virtual ~MatrixInterface() = 0;
+};
+
+#endif

MatrixWrapper.h

@@ -0,0 +1,48 @@
+#ifndef __MATRIX_WRAPPER__
+#define __MATRIX_WRAPPER__
+
+#include <iostream>
+#include <Eigen/Sparse>
+#include "MatrixInterface.h"
+
+using Eigen::MatrixXd;
+using Eigen::SparseMatrix;
+
+class DenseMatrixWrapper : public MatrixInterface
+{
+private:
+    MatrixXd m_dense_matrix;
+public:
+    DenseMatrixWrapper(int row, int col) : m_dense_matrix(row, col) { }
+
+    void setNum(int row, int col, double value)
+    {
+        m_dense_matrix(row, col) = value;
+    }
+
+    MatrixXd& getMatrix()
+    {
+        return m_dense_matrix;
+    }
+};
+
+class SparseMatrixWrapper : public MatrixInterface
+{
+private:
+    SparseMatrix<double> m_sparse_matrix;
+public:
+    SparseMatrixWrapper(int row, int col) : m_sparse_matrix(row, col) { }
+
+    void setNum(int row, int col, double value)
+    {
+        // ref. http://eigen.tuxfamily.org/dox/group__TutorialSparse.html - Filling a sparse matrix
+        m_sparse_matrix.insert(row, col) = value; // alternative: A.coeffRef(i,j) += v_ij;
+    }
+
+    SparseMatrix<double>& getMatrix()
+    {
+        return m_sparse_matrix;
+    }
+};
+
+#endif