现在的位置: 首页 > 综合 > 正文

第一类边界条件,三角单元剖分,线性插值的位场延拓,LDLT高效求解

2019年03月16日 ⁄ 综合 ⁄ 共 6136字 ⁄ 字号 评论关闭
</pre><pre name="code" class="cpp">/*函数名称:三角单元,线性插值,计算单元系数矩阵ke[3][3]////////////////////////
 *xco[3]:               double,input,store the x coordination,xi,xj,xm
 *yco[3]:               double,input,store the y coordination,yi,yj,ym
 *ke[3][3]:             double,output,store the element coefficient matrix
 */
void uke1 (double xco[], double yco[], double ke[][3])
{   int ce=0;
    double a[3], b[3];
    a[0] = yco[1] - yco[2];
    a[1] = yco[2] - yco[0];
    a[2] = yco[0] - yco[1];
    b[0] = xco[2] - xco[1];
    b[1] = xco[0] - xco[2];
    b[2] = xco[1] - xco[0];
    double deltas = 2.0 * (a[0] * b[1] - a[1] * b[0]);
    for(int i=0; i<3; i++)
    {
        for(int j=0; j<=i; j++)
            ke[i][j]=0;
    }
    for (int i = 0; i < 3; i++)
        for (int j = 0; j <= i; j++)
            ke[i][j] = (a[i] * a[j] + b[i] * b[j]) / deltas;
}
/*函数名称:mbw()///////////////////////////////////////////////////////////////
 *函数目的:计算总体系数矩阵的半带宽,计算每一个单元节点编号两两之差值的最大值,
 *          对得到的所有单元的最大值+1作为半带宽
 *参数说明:ne:                 int型,输入参数,节点总数
 *          index3[3][ne]:      int型,输入参数,单元的节点编号,3*ne
 */
int mbw (int ne, int index3[][16])
{
    int bbw = 0;
    for (int i = 0; i < ne; i++)
    {
        int index1_2 = abs (index3[0][i] - index3[1][i]);
        int index2_3 = abs (index3[1][i] - index3[2][i]);
        int index3_1 = abs (index3[2][i] - index3[0][i]);
        int bbwm = max (max (index1_2, index2_3), index3_1);
        if (bbwm + 1 > bbw)
            bbw = bbwm + 1;
    }
    return bbw;
}

/*函数名称:uk1(int nd,int ne,int bbw,int index3[][ne],double xy[][nd],double sk[nd][bbw])
 *函数功能:三角单元,线性插值时用定带宽存储的方法集成总体矩阵,注意,只存储下三角部分
 *nd:                   int,input,the total number of the node
 *ne:                   int,input,the total number of the element
 *bbw:                  int,input,the bibandwidth of the total matrix
 *index3[3][ne]:        double,input,store the node index of the element,i1,j1,m1,...,ine,jne,mne. 3Xne
 *XY[2][nd]:            double,input,store the x and y coordination,the format is x1,y1,...,xnd,ynd.
 *sk[nd][bbw]:          double,output,store the fixed band coefficient total matrix
 */
void uk1(int nd, int ne, int bbw, int index3[][16], double xy[][15], double sk[][5])
{
    double xco[3], yco[3];      //xco,xcoordination;yco,ycoordination
    double ke[3][3];            // the element stiffness matrix

    for (int l = 0; l < ne; l++)
    {
        for (int j = 0; j < 3; j++)
//对单元的节点作循环,三角形单元有3个节点,取出某个单元的3个节点坐标
        {
            xco[j] = xy[0][index3[j][l]-1];
            yco[j] = xy[1][index3[j][l]-1];
        }
////////////////////////////////////////////////////////////////////////////////
//调用单元矩阵计算函数,产生l单元的单元矩阵,得到单元矩阵ke[3][3]
        uke1 (xco, yco, ke);
///对第l个单元的各个节点循环,把单个单元矩阵叠加到整体矩阵上去。
        for (int j = 0; j < 3; j++)
        {
            int nj = index3[j][l]-1;
            for (int k = 0; k <= j; k++)
            {
                int nk = index3[k][l]-1;
                if (nj >= nk)
                {
                    nk   = nk - nj + bbw-1;
                    sk[nj][nk] = sk[nj][nk] + ke[j][k];
                }
                else
                {
                    nj   = nj - nk + bbw-1;
                    sk[nk][nj] = sk[nk][nj] + ke[j][k];
                    nj   = nj + nk - bbw+1;
                }
            }
        }
    }
}


