.Net ( c# ) 与 Fortran 混合编程实例（二）：杆系结构有限元法—

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

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、ypublic int rowsBars = ClassBasicInfo.BarsNodes.GetLength(0);//杆件数public int rowsNodes = ClassBasicInfo.Coordinate.GetLength(0);//节点数//构造函数取得各杆件x1，x2，y1，y2public ClassCalculation(){XY = new Single[rowsBars, 4];for (int i = 0; i < rowsBars; i++){XY[i, 0] = ClassBasicInfo.Coordinate[ClassBasicInfo.BarsNodes[i, 0] - 1, 0];//起点xXY[i, 1] = ClassBasicInfo.Coordinate[ClassBasicInfo.BarsNodes[i, 0] - 1, 1];//起点yXY[i, 2] = ClassBasicInfo.Coordinate[ClassBasicInfo.BarsNodes[i, 1] - 1, 0];//终点xXY[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,Kreal(4)::K(4,4)real(4)::EA,x1,x2,y1,y2real(4)::Lreal(4)::Ireal(4)::cosa,sinaL = ((x1-x2)**2 + (y1-y2)**2)**0.5I = EA/Lcosa = (x2-x1)/Lsina = (y2-y1)/Lcall InitK(K,I,cosa,sina)end subroutinesubroutine InitK(K,I,cos,sin)implicit nonereal(4)::K(4,4)real(4)::Ireal(4)::cos,sinK(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,numinteger::numreal(4)::a(num,num)integer::i,j,kinteger::nn = numdo k=1,Na(k,k) = 1.d0/a(k,k)do j=1,Nif(j/=k) a(j,k) = a(k,k)*a(j,k)end dodo i=1,Ndo j=1,Nif(i/=k.and.j/=k) a(j,i) = a(j,i) - a(k,i)*a(j,k)end doif(i/=k) a(k,i) = -a(k,i)*a(k,k)end doend doend subroutine

subroutine MatrixJie(a,num,b,x)implicit none!dec$ attributes dllexport::MatrixJie!dec$ attributes alias:"MatrixJie"::MatrixJie!dec$ attributes reference::a,num,b,xinteger::numreal(4)::a(num,num)real(4)::b(num)real(4)::x(num)integer::n,i,jn = numdo i = 1,Ndo j = 1,Nx(i) = x(i) + a(j,i) * b(j)end doend doend subroutine

subroutine MatrixChen(a,num,b,x)implicit none!dec$ attributes dllexport::MatrixChen!dec$ attributes alias:"MatrixChen"::MatrixChen!dec$ attributes reference::a,num,b,xinteger::numreal(4)::a(num,num)real(4)::b(num)real(4)::x(num)integer::n,i,jn = numb(:) = 0do i = 1,ndo j = 1,nb(i) = b(i) + a(i,j) * x(j)end doend doend subroutine

至此，程序已基本完成。

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

.Net ( c# ) 与 Fortran 混合编程实例（二）：杆系结构有限元法——平面桁架解答(2)...
第二节桁架结构基本信息录入与保存此节之后会结合一个算例.结构如图,EA = 300000 kN 杆件编号.节点编号.坐标系: 2.1 文本信息 bars.txt 1,1,3,300000 2, ...
QML和C++混合编程--（二）
文章目录 QML和C++混合编程--(二)之QML连接C++方法和信号简介: 一.修改app.qml 文件二.修改c++中的类 QML和C++混合编程–(二)之QML连接C++方法和信号简介: ...
MFC编程实例二：进度条的使用
MFC编程实例二:进度条的使用 (以下蓝色字体为手动添加的代码) (1) 创建一个基于对话框的应用程序. (2) 在对话框界面添加一个文本框.一个进度条控件(Progress控件)和一个按钮. ...
python和C语言混合编程实例
转载自:http://www.jb51.net/article/50633.htm 这篇文章主要介绍了python和C语言混合编程实例,文中开发了一个tcp端口ping程序来介绍混合编程,需要的朋友可 ...
平面/空间杆系结构有限元编程计算（MATLAB）
GitHub链接:MATLAB源码博主QQ:915339719,有问题可以随时交流程序简介本程序可以对绝大多数的平面/空间杆系结构进行静力学计算.自振频率计算,以及平面杆系结构的稳定性计算,具有 ...
python fortran混合编程_python fortran c 混合编程
背景因为编程需要使用python和fortran(或者c,c++)混合编程,所以进行了一些资料查询和实验.首先尝试了f2py,但是这个项目已经很久没有更新,为了以后兼容性考虑(后期可能会使用到c++ ...
ESP8266-Arduino编程实例-二维码（QR Code）生成与TFT LCD显示
二维码(QR Code)生成与TFT LCD显示 QR 码(快速响应码的缩写)是一种矩阵条码(或二维条码,由日本汽车公司 Denso Wave 于 1994 年发明.条形码是一种机器可读的光学标签,可 ...
MPI和OpenMP混合编程实例
MPI和OpenMP混合编程计算pi值大体思路就是先将计算区域划分成100000 份,然后将计算分成进程数nproc(nproc)份,每一份有nbin个计算. 步长为step = 1.0 / NB ...
windows机制下的游戏编程实例二
通过调用DirectX SDK库实现游戏空间的三维化,也就是我们通常玩的网游,鼠标拖着主角到处跑去完成任务(个人觉得这种游戏超无聊).下面先给出效果图: 还是先给出其文件目录吧,从总体把握实例,也是很 ...

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

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

最新文章

热门文章