摘译-使用有限元方法求解泊松方程（二）

本内容翻译自《Plasma Simulations by Example》(Lubos Brieda，2019），用于在三维空间内求解泊松方程。本部分基于前一篇文章内容，对有限元代码进行解析。

1 初始化

构造函数定义：

FESolver::FESolver(World &world, int max_it, double tol):
    world{world}, vm{world.vm}, max_solver_it(max_it), tolerance(tol)
{
    /*count number of unknowns*/
    neq = 0;
    /*OPEN nodes are "h" nodes*/
    for (size_t i=0;i<vm.nodes.size();i++)
        if (vm.nodes[i].type==NORMAL ||vm.nodes[i].type==OPEN) 
            neq++;
    cout<<"There are "<<neq<<" unknowns"<<endl;
    /*allocate neq*neq K matrix*/
    K = new double*[neq];
    for (int i=0;i<neq;i++) 
        K[i] = new double[neq];
    cout<<"Allocated "<<neq<<"x"<<neq<<" stiffness matrix"<<endl;

    /*allocate neq*neq J matrix*/
    J = new double*[neq];
    for (int i=0;i<neq;i++) 
        J[i] = new double[neq];
    cout<<"Allocated "<<neq<<"x"<<neq<<" Jacobian matrix"<<endl;

    /*allocate F0 and F1 vectors*/
    F0.reserve(neq);
    F1.reserve(neq);
    cout<<"Allocated two "<<neq<<"x1 force vectors"<<endl;

    n_nodes = vm.nodes.size();
    n_elements = vm.tets.size();

    /*allocate ID vector*/
    ID.reserve(n_nodes);
    cout<<"Allocated "<<n_nodes<<"x1 ID vector"<<endl;

    /*allocate location matrix, n_elements*4 */
    LM = new int*[n_elements];
    for (int e=0;e<n_elements;e++) 
        LM[e] = new int[4];
    cout<<"Allocated "<<n_elements<<"x4 location matrix"<<endl;

    /*allocate NX matrix*/
    NX = new double**[n_elements];
    for (int e=0;e<n_elements;e++) 
    {   
        NX[e] = new double*[4];
        for (int a=0;a<4;a++) 
        NX[e][a] = new double[3];
    }
    cout<<"Allocated "<<n_elements<<"x4x3 NX matrix"<<endl;

    /*solution array*/
    d.reserve(neq);     /*initialized to zero*/

    /*allocate memory for g and uh arrays*/
    g.reserve(n_nodes);
    uh.reserve(n_nodes);
    cout<<"Allocated "<<n_nodes<<"x1 g and uh vector"<<endl;

    detJ.reserve(n_elements);
    cout<<"Allocated "<<n_elements<<"x1 detJ vector"<<endl;

    /*electric field*/
    ef.reserve(n_elements);

    /*set up the ID array note valid values are 0 to neq-1 and -1 indicates "g" node*/
    int P=0;
    for (int n=0;n<n_nodes;n++)
        if (vm.nodes[n].type==NORMAL ||
            vm.nodes[n].type==OPEN) {ID[n]=P;P++;}
        else ID[n]=-1;  /*dirichlet node*/

    /*now set up the LM matrix*/
    for (int e=0;e<n_elements;e++)
        for (int a=0;a<4;a++)   /*tetrahedra*/
            LM[e][a] = ID[vm.tets[e].con[a]];
    cout<<"Built ID and LM matrix"<<endl;

    /*set quadrature points*/
    l[0]=-sqrt(1.0/3.0); l[1]=sqrt(1.0/3.0);
    W[0]=1; W[1]=1;
    n_int = 2;

    /*initialize solver "g" array*/
    for (int n=0;n<n_nodes;n++)
    {
        if (vm.nodes[n].type == INLET) g[n]=0;  /*phi_inlet*/
        else if (vm.nodes[n].type == SPHERE) g[n]=-100; /*phi_sphere*/
        else g[n] = 0;  /*default*/
    }

    /*compute NX matrix*/
    computeNX();

    /*sample assembly code*/
    startAssembly();
    preAssembly();  /*this will form K and F0*/
}