/*函数名称:ub1(int nd1,int nb1,double u1,int nd,int bbw,double sk[nd][bbw],double u)
 *函数功能: 在定带宽储存的总体系数矩阵和右侧列向量上加上第一类边界条件
 *函数参数:
 *              nd1:            int,input,the node count of the first type boundary condition
 *              nb1[nd1]:       int,input,the node number of the first type boundary condition
 *              u1[nd1]:        double,input,the node field value of the first type boundary condition
 *              nd:             int,input,the total count of the node
 *              bbw:            int,input,bibandwidth
 *              sk[nd][bbw]:    double,input/output,the total coefficient matrix formating the fixed band width, store the first type boundary condition when outputing
 *              u[nd]:          double,output,输出加入第一类边界条件后的(3.3.10)式的右端列向量
 */

void ub1 (int nd1, int nb1[], double u1[], int nd, int bbw, double sk[][5], double u[])
{
    for (int i = 0; i < nd; i++)
        u[i] = 0;
    for (int i = 0; i < nd1; i++)
    {
        int j = nb1[i]-1;
        sk[j][bbw-1] = sk[j][bbw-1] * 1.0E10;
        u[j] = sk[j][bbw-1] * u1[i];
    }
}

/*函数名称:LDLT(double ASubTri[n][bbw],int n,int bbw,double p[n],int sigular)
 *函数功能:对称带型线性方程组的系数矩阵的下三角部分被定带宽储存在矩阵数组A中,利用数组A来解方程组
 *函数参数:
 *a[n][bbw]:            double,input,存放对称带型方程组的系数矩阵的下三角部分
 *n:                    int,input,方程组的阶数
 *bbw:                 int,input,对称带型方程组的半带宽
 *p[n]:                 double,input/output,开始存放方程组的右侧列向量,工作结束时,存放解向量
 *sigular:             int,output,lable symbol,sigular=0,正常结束,sigular=1,系数矩阵奇异
 */
void ldlt (double a[][5], int n, int iw, double p[], int &sigular)
{
    int it;
    sigular = 0;
    for (int i = 0; i < n; i++)
    {
        if (i > iw-1)
        {   it = i - (iw-1);

        }
        else
        {
            it = 0;
        }
        int k = i-1;
        if (i != 0)
        {
            for (int l = it; l <= k; l++)
            {
                int il=(iw-1)-(i-l);
                double b= a[i][il];
                a[i][il]= b / a[l][iw-1];
                p[i]= p[i] - a[i][il] * p[l];
                int mi = l + 1;
                for (int j = mi; j <= i; j++)
                {
                    int ij = j + (iw-1) - i;
                    int jl = l + (iw-1) - j;
                    a[i][ij]=a[i][ij]-a[j][jl]*b;
                }
            }
        }

        if (a[i][iw-1] ==0)
        {
            sigular = 1;
        }
    }
    if (sigular==1)
        return;

    for (int j = 0; j < n; j++)
    {
        if (j > iw-1)
            it = (n-1) - j + iw-1;
        else
            it = n-1;

        int i = (n-1) - j;
        p[i] = p[i] / a[i][iw-1];
        if (j != 0)
        {
            int k = i + 1;
            for (int mj = k; mj <= it; mj++)
            {
                int ij=i-mj+iw-1;
                p[i] = p[i] - p[mj] * a[mj][ij];
            }
        }
    }
    sigular = 0;
    return;
}

