.Net ( c# ) 与 Fortran 混合编程实例（二）：杆系结构有限元法——平面桁架解答(3)

第三节  构造有限元法基本方程

3.1  形成未引入边界的总刚度矩阵、总荷载列阵、总边界列阵

新建类，命名为 ClassCalculation，贴入以下代码：

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

using System.Runtime.InteropServices;

namespace Business
{
    public class ClassCalculation
    {
        public Single[,] K = new Single[4, 4];//存放临时单元刚度矩阵
        public Single[,] XY;//杆件（起点、终点）节点的坐标  x、y

        public int rowsBars = ClassBasicInfo.BarsNodes.GetLength(0);//杆件数
        public int rowsNodes = ClassBasicInfo.Coordinate.GetLength(0);//节点数

        //构造函数取得各杆件x1，x2，y1，y2
        public ClassCalculation()
        {
            XY = new Single[rowsBars, 4];

            for (int i = 0; i < rowsBars; i++)
            {
                XY[i, 0] = ClassBasicInfo.Coordinate[ClassBasicInfo.BarsNodes[i, 0] - 1, 0];//起点x
                XY[i, 1] = ClassBasicInfo.Coordinate[ClassBasicInfo.BarsNodes[i, 0] - 1, 1];//起点y
                XY[i, 2] = ClassBasicInfo.Coordinate[ClassBasicInfo.BarsNodes[i, 1] - 1, 0];//终点x
                XY[i, 3] = ClassBasicInfo.Coordinate[ClassBasicInfo.BarsNodes[i, 1] - 1, 1];//终点y
            }
        }

        [DllImport("FORTRANCAL/MatrixCal.dll")]
        public static extern void UnitStiffnessMatrix(ref Single EA, ref Single x1, ref Single x2, ref Single y1, ref Single y2, ref Single K);

        //取得整体坐标下的单元刚度矩阵
        public Single[,] GetK(int i)
        {
            UnitStiffnessMatrix(ref ClassBasicInfo.LinearStiffness[i], ref XY[i, 0], ref XY[i, 2], ref XY[i, 1], ref XY[i, 3], ref K[0,0]);
            return K;
        }

        //取得总刚度矩阵
        public void GetTatalStiffnessMatrix()
        {
            Single[,] K;

            for (int i = 0; i < rowsBars; i++)
            {
                int A = ClassBasicInfo.BarsNodes[i, 0];//起点编号
                int B = ClassBasicInfo.BarsNodes[i, 1];//终点编号

                K = this.GetK(i);//取得单元刚度矩阵
                //分块进行赋值
                //左上块
                ClassBasicInfo.TatalStiffnessMatrix[2 * A - 2, 2 * A - 2] += K[0, 0];
                ClassBasicInfo.TatalStiffnessMatrix[2 * A - 2, 2 * A - 1] += K[0, 1];
                ClassBasicInfo.TatalStiffnessMatrix[2 * A - 1, 2 * A - 2] += K[1, 0];
                ClassBasicInfo.TatalStiffnessMatrix[2 * A - 1, 2 * A - 1] += K[1, 1];
                //右上块
                ClassBasicInfo.TatalStiffnessMatrix[2 * A - 2, 2 * B - 2] += K[0, 2];
                ClassBasicInfo.TatalStiffnessMatrix[2 * A - 2, 2 * B - 1] += K[0, 3];
                ClassBasicInfo.TatalStiffnessMatrix[2 * A - 1, 2 * B - 2] += K[1, 2];
                ClassBasicInfo.TatalStiffnessMatrix[2 * A - 1, 2 * B - 1] += K[1, 3];
                //左下块
                ClassBasicInfo.TatalStiffnessMatrix[2 * B - 2, 2 * A - 2] += K[2, 0];
                ClassBasicInfo.TatalStiffnessMatrix[2 * B - 2, 2 * A - 1] += K[2, 1];
                ClassBasicInfo.TatalStiffnessMatrix[2 * B - 1, 2 * A - 2] += K[3, 0];
                ClassBasicInfo.TatalStiffnessMatrix[2 * B - 1, 2 * A - 1] += K[3, 1];
                //右下块
                ClassBasicInfo.TatalStiffnessMatrix[2 * B - 2, 2 * B - 2] += K[2, 2];
                ClassBasicInfo.TatalStiffnessMatrix[2 * B - 2, 2 * B - 1] += K[2, 3];
                ClassBasicInfo.TatalStiffnessMatrix[2 * B - 1, 2 * B - 2] += K[3, 2];
                ClassBasicInfo.TatalStiffnessMatrix[2 * B - 1, 2 * B - 1] += K[3, 3];
            }
        }

        //取得总荷载列阵（未知支座反力以0代替）
        public void GetTatalLoads()
        {
            for (int i = 0; i < rowsNodes; i++)
            {
                ClassBasicInfo.TatalLoads[2 * i] = ClassBasicInfo.Loads[i, 0];
                ClassBasicInfo.TatalLoads[2 * i + 1] = ClassBasicInfo.Loads[i, 1];
            }
        }

        //取得总边界条件
        public void GetTatalRestraint()
        {
            for (int i = 0; i < rowsNodes; i++)
            {
                ClassBasicInfo.TatalRestraint[2 * i] = ClassBasicInfo.Restraint[i, 0];
                ClassBasicInfo.TatalRestraint[2 * i + 1] = ClassBasicInfo.Restraint[i, 1];
            }
        }

        //接口
        public void GetAll()
        {
            this.GetTatalStiffnessMatrix();
            this.GetTatalLoads();
            this.GetTatalRestraint();
        }
    }
}

取得总体坐标下单元刚度矩阵函数中所用 Fortran 程序为：

subroutine UnitStiffnessMatrix(EA,x1,x2,y1,y2,K)
	implicit none

	!dec$ attributes dllexport::UnitStiffnessMatrix
	!dec$ attributes alias:"UnitStiffnessMatrix"::UnitStiffnessMatrix
	!dec$ attributes reference::EA,x1,x2,y1,y2,K

	
	real(4)::K(4,4)
	real(4)::EA,x1,x2,y1,y2

	real(4)::L
	real(4)::I
	real(4)::cosa,sina

	L = ((x1-x2)**2 + (y1-y2)**2)**0.5
	I = EA/L
	cosa = (x2-x1)/L
	sina = (y2-y1)/L

	call InitK(K,I,cosa,sina)

end subroutine

subroutine InitK(K,I,cos,sin)
	implicit none

	real(4)::K(4,4)
	real(4)::I
	real(4)::cos,sin

	K(1,1) = I * (cos**2)
	K(1,2) = I * (cos*sin)
	K(1,3) = -1.0 * K(1,1)
	K(1,4) = -1.0 * K(1,2)

	K(2,2) = I * (sin**2)
	K(2,3) = K(1,4)
	K(2,4) = -1.0 * K(2,2)

	K(3,3) = K(1,1)
	K(3,4) = K(1,2)

	K(4,4) = K(2,2)

	K(2,1) = K(1,2)
	K(3,1) = K(1,3)
	K(3,2) = K(2,3)
	K(4,1) = K(1,4)
	K(4,2) = K(2,4)
	K(4,3) = K(3,4)

end subroutine

至此，总刚度矩阵、总荷载列阵、总边界列阵已存入静态数组中。

3.2  引入边界条件，求解方程，得到各节点位移、力

新建类，命名为 ClassBoundaryIn，贴入以下代码：

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

using System.Runtime.InteropServices;

namespace Business
{
    public class ClassBoundaryIn
    {
        public Single[,] A;
        public Single[] x;
        public Single[] b;

        int n = ClassBasicInfo.TatalLoads.GetLength(0);

        int num = 0;

        public void Init()
        {
            for (int i = 0; i < n; i++)
            {
                if (ClassBasicInfo.TatalRestraint[i] == false)//自由端
                {
                    num = num + 1;
                }
            }

            A = new Single[num, num];
            x = new Single[num];
            b = new Single[num];
        }

        //整合总刚度矩阵，引入边界条件
        public Single[,] GetA()
        {
            int j = 0;
            for (int i = 0; i < n; i++)
            {
                if (ClassBasicInfo.TatalRestraint[i] == false)//自由端
                {
                    int m = j;
                    for (int k = i; k < n; k++)
                    {
                        if (ClassBasicInfo.TatalRestraint[k] == false)
                        {
                            A[j, m] = ClassBasicInfo.TatalStiffnessMatrix[i, k];
                            A[m, j] = ClassBasicInfo.TatalStiffnessMatrix[k, i];
                            m = m + 1;
                        }
                    }
                    j = j + 1;
                }
            }
            return A;
        }

        //整合总荷载列阵
        public Single[] Getb()
        {
            int j = 0;
            for (int i = 0; i < n; i++)
            {
                if (ClassBasicInfo.TatalRestraint[i] == false)//自由端
                {
                    b[j] = ClassBasicInfo.TatalLoads[i];
                    j = j + 1;
                }
            }
            return b;
        }

        //取得位移列阵
        [DllImport("FORTRANCAL/MatrixCal.dll")]
        public static extern void MatrixNi(ref Single A, ref int n);
        public Single[,] GetA_()
        {
            MatrixNi(ref A[0, 0], ref num);
            return A;
        }

        [DllImport("FORTRANCAL/MatrixCal.dll")]
        public static extern void MatrixJie(ref Single A_, ref int n, ref Single b, ref Single x);
        public Single[] Getx()
        {
            MatrixJie(ref this.GetA_()[0, 0], ref num, ref b[0], ref x[0]);
            return x;
        }

        //形成完整总位移列阵
        public void GetTatalDisplacement()
        {
            int j = 0;
            for (int i = 0; i < n; i++)
            {
                if (ClassBasicInfo.TatalRestraint[i] == false)
                {
                    ClassBasicInfo.TatalDisplacement[i] = x[j];
                    j = j + 1;
                }
                else
                {
                    ClassBasicInfo.TatalDisplacement[i] = 0;
                }
            }
        }

        //形成完整总荷载列阵（包括未知反力）
        [DllImport("FORTRANCAL/MatrixCal.dll")]
        public static extern void MatrixChen(ref Single A, ref int n, ref Single b, ref Single x);
        public void GetTatalLoads()
        {
            MatrixChen(ref ClassBasicInfo.TatalStiffnessMatrix[0, 0], ref n, ref ClassBasicInfo.TatalLoads[0], ref ClassBasicInfo.TatalDisplacement[0]);
        }

        //接口
        public void GetAll()
        {
            this.Init();
            this.GetA();
            this.Getb();
            this.Getx();
            this.GetTatalDisplacement();
            this.GetTatalLoads();
        }
    }
}

三个 Fortran 子例程分别为：

subroutine MatrixNi(a,num)
	implicit none

	!dec$ attributes dllexport::MatrixNi
	!dec$ attributes alias:"MatrixNi"::MatrixNi
	!dec$ attributes reference::a,num

	integer::num
	real(4)::a(num,num)
	integer::i,j,k
	integer::n
	
	n = num
	
	do k=1,N
        a(k,k) = 1.d0/a(k,k)
        do j=1,N
            if(j/=k) a(j,k) = a(k,k)*a(j,k)
        end do
		
		do i=1,N
            do j=1,N
                if(i/=k.and.j/=k) a(j,i) = a(j,i) - a(k,i)*a(j,k)
            end do
            if(i/=k) a(k,i) = -a(k,i)*a(k,k)
        end do
	end do

end subroutine

subroutine MatrixJie(a,num,b,x)
	implicit none

	!dec$ attributes dllexport::MatrixJie
	!dec$ attributes alias:"MatrixJie"::MatrixJie
	!dec$ attributes reference::a,num,b,x

	integer::num
	real(4)::a(num,num)
	real(4)::b(num)
	real(4)::x(num)

	integer::n,i,j

	n = num

	do i = 1,N
		do j = 1,N
			x(i) = x(i) + a(j,i) * b(j)
		end do
	end do
	
end subroutine	 

subroutine MatrixChen(a,num,b,x)
    implicit none

    !dec$ attributes dllexport::MatrixChen
    !dec$ attributes alias:"MatrixChen"::MatrixChen
    !dec$ attributes reference::a,num,b,x

    integer::num
    real(4)::a(num,num)
    real(4)::b(num)
    real(4)::x(num)

    integer::n,i,j

    n = num
    b(:) = 0

    do i = 1,n
        do j = 1,n
            b(i) = b(i) + a(i,j) * x(j)
        end do
    end do

end subroutine

至此，程序已基本完成。