int main (int argc, char *argv[])
{
//input parameter
    int nd;             //节点总数
    int ne;             //单元总数
    int nd1;            //第一类边界点的总数
    int index3[3][16];  //每个单元的节点编号,三角形单元,每个单元有3个节点
    int nb1[12];        //第一类边界节点号,共有12个第一类节点

    double xy[2][15];   //每个节点的坐标,共有15个节点,二维空间上每个节点有x和y两个坐标
    double u[15];       //加入第一类边界节点场值后的所有15个节点的场值
    double u1[12];      //第一类边界节点场值
    double sk[15][5];   //存放整体系数矩阵
    double ke[3][3];    //单元系数矩阵

    int sigular;        //整体系数矩阵是否奇异的标志,为0,表示非奇异,为1,表示奇异,方程无解。可以作为return的返回值

    ifstream infile;
    infile.open("input.data",ios_base::in);
    infile>>nd>>ne>>nd1;
    cout<<"nd = "<<nd<<"; ne = "<<ne<<"; nd1 = "<<nd1<<endl;
///////////////////////////////////////////////////////////
    for(int i=0; i<15; i++)
    {
        for(int j=0; j<2; j++)
        {
            infile>>xy[j][i];
        }
    }
///////////////////////////////////////////////////////////
    for(int i=0; i<16; i++)
    {
        for(int j=0; j<3; j++)
        {
            infile>>index3[j][i];
        }
    }
///////////////////////////////////////////////////////////
    for(int i=0; i<12; i++)
    {   infile>>nb1[i]>>u1[i];
    }

    infile.close();

    int bbw=mbw (ne, index3);

    uk1(nd, ne, bbw, index3, xy, sk);
///////////////////////////////////////////////////////////
    for(int i=0; i<nd; i++)
    {
        for(int j=0; j<bbw; j++)
        {
            cout<<"sk["<<setw(2)<<i<<"]["<<setw(2)<<j<<"]="<<setw(8)<<sk[i][j]<<"\t";
        }
        cout<<endl;
    }
///////////////////////////////////////////////////////////

//功能:在定带宽存储的总体系数矩阵和右侧列向量上加上第一类边界条件。
    ub1 (nd1, nb1, u1, nd, bbw, sk, u);
    cout<<endl<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
    for(int i=0; i<nd; i++)
    {
        for(int j=0; j<bbw; j++)
        {
            cout<<"sk["<<setw(2)<<i<<"]["<<setw(2)<<j<<"]="<<setw(8)<<sk[i][j]<<"\t";
        }
        cout<<endl;
    }

    cout<<endl<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
    for (int i=0; i<nd; i++)
    {
        cout<<"u["<<i<<"]="<<setw(5)<<u[i]<<endl;
    }

// 功能:对称带型线性方程组的系数矩阵的下三角部分被定带宽地储存杂矩阵数组SK中,利用SK数组来解方程组。
    ldlt (sk,nd,bbw,u,sigular);
    cout<<endl<<"*************************************************"<<endl;
    for (int i=0; i<nd; i++)
    {
        cout<<"u["<<setw(2)<<i<<"]= "<<setw(12)<<u[i];
        if(i==5||i==8||i==11)cout<<"\t***";
        cout<<endl;
    }
    cout<<endl<<"sigular = "<<sigular<<endl;
}

以下是input.data文件

  15 16 12
  0.0  0.0
  0.0  1.0
  0.0  2.0
  1.0  0.5
  1.0  1.2
  1.0  2.0
  2.0  1.0
  2.0  1.5
  2.0  2.0
  3.0  0.7
  3.0  1.3
  3.0  2.0
  4.0  0.5
  4.0  1.2
  4.0  2.0
  1  4  2
  5  2  4
  5  3  2
  5  6  3
  5  4  7
  5  7  8
  5  8  6
  8  9  6
 11  7 10
  8  7 11
  8 11 12
  8 12  9
 10 13 14
 11 10 14
 11 14 15
 11 15 12
  1 1.00
  2 0.50
  3 0.33
  4 0.46
  6 0.30
  7 0.25
  9 0.23
 10 0.14
 12 0.17
 13 0.08
 14 0.11
 15 0.12

参考文献:

1、Fem in Geophysics,p94-108,科学出版社

2、FORTRAN 算法汇编 第一分册 刘德贵 费景高 p204,

3、电子计算机常用算法,中国科学院沈阳计算技术研究所等编科学出版社,p132~136,对深入理解ldlt函数有特别帮助

抱歉!评论已关闭.