有限元编程的c++实现算例
有限元作业
有限元作业1.输出1-100之间能被3整除的数,并求这些数的累计之和。
解答：在workplace中输入如下代码：READ*,NN=100SUM=0.0DO I=1,NIF (MOD(I,3)==0) THENSUM=SUM+IWRITE(6,*)IELSESUM=SUM+0END IFEND DOWRITE(6,*) SUMEND按Ctrl+F5键，输入任意一个N，因为N为给定的100，无论输入任意数据都行，在界面中输入数字1，得出结果为：2.某大学按学分收费,如果不超过12学分,应交纳学费4000元,如果超过12学分,则每超过1个学分,加收500元学费,输出学分和应缴纳的学费.解答：在workplace中输入如下代码：READ*,NSUM=4000DO I=1,NIF(I>=12)THENSUM=4000+(I-12)*500ENDIFWRITE(6,*)I,SUMENDDOEND按Ctrl+F5键，输入任意一个N ，在界面中输入数字20，得出结果为：3.用迭代法求解方程：038223=++-x x x 解答：在workplace 中输入如下代码：READ *,X0 N=0 X0=3 DO I=1,NF1=X*X*X-2*X*X+8*X+3 F2=3*X*X-4*X+8 X=X0-F1/F2 N=N+1IF (ABS (X-X0)>=0.00001)THEN X=X ENDIF ENDDO WRITE (6,*)X,F1 END按Ctrl+F5键，输入任意一个N ，在界面中输入数字0.02，得出结果为：4. 编写任意大小的两个矩阵相乘的程序，并算一具体实例5. 哥德巴赫提出，一个不小于6的偶数必定能表示为两个素数之和，请将6-100之间的全部偶数表示为两个素数之和6. 已知2)sinh(xx e e x --=利用级数展开式求出x=1~5时的sinh(x)7.已知矩阵，求出该矩阵对角线上元素的和⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=359162173895079684835296153********2A 8. 利用Fortran 语言编译一程序打印乘法口诀表 9. 设计一个程序计算：S=9!+ 7!+ 5!+3!+1! 10. 编译一程序求逆矩阵并用一实例验证 11. 不限方法求积分⎰10dxe x12. 已知杆的E I =6，画连续梁的内力图kN.m313.画刚架的内力图14.已知各杆的弹性模量和横截面面积均相同，为E=29.5´104N/mm2, A=100mm2,求各杆的内力。
有限元方法编程
有限元方法编程简介有限元方法（Finite Element Method, FEM）是一种数值分析方法，用于求解连续介质力学问题。
它将复杂的物理问题离散化为有限数量的简化子问题，并通过构建逼近函数来近似求解。
有限元方法广泛应用于结构力学、流体力学、电磁场等领域。
在本文中，我们将探讨有限元方法的编程实现。
我们将介绍有限元模型的建立、离散化、边界条件的处理以及求解过程等关键步骤。
有限元模型的建立在使用有限元方法求解物理问题之前，首先需要建立一个合适的有限元模型。
这包括选择合适的几何形状和材料属性，并对其进行离散化。
几何形状几何形状是指待求解区域的外部形状，可以是简单几何形状如矩形、圆形，也可以是复杂几何形状如曲线、曲面。
在编程中，我们可以使用各种数学函数或CAD软件来表示几何形状，并将其转换为计算机可识别的格式。
材料属性材料属性是指待求解区域内物质的力学性质，如弹性模量、泊松比等。
在编程中，我们需要将这些属性赋予模型，并将其储存在合适的数据结构中。
离散化离散化是将连续的物理问题转化为离散的子问题。
在有限元方法中，我们通常使用三角形或四边形网格来离散化几何形状。
这些网格被称为有限元网格，每个单元代表一个子问题。
有限元模型的编程实现有限元模型的编程实现主要包括以下几个关键步骤：网格生成、单元属性赋值、边界条件处理以及求解过程。
网格生成在编程中，我们可以使用各种算法来生成有限元网格。
其中最常用的算法是Delaunay三角剖分算法和四边形剖分算法。
这些算法可以根据几何形状和所需精度生成合适的网格，并将其储存在数据结构中。
单元属性赋值每个有限元单元都具有一组属性，包括几何信息和材料信息。
在编程中，我们需要将这些属性赋予每个单元，并将其储存在合适的数据结构中。
这些数据结构通常是数组或矩阵。
边界条件处理边界条件是指在求解过程中所施加的约束条件。
它们可以是位移边界条件、力边界条件或温度边界条件等。
在编程中，我们需要将这些边界条件应用于有限元模型，并将其储存在合适的数据结构中。
有限元编程算例(fortran)
有限元编程算例(Fortran)本程序通过Fortran语言编写，程序在Intel Parallel Studio XE 2013 with VS2013中成功运行，程序为《计算力学》（龙述尧等编）一书中的源程序，仅作研究学习使用，省去了敲写的麻烦。
源程序为：!Page149COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200) OPEN(5,FILE='DATAIN')!OPEN(6,FILE='DATAOUT',STATUS='NEW')CALL DATAIF(IND.EQ.0)GOTO 10EO=EO/(1.0-UN*UN)UN=UN/(1.0-UN)10 CALL TOTSTICALL LOADCALL SUPPORCALL SOLVEQCALL STRESSPAUSE!STOPENDSUBROUTINE DATACOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)READ(5,*)NJ,NE,NZ,NDD,NPJ,INDNJ2=NJ*2NPJ1=NPJ+1READ(5,*)EO,UN,GAMA,TEREAD(5,*)((JM(I,J),J=1,3),I=1,NE)READ(5,*)((CJZ(I,J),J = 1,2),I=1,NJ)!Page150READ(5,*)(NZC(I),I=1,NZ)READ(5,*)((PJ(I,J),J=1,2),I=1,NPJ1)WRITE(6,10)(I,(CJZ(I,J),J=1,2),I=1,NJ)10 FORMA T(4X,2HNO,6X,1HX,6X,1HY/(I6,2X,F7.2,F7.2))RETURNENDSUBROUTINE ELEST(MEO,IASK)COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200) IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)CM=CJZ(JE,1)-CJZ(IE,1)BM=CJZ(IE,2)-CJZ(JE,2)CJ=CJZ(IE,1)-CJZ(ME,1)BJ=CJZ(ME,2)-CJZ(IE,2)AE=(BJ*CM-BM*CJ)/2.0IF(IASK.LE.1) GOTO 50DO 10 I=1,3DO 10 J=1,6B(I,J)=0.010 CONTINUEB(1,1)=-BJ-BMB(1,3)=BJB(1,5)=BMB(2,2)=-CJ-CMB(2,4)=CJB(2,6)=CMB(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)!Page151B(3,6)=B(1,5)DO 20 I=1,3DO 20 J=1,6B(I,J)=B(I,J)/(2.0*AE) 20 CONTINUED(1,1)=EO/(1.0-UN*UN)D(1,2)=EO*UN/(1.0-UN*UN) D(2,1)=D(1,2)D(2,2)=D(1,1)D(1,3)=0.0D(2,3)=0.0D(3,1)=0.0D(3,2)=0.0D(3,3)=EO/(2.0*(1.0+UN))DO 30 I=1,3DO 30 J=1,6S(I,J)=0.0DO 30 K=1,3S(I,J)=S(I,J)+D(I,K)*B(K,J)30 CONTINUEIF(IASK.LE.2) GOTO 50DO 40 I=1,6DO 40 J=1,6EKE(I,J)=0.0DO 40 K=1,3!**********************************Exchange B And S***********************************************EKE(I,J)=EKE(I,J)+B(K,I)*S(K,J)*AE*TE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE TOTSTICOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)!Page152DO 20 I=1,NJ2DO 20 J=1,NDDTKZ(I,J)=0.020 CONTINUE!*************Not Understanded*****************************DO 30 MEO=1,NECALL ELEST(MEO,3)DO 30 I=1,3DO 30 II=1,2LH=2*(I-1)+IILDH=2*(JM(MEO,I)-1)+IIDO 30 J=1,3DO 30 JJ=1,2L=2*(J-1)+JJLZ=2*(JM(MEO,J)-1)+JJLD=LZ-LDH+1IF(LD.LE.0) GOTO 30TKZ(LDH,LD)=TKZ(LDH,LD)+EKE(LH,L)30 CONTINUERETURNENDSUBROUTINE LOADCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 10 I=1,NJ2P(I)=0.010 CONTINUEIF(NPJ.EQ.0) GOTO 30DO 20 I=1,NPJI1=I+1J=IFIX(PJ(I1,2))P(J)=PJ(I1,1)20 CONTINUE30 IF(GAMA.LE.0.0) GOTO 50!Page153DO 40 MEO=1,NECALL ELEST(MEO,1)PE=-GAMA*AE*TE/3.0IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)P(2*IE)=P(2*IE)+PEP(2*JE)=P(2*JE)+PEP(2*ME)=P(2*ME)+PE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE SUPPORCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 60 I=1,NZMZ=NZC(I)TKZ(MZ,1)=1.0DO 10 J=2,NDDTKZ(MZ,J)=0.010 CONTINUEIF(MZ-NDD)20,20,3020 JO=MZGOTO 4030 JO=NDD40 DO 50 J = 2,JOJ1=MZ-JTKZ(J1+1,J)=0.050 CONTINUEP(MZ)=0.060 CONTINUERETURNEND!Page154SUBROUTINE SOLVEQCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)NJ1=NJ2-1DO 50 K=1,NJ1IF(NJ2-K-NDD+1)10,10,2010 IM=NJ2GOTO 3020 IM=K+NDD-130 K1=K+1DO 50 I=K1,IML=I-K+1C=TKZ(K,L)/TKZ(K,1)LD1=NDD-L+1DO 40 J=1,LD1M=J+I-KTKZ(I,J)=TKZ(I,J)-C*TKZ(K,M)40 CONTINUEP(I)=P(I)-C*P(K)50 CONTINUEP(NJ2)=P(NJ2)/TKZ(NJ2,1)DO 100 I1 = 1,NJ1I=NJ2-I1!************************************************************************下面一行可能出错IF(NDD-NJ2+I-1)60,60,7060 JO=NDDGOTO 8070 JO=NJ2-I+180 DO 90 J=2,JOLH=J+I-1P(I)=P(I)-TKZ(I,J)*P(LH)90 CONTINUEP(I)=P(I)/TKZ(I,1)100 CONTINUE!Page155WRITE(6,110)(I,P(2*I-1),P(2*I),I=1,NJ)!************************************************************************************ 110 FORMA T(2X,3HJD=,3X,2HU=,12X,2HV=/(I4,3X,F16.7,3X,F16.7))RETURNENDSUBROUTINE STRESSCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200) DIMENSION WY(6),YL(3)DO 60 MEO=1,NECALL ELEST(MEO,2)DO 10 I=1,3DO 10 J=1,2LH=2*(I-1)+JLDH=2*(JM(MEO,I)-1)+JWY(LH)=P(LDH)10 CONTINUEDO 20 I=1,3YL(I)=0.0DO 20 J=1,6YL(I)=YL(I)+S(I,J)*WY(J)20 CONTINUESIGX=YL(1)SIGY=YL(2)TOXY=YL(3)PYL=(SIGX+SIGY)/2.0SIG=(SIGX-SIGY)**2/4.0+TOXY*TOXYRYL=SQRT(SIG)SIG1=PYL+RYLSIG2=PYL-RYLIF(SIGY.EQ.SIG2) GOTO 30CETA1=TOXY/(SIGY-SIG2)CETA=90.0-57.29578*ATAN(CETA1)GOTO 40!Page15630 CETA=0.040 WRITE(6,50)MEO,SIGX,SIGY,TOXY,SIG1,SIG2,CETA50FORMA T(4X,2HE=,I3/2X,3HSX=,F11.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HS1=,F11.3,3X,3HS2=,F11. 3,3X,4HCET=,F11.3)!50FORMA T(4X,2HE=,I3/2X,3HSX=,Fll.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HSl=,Fll.3,3X,3HS2=,F11.3,3 X,4HCET=,F11.3)60 CONTINUERETURNEND输入文件为datain28,36,9,10,4,01,0.17,0,11,5,22,5,62,6,33,6,73,7,44,7,85,9,66,9,106,10,77,10,117,11,88,11,129,13,1010,13,1410,14,11 11,14,15 11,15,12 12,15,16 13,17,14 14,17,18 14,18,15 15,18,19 15,19,16 16,19,20 17,21,18 18,21,22 18,22,19 19,22,23 19,23,20 20,23,24 21,25,22 22,25,26 22,26,23 23,26,27 23,27,24 24,27,28 0,61,62,63,60,51,52,53,50,41,42,43,40,31,32,33,30,21,22,23,20,11,12,13,10,01,02,03,07,15,23,31,39,47,49,50,550,0-5E4,2-10E4,4-10E4,6-5E4,8输出结果为：DATAOUTNO X Y1 0.00 6.002 1.00 6.003 2.00 6.004 3.00 6.005 0.00 5.006 1.00 5.007 2.00 5.008 3.00 5.009 0.00 4.0010 1.00 4.0011 2.00 4.0012 3.00 4.0013 0.00 3.0014 1.00 3.0015 2.00 3.0016 3.00 3.0017 0.00 2.0018 1.00 2.0019 2.00 2.0020 3.00 2.0021 0.00 1.0022 1.00 1.0023 2.00 1.0024 3.00 1.0025 0.00 0.0026 1.00 0.0027 2.00 0.0028 3.00 0.00JD= U= V=1 -29766.873 -1173917.7502 -14003.185 -1174018.8753 -3753.270 -1179518.1254 0.000 -1181719.7505 -26382.471 -1072681.5006 -10746.993 -1073615.0007 -2064.593 -1082360.7508 0.000 -1085873.2509 -13536.995 -964010.12510 3372.794 -970055.12511 7268.415 -989269.12512 0.000 -998401.81213 7816.581 -835383.43814 27176.234 -861713.93815 22063.230 -905726.12516 0.000 -927165.18817 29514.479 -665602.87518 53419.637 -747340.43819 34876.832 -839806.81220 0.000 -881219.12521 29580.273 -416288.71922 52944.918 -632601.12523 17504.195 -803765.68824 0.000 -859481.93825 0.000 0.00026 -120102.820 -583505.37527 -76202.375 -787347.18828 0.000 -829170.812E= 1SX= -1489.530 SY=-101489.383 TAU= -1489.531 S1= -1467.348 S2=-101511.562 CET= 179.147 E= 2SX= -1475.844 SY=-100654.875 TAU= -1790.500 S1= -1443.531 S2=-100687.188 CET= 178.966 E= 3SX= -7021.670 SY=-101597.672 TAU= -3741.688 S1= -6873.875 S2=-101745.469 CET= 177.738 E= 4SX= -8067.500 SY= -98528.750 TAU= -4459.156 S1= -7848.227 S2= -98748.023 CET= 177.185 E= 5SX= -13143.328 SY= -99391.750 TAU= -1662.500 S1= -13111.293 S2= -99423.781 CET= 178.896 E= 6SX= -14652.781 SY= -98337.500 TAU= -1501.062S1= -14625.867 S2= -98364.414 CET= 178.973 E= 7SX= -2923.122 SY=-109168.297 TAU= -5888.469 S1= -2597.762 S2=-109493.656 CET= 176.837 E= 8SX= -716.078 SY=-103681.562 TAU= -8617.406 S1= 0.148 S2=-104397.789 CET= 175.249 E= 9SX= -9188.316 SY=-105121.867 TAU= -9771.594 S1= -8203.125 S2=-106107.062 CET= 174.243 E= 10SX= -12285.000 SY= -95180.250 TAU= -12199.594 S1= -10526.887 S2= -96938.359 CET= 171.799 E= 11SX= -14170.516 SY= -95500.750 TAU= -5489.531 S1= -13801.664 S2= -95869.602 CET= 176.156 E= 12SX= -22797.406 SY= -91347.000 TAU= -3902.844 S1= -22575.914 S2= -91568.492 CET= 176.752 E= 13SX= -5104.269 SY=-129494.438 TAU= -11708.750 S1= -4011.727 S2=-130586.977 CET= 174.669 E= 14SX= 969.672 SY=-108176.375 TAU= -21424.750 S1= 5024.582 S2=-112231.281 CET= 169.283 E= 15SX= -14954.572 SY=-110883.469 TAU= -18383.531 S1= -11552.273 S2=-114285.766 CET= 169.515 E= 16SX= -19890.141 SY= -86924.312 TAU= -25131.188 S1= -11514.844 S2= -95299.609 CET= 161.569 E= 17SX= -22109.688 SY= -87301.625 TAU= -10225.406 S1= -20543.453 S2= -88867.859 CET= 171.292 E= 18SX= -35190.453 SY= -77219.000 TAU= -9162.000 S1= -33280.023 S2= -79129.430 CET= 168.222 E= 19SX= -9785.850 SY=-171444.172 TAU= -20524.969 S1= -7220.594 S2=-174009.422 CET= 172.876 E= 20SX= 4594.438 SY=-113592.375 TAU= -46145.688 S1= 20477.398 S2=-129475.336 CET= 161.007 E= 21SX= -25287.307 SY=-118672.312 TAU= -30023.750 S1= -16467.512 S2=-127492.109 CET= 163.629 E= 22SX= -30634.422 SY= -71127.188 TAU= -44991.469 S1= -1543.715 S2=-100217.891 CET= 147.114 E= 23SX= -34259.609 SY= -71743.438 TAU= -14637.906 S1= -29220.699 S2= -76782.344 CET= 161.005 E= 24SX= -43958.047 SY= -53418.938 TAU= -17697.562 S1= -30369.627 S2= -67007.359 CET= 142.482 E= 25SX= -19028.160 SY=-252549.000 TAU= -34958.688 S1= -13907.055 S2=-257670.094 CET= 171.666 E= 26SX= 3973.812 SY=-114063.750 TAU= -92238.344 S1= 54459.047 S2=-164548.984 CET= 151.307 E= 27SX= -39180.809 SY=-121400.055 TAU= -39312.688 S1= -23409.074 S2=-137171.781 CET= 158.140 E= 28SX= -42804.766 SY= -43317.938 TAU= -65723.062 S1= 22662.211 S2=-108784.914 CET= 135.112 E= 29SX= -42224.094 SY= -43219.188 TAU= -10273.375 S1= -32436.225 S2= -53007.055 CET= 136.386 E= 30SX= -21830.422 SY= -25448.312 TAU= -23810.344 S1= 239.594 S2= -47518.328 CET= 137.172 E= 31SX= -48815.199 SY=-424587.344 TAU= -79800.078 S1= -32570.844 S2=-440831.688 CET= 168.494 E= 32SX=-132271.750 SY= -71582.000 TAU=-175409.250 S1= 76087.781 S2=-279941.531 CET= 130.093 E= 33SX= -45090.102 SY= -56761.105 TAU= 804.781 S1= -45034.867 S2= -56816.336 CET= 3.926 E= 34SX= 42332.711 SY= -9221.938 TAU= -47066.344 S1= 70218.328 S2= -37107.555 CET= 149.354 E= 35SX= -20899.344 SY= -19971.375 TAU= 16235.219 S1= -4193.512 S2= -36677.207 CET= 45.819E= 36SX= 73163.914 SY= -17873.250 TAU= -17873.344 S1= 76547.250 S2= -21256.586 CET= 169.281。
有限元法编程
有限元法编程概述有限元法编程是一种在工程领域广泛应用的数值计算方法，它通过将复杂的连续体问题离散为有限个小单元，然后通过对这些小单元进行数值计算，得到整体问题的近似解。
本文将详细介绍有限元法编程的基本原理、步骤以及在实际工程问题中的应用。
基本原理有限元法编程的基本原理是将连续体分割为若干个小单元，每个小单元称为有限元。
这些有限元通过节点相连形成一个离散网格，然后通过对每个有限元的本构关系进行数值计算，得到整个连续体的力学行为。
有限元法编程的基本步骤如下：1.网格生成：通过一定的方法将连续体分割为有限元网格；2.选取插值函数：由于力学场在每个有限元上是未知的，为了对力学场进行数值计算，需要对其进行插值。
常用的插值函数有线性插值和二次插值等；3.设置本构关系：选择适当的本构关系来描述材料的力学性能；4.形成方程组：通过应力-应变关系和边界条件，列出表示力学问题的线性方程组；5.求解方程组：利用数值方法求解线性方程组，得到力学场的数值解；6.后处理：根据实际问题的需要，对数值解进行进一步的处理和分析。
网格生成网格生成是有限元法编程的第一步，它的目的是将连续的问题离散为有限个小单元。
常用的网格生成方法有四边形单元法、三角形剖分法和网格生成软件等。
其中，四边形单元法适用于二维问题，三角形剖分法适用于平面或曲面问题，而网格生成软件则可以用于生成复杂的三维网格。
选取插值函数选取插值函数是为了对力学场进行数值计算。
常用的插值函数有线性插值和二次插值。
线性插值函数适用于简单的问题，而二次插值函数则适用于更复杂的情况。
通过对节点位置和节点上的数值进行插值，可以在整个有限元中得到力学场的数值近似解。
设置本构关系本构关系是材料力学性能的数学表达式。
根据不同的材料特性和力学问题，可以选择适当的本构关系。
常见的本构关系有线弹性模型、非线性弹性模型和塑性模型等。
形成方程组通过应力-应变关系和边界条件，可以得到力学问题的线性方程组。
有限元实例教程
• 9.修改【优化设置】，对话框中的【最大迭 代次数】为10，其它参数值均为默认。 • （2）优化求解及其结果查看 • 1）右键单击【Setup1】节点，在快捷菜单 中单击【求解】命令，如下图 即可提交作 业解算。 • 2）稍等完成计算作业，并切换到Microsoft Excel工作程序显示优化结果，其中该表包 括【Optimization】,【Objective】,和 【Link的表达式】3个工作表格，其中 【Optimization】工作表主要显示设计目标
• 3）在【优化设置】对话框中单击【定义目标】按 钮，弹出图框如下图所示，默认目标对象为【模 型】，在图形窗口中选中连杆机构，默认【类型】 为【重量】，默认目标为【最小化】， • 4）在【优化设置】对话框中单击【定义约束】， 弹出【约束】对话框，默认类型为【位移】，选 择【平移X轴】，默认【限制类型】为【上限】， 【限制值】为0.04。 • 5）继续进行【约束】定义,选择【3D对象】，在 【类型】中设为【应力】，应用于【面】，单击 【VonMises】,单击连杆模型小孔内表面的右侧 面，默认类型为【上限】，【限制值】为225， 【下限】系统默认为0
4）在图形窗口中单击创建的基准平面，单击【确定】按钮即可将连杆小孔分割 为两个半圆柱面，为施加单侧力载荷提供了条件。 5）右键单击【仿真导航器】窗口中的【fem-link-.prt】节点，在弹出的快捷菜单 中单击【显示FEM】→【fem-link。Femt】命令，返回到有限元模型环境。 6）单击工具栏中的【更新有限元模型】图标，等待模型网络重新划分完成后， 右键单击【仿真导航器】窗口中的【fem-link.femt】节点，在弹出的快捷菜 单中【显示仿真】→【fem-link.sim器】图标，弹出对话框，在【单元族】选项中选 择【3D】，在【集合类型】选项中选择【实体】，在【物理属性】的【solid property】选项中选取上述设置的【PSOILID1】，默认【网格补集器】的名 称为【solid（1）】，单击【确定】。 6）单击工具栏中的【3D四面体网络】图标右侧的小三角形，弹出对话框，在窗 口中选择连杆三维模型，默认【单元属性】的【类型】为【CTETRA （10）】，单击【单元大小】右侧的【自动单元大小】图标，对话框中出现 【5.78】，手动修改为【2】，【目标补集器】中【meshcollector】选项为上 述设置的【solid（1）】，其他选项均为默认值，单击【确定】按钮。 7）单击工具栏中的【有限元模型检查】，弹出【模型检查】对话框，单击【应 用】按钮，在弹出的【信息】中出现【Number Fauled】信息，发现模型正 常，没有出现划分失败的网络。
有限元C++编程实践
基于有限元算法的编程实践学号：2011043010031姓名：廖校毅电子科技大学物理电子学院【摘要】本文就有限元算法在电磁场分析中的应用展开研究与实践，从电磁场的Maxwell方程出发,根据电磁场的边值问题及变分公式建立了有限元方程组。
具体在实践中，将这些知识运用到C++语言和Matlab中，并将这两种语言有机结合，编程并实现二维FEM。
程序最后通过计算含两种介质的电位槽电位分布来验证其正确性。
关键词: 有限元变分法C++ MatlabThe Programming Practice Based on The Finite Element Algorithm Student ID：2011043010031Name：Liao Xiaoyi University of Electronic Science and technology &School of PhysicalElectronicsAbstract In this paper, we take the application of finite element method in electromagnetic field analysis into research and practice. Starting from Maxwell equations of electromagnetic field, the electromagnetic field boundary value problems and variational formula established the finite element equations. In specific practice, this knowledge will be applied to the C++ language and Matlab, and the organic combination of two languages, programming and implementation of two-dimensional FEM. Finally, through the program to verify the validity of the calculation of potential distribution in channel potential containing two kinds of medium.Key words The finite element method The variational method C++ Matlab一、前言在数学中，有限元法（FEM，Finite Element Method）是一种为求得偏微分方程边值问题近似解的数值技术。
有限单元法程序代码
弹性力学平面问题有限元程序#define NODE_NUM 300#define ELE_NUM 200#define NODE2 NODE_NUM*2#include<stdlib.h>#include<stdio.h>#include<math.h>int nj,ne,nm,nz,npj,jm[ELE_NUM][4],nzc[60],nj2,bw,lxm;float ea,e,um,oul,te,xy[NODE_NUM][4],EA[5],p[NODE2],pj,F[ELE_NUM],etm[5][5]; float ek[7][7];double **K,*f;int BAND2(int n,int m, double **K,double *f);double **TwoArrayAlloc(int,int);void TwoArrayFree(double **);FILE *outf;/*-------------------------------------------*/double **TwoArrayAlloc(int r,int c){double *x,**y;int n;x=(double*)calloc(r*c,sizeof(double));y=(double**)calloc(r,sizeof(double*));for(n=0;n<=r-1;++n)y[n]=&x[c*n];}/*-------------------------------------------*/void TwoArrayFree(double **x){free(x[0]);free(x);}/*------------------------------------------*/int BAND2(int n,int m, double **K,double *f) {int i,j,t,ij,ji,it,jt,tm,m1;double s,w;m1=m+1;for(i=1;i<=n;i++){if(K[i-1][m1-1]<=0) return(0);w=0.0;if(i>m1) tm=i-m;else tm=1;for(j=tm;j<=i;j++){s=0.0;ij=j-i+m1;for(t=tm;t<=j-1;t++){jt=t-j+m1;s=s+K[i-1][it-1]*K[j-1][jt-1]/K[t-1][m1-1]; }K[i-1][ij-1]=K[i-1][ij-1]-s;if(j==i)f[i-1]=f[i-1]-w;elsew=w+K[i-1][ij-1]*f[j-1]/K[j-1][m1-1];}}for(i=n;i>=1;i--){s=0.0;if(i>n-m1) tm=n; else tm=i+m;for(j=i+1;j<=tm;j++){ji=i-j+m1;s=s+K[j-1][ji-1]*f[j-1];}f[i-1]=(f[i-1]-s)/K[i-1][m1-1];}return 1;}/*-----------------------------------------*/void input(){int jj,j,i,nh,nl;float dx,dy;FILE *infile;char name1[30],name2[30];R:printf("please enter data-filename\n");scanf("%s",name1);if((infile=fopen(name1,"r"))==NULL){ printf("the data-file not exit!");goto R;}printf("Please enter result-filename\n");scanf("%s",name2);outf=fopen(name2,"w");fscanf(infile,"%d%d%d%d%d%d5d",&nj,&ne,&nm,&nz,&npj,&lxm); for(i=1;i<=nj;i++){for(j=1;j<3;j++)fscanf(infile,"%f",&xy[i][j]);}for(i=1;i<=ne;i++){for(j=0;j<4;j++)fscanf(infile,"%d",&jm[i][j]);}for(i=1;i<=nm;i++){for(j=1;j<=4;j++)fscanf(infile,"%f",&etm[i][j]);}for(i=1;i<=nz;i++)fscanf(infile,"%d",&nzc[i]);nj2=nj*2;for(i=0;i<nj2;i++) p[i]=0.0;for(i=0;i<npj;i++){fscanf(infile,"%d%f",&jj,&pj);p[jj]=pj;}fprintf(outf,"The Num. Of Nodes: %3d\n",nj);fprintf(outf,"The Num. Of Elements.: %3d\n",ne);fprintf(outf,"The Num. Of Type Of Section Characteristic: %3d\n",nm); fprintf(outf,"The Num. Of Restriction: %3d\n",nz);fprintf(outf,"The Num. Of Nodal Loads: %3d\n",npj);fprintf(outf,"LXM= %d\n",lxm);fprintf(outf,"Coordinates x and y Of Nodes:\n");fprintf(outf," Node x y\n");for(i=1;i<=nj;i++){fprintf(outf,"%10d%10.2f%10.2f\n",i,xy[i][1],xy[i][2]);}fprintf(outf,"The Nodes Num. Of Mem:\n");fprintf(outf," Mem. Type Left Right\n");for(i=1;i<=ne;i++){fprintf(outf,"%10d%10d%10d%10d%10d\n",i,jm[i][0],jm[i][1],jm[i][2],jm[i][3]);}for(i=1;i<=nm;i++){fprintf(outf,"nm=%d E=%f,UM=%f,OUL=%f,T=%f\n",i,etm[i][1],etm[i][2],etm[i][3],etm[i][ 4]);}fprintf(outf,"The Position Of Retriction:");for(i=1;i<=nz;i++){fprintf(outf,"%d ",nzc[i]);}fprintf(outf,"\n===========loads of nodes==============\n");for(i=1;i<=nj2;i++){fprintf(outf,"%12d %12.2f\n",i,p[i]);}if(lxm>0){for(i=1;i<=nm;i++){etm[i][1]=etm[i][1]/(1-etm[i][2]*etm[i][2]); }}fclose(infile);}/*------------计算半带宽函数----------------*/ int bwidth(){int ie,i,j,min,max,jj,ib,lk,m,m1,m2,t,nn[8]; ib=0;for(ie=1;ie<=ne;ie++){m=abs(jm[ie][1]-jm[ie][2]);m1=abs(jm[ie][1]-jm[ie][3]);m2=abs(jm[ie][2]-jm[ie][3]);if(m>ib) ib=m;if(m1>ib) ib=m1;if(m2>ib)ib=m2;}bw=2*ib+2;fprintf(outf," bw=%d\n ",bw);return bw;}/*--------------形成单元刚度矩阵函数-----------*/void stiff(int ie){int i,j,k,m,n,i1,j1,i2,j2,ii,jj;float cc,bb,bc,cb,a1,u1,c[4],b[4];i=jm[ie][1];j=jm[ie][2];m=jm[ie][3]; /*取出单元节点号*/n=jm[ie][0];e=etm[n][1];um=etm[n][2];oul=etm[n][3];te=etm[n][4]; /*取出材料参数*/c[3]=xy[j][1]-xy[i][1];b[3]=xy[i][2]-xy[j][2];c[2]=xy[i][1]-xy[m][1];b[2]=xy[m][2]-xy[i][2];c[1]=-c[2]-c[3];b[1]=-b[2]-b[3]; /*计算[B]阵中的bi,ci,bj,cj,bm,cm*/ ea=(b[2]*c[3]-b[3]*c[2])/2.0;/*计算单元面积*/a1=e*te/4.0/(ea-um*um*ea);u1=(1.0-um)/2.0;for(ii=1;ii<=3;ii++) /*按2子块循环求出单刚计算单元面积*/ {for(jj=1;jj<=3;jj++){bb=b[ii]*b[jj];cc=c[ii]*c[jj];cb=c[ii]*b[jj];bc=b[ii]*c[jj];i2=2*ii;j2=2*jj;i1=i2-1;j1=j2-1;ek[i1][j1]=a1*(bb+u1*cc);ek[i1][j2]=a1*(um*bc+u1*cb);ek[i2][j1]=a1*(um*cb+u1*bc);ek[i2][j2]=a1*(cc+u1*bb);}}}/*-------集成总刚函数，存在刚下三角--------------------------*/int ekzk(int ie){int i1,j1,i,j,i2,j2,ii,jj,ji;for(i1=1;i1<=3;i1++){for(i2=1;i2<3;i2++){i=2*(i1-1)+i2;ii=2*(jm[ie][i1]-1)+i2;for(j1=1;j1<=3;j1++){for(j2=1;j2<3;j2++){j=2*(j1-1)+j2;jj=2*(jm[ie][j1]-1)+j2;ji=bw+jj-ii;if(ji<=bw)K[ii-1][ji-1]=K[ii-1][ji-1]+ek[i][j];}}}}}/*--------用0，1置换法进行约束处理-----------------------------*/ void restrict(){int i,j,i1,j1,ii;for(i=1;i<=nz;i++){j=nzc[i]-1;for(i1=0;i1<bw-1;i1++){K[j][i1]=00.00;ii=j+i1+1;if(ii<nj2)K[ii][bw-i1-2]=00.00;}K[j][bw-1]=1.00;f[j]=0;}}/*----------计算单元应力函数-----------------*/void force(){int i,j,ie,m,n,j1,j2,ii,jj;float aa,pyl,ryl,w[7],c[4],b[4],s1[4][7],yl[4],yl1[4];fprintf(outf,"==========stresses of elements ==========\n"); fprintf(outf,"NO SX SY ST S1 S2 CETA\n\n"); for(ie=1;ie<=ne;ie++){i=jm[ie][1];j=jm[ie][2];m=jm[ie][3];n=jm[ie][0];e=etm[n][1];um=etm[n][2];oul=etm[n][3];te=etm[n][4];c[3]=xy[j][1]-xy[i][1];b[3]=xy[i][2]-xy[j][2];c[2]=xy[i][1]-xy[m][1];b[2]=xy[m][2]-xy[i][2];c[1]=-c[2]-c[3];b[1]=-b[2]-b[3];ea=(b[2]*c[3]-b[3]*c[2])/2.0; aa=e/(2.0*ea*(1-um*um)); w[1]=f[2*i-2];w[2]=f[2*i-1];w[3]=f[2*j-2];w[4]=f[2*j-1];w[5]=f[2*m-2];w[6]=f[2*m-1];for(ii=1;ii<=3;ii++){j1=2*ii-1;j2=j1+1;s1[1][j1]=b[ii]*aa;s1[2][j1]=s1[1][j1]*um;s1[3][j2]=s1[1][j1]*(1-um)/2.0;s1[2][j2]=c[ii]*aa;s1[1][j2]=s1[2][j2]*um;s1[3][j1]=s1[2][j2]*(1-um)/2.0;}for(ii=1;ii<=3;ii++){yl[ii]=0.0;for(jj=1;jj<=6;jj++){yl[ii]=yl[ii]+s1[ii][jj]*w[jj];}}pyl=(yl[1]+yl[2])/2.0;ryl=sqrt(((yl[1]-yl[2])/2.0)*((yl[1]-yl[2])/2.0)+yl[3]*yl[3]); yl1[1]=pyl+ryl;yl1[2]=pyl-ryl;aa=abs(yl[2]-yl1[2]);if(abs(yl[2]-yl1[2])<6.1e-2){yl1[3]=0.0;}else yl1[3]=90.0-57.29578*atan(yl[3]/(yl[2]-yl1[2]));fprintf(outf,"%2d%10.2f%10.2f%10.2f%10.2f%10.2f%10.2f\n",ie,yl[1],yl[2],yl[3],yl1[1],yl1[2 ],yl1[3]);}}/*---------------------------------------------------------*/void main(){int ie,i,j,ii,jj,i1,j1;float p1,p2,p3;input();bw=bwidth();K=TwoArrayAlloc(nj2,bw);f=(double*)calloc(nj2,sizeof(double));if(f==NULL) exit(1);for(i=0;i<nj2;i++) f[i]=p[i+1];for(i=0;i<nj2;i++)for(j=0;j<bw;j++) K[i][j]=0;for(ie=1;ie<=ne;ie++){stiff(ie);ekzk(ie);}restrict();if(!BAND2(nj2,bw-1,K,f)) {printf("Failed !\n"); exit(0);}fprintf(outf,"=============== displacements of nodes ================\n");fprintf(outf,"No x y\n");for(ii=0;ii<nj;ii++){p1=f[2*ii];p2=f[2*ii+1];fprintf(outf,"%2d %14.5f %14.5f\n",ii+1,p1,p2);}force();}平面桁架静力分析有限元程序#define NODE_NUM 300#define ELE_NUM 200#include<stdlib.h>#include<stdio.h>#include<math.h>int nj,ne,nm,nz,npj,jm[ELE_NUM][3],nzc[60],nj2,bw;float xy[NODE_NUM][4],EA[5],p[NODE_NUM],pj,F[ELE_NUM];float ek[5][5];double **K,*f;int BAND2(int n,int m,double **K,double *f);double **TwoArrayAlloc(int,int);void TwoArrayFree(double**);FILE *outfile;/*-----------------------------*/double **TwoArrayAlloc(int r,int c){double *x,**y;int n;x=(double*)calloc(r*c,sizeof(double));y=(double**)calloc(r,sizeof(double*));for(n=0;n<=r-1;++n)y[n]=&x[c*n];return(y);}/*--------------------------------*/void TwoArrayFree(double **x){free(x[0]);free(x);}/*--------------------------------*/int BAND2(int n,int m,double **K,double *f) /*解方程函数*/ {int i,j,t,ij,ji,it,jt,tm,m1;double s,w;m1=m+1;for(i=1;i<=n;i++){if(K[i-1][m1-1]<=0) return(0);w=0.0;if(i>m1) tm=i-m;else tm=1;for(j=tm;j<=i;j++){s=0.0;ij=j-i+m1;for(t=tm;t<=j-1;t++){it=t-i+m1;jt=t-j+m1;s=s+K[i-1][it-1]*K[j-1][jt-1]/K[t-1][m1-1];}K[i-1][ij-1]=K[i-1][ij-1]-s;if(j==i)f[i-1]=f[i-1]-w;elsew=w+K[i-1][ij-1]*f[j-1]/K[j-1][m1-1];}}for(i=n;i>=1;i--){s=0.0;if(i>n-m1) tm=n; else tm=i+m;for(j=i+1;j<=tm;j++){ji=i-j+m1;s=s+K[j-1][ji-1]*f[j-1];}f[i-1]=(f[i-1]-s)/K[i-1][m1-1];}return 1;}/*--------------------------------------------*/void input()/*输入函数*/{int jj,j,i,nh,nl;float dx,dy;FILE *infile;infile=fopen("dat.txt","r");outfile=fopen("out.txt","w");fscanf(infile,"%d%d%d%d%d",&nj,&ne,&nm,&nz,&npj); /*读入控制数据*/fprintf(outfile,"The Num. Of Nodes: %3d\n",nj);/*输出节点数*/fprintf(outfile,"The Num. Of Mem.: %3d\n",ne);/*输出单元数*/fprintf(outfile,"The Num. Of Type Of Section Characteristic: %3d\n",nm);/*输出材料类型数*/fprintf(outfile,"The Num. Of Restriction: %3d\n",nz);/*输出荷载数*/fprintf(outfile,"The Num. Of Nodal Loads: %3d\n",npj);/*输出荷载数*/for(i=1;i<=nj;i++){for(j=1;j<3;j++)fscanf(infile,"%f",&xy[i][j]);/*读入节点坐标*/}fprintf(outfile,"Coordinates x and y Of Nodes:\n");fprintf(outfile," Node x y\n");for(i=1;i<=nj;i++){fprintf(outfile,"%10d%10.2f%10.2f\n",i,xy[i][1],xy[i][2]);/*输出节点坐标*/}for(i=1;i<=ne;i++){for(j=0;j<3;j++)fscanf(infile,"%d",&jm[i][j]);/*读入单元信息（材料类型、单元左右节点码）*/}fprintf(outfile,"The Nodes Num. Of Mem.:\n");fprintf(outfile," Mem. Type Left Right\n");for(i=1;i<=ne;i++){fprintf(outfile,"%10d%10d%10d%10d\n",i,jm[i][0],jm[i][1],jm[i][2]);/*输出单元信息*/}for(i=1;i<=nm;i++) fscanf(infile,"%f",&EA[i]);/*读入材料刚度*/for(i=1;i<=nm;i++){fprintf(outfile,"nm=%d EA=%f\n",i,EA[i]);/*输出材料刚度*/}for(i=1;i<=nz;i++) fscanf(infile,"%d",&nzc[i]);/*读入约束位置*/fprintf(outfile,"The Position Of Restriction: ");for(i=1;i<=nz;i++){fprintf(outfile,"%d",nzc[i]);/*输出约束位置*/}fprintf(outfile,"\n");nj2=nj*2;/*形成荷载*/for(i=0;i<nj2;i++) p[i]=0.0;for(i=0;i<npj;i++){fscanf(infile,"%d%f",&jj,&pj);/*读入荷载位置、荷载大小*/p[jj]=pj;}fclose(infile);}/*-------------------------------------------------------------------------*/int bwidth()/*求半带宽函数*/{int ie,i,j,min,max,jj,ib,lk,m,nn[8];ib=0;for(ie=1;ie<=ne;ie++){m=abs(jm[ie][1]-jm[ie][2]);if(m>ib) ib=m;/*找出单元左右节点码之差的最大值*/}return 2*ib+2;/*返回半带宽*/}/*-------------------------------------------------------------------------------*/ void stiff(int ie) /*形成单刚函数*/{int i,j,k,m,n;float dx,dy,dz,l,cx,cy,cz,ea;i=jm[ie][1];/*单元左节点号*/j=jm[ie][2];/*单元右节点号*/m=jm[ie][0];/*单元材料类型号*/dx=xy[j][1]-xy[i][1];/*单元左右节点横坐标之差*/dy=xy[j][2]-xy[i][2];/*单元左右节点纵坐标之差*/l=sqrt(dx*dx+dy*dy);/*单元长度*/cx=dx/l;/*求余弦*/cy=dy/l;/*求正弦*/ea=EA[m]/l;/*fprintf(outfile,"%d%10.2f\n",ie,l);*/ek[1][1]=ek[3][3]=cx*cx*ea;/*求单刚*/ek[2][2]=ek[4][4]=cy*cy*ea;ek[2][1]=ek[1][2]=ek[3][4]=ek[4][3]=cx*cy*ea;ek[4][1]=ek[1][4]=ek[3][2]=ek[2][3]=-cx*cy*ea;ek[1][3]=ek[3][1]=-cx*cx*ea;ek[2][4]=ek[4][2]=-cy*cy*ea;}/*-----------------------------------------------------------------------------*/ int ekzk(int ie)/*集成总刚度*/{int i1,j1,i,j,i2,j2,ii,jj,ji;for(i1=1;i1<=2;i1++){for(i2=1;i2<=2;i2++){i=2*(i1-1)+i2;ii=2*(jm[ie][i1]-1)+i2;for(j1=1;j1<=2;j1++){for(j2=1;j2<=2;j2++){j=2*(j1-1)+j2;jj=2*(jm[ie][j1]-1)+j2;ji=bw+jj-ii;if(ji<=bw) K[ii-1][ji-1]=K[ii-1][ji-1]+ek[i][j];}}}}}/*---------------------------------------------------------------------------*/ void force()/*求单元轴力函数*/{int i,j,ie,m;float dx,dy,dz,l,cx,cy,cz,ea,w[7];fprintf(outfile,"============dyzl==================\n"); for(ie=1;ie<=ne;ie++){i=jm[ie][1];j=jm[ie][2];m=jm[ie][0];w[1]=f[2*i-2];w[2]=f[2*i-1];w[3]=f[2*j-2];w[4]=f[2*j-1];dx=xy[j][1]-xy[i][1];dy=xy[j][2]-xy[i][2];l=sqrt(dx*dx+dy*dy);cx=dx/l;cy=dy/l;ea=EA[m]/l;dx=w[3]-w[1];/*左右节点水平位移之差*/dy=w[4]-w[2];/*左右节点竖直位移之差*/l=ea*(cx*dx+cy*dy);/*求单元轴力*/F[ie]=1;fprintf(outfile,"%d%10.4f\n",ie,l);/*输出单元号、单元轴力*/}}/*--------------------------------------------------------------------------------------*/void main(){int ie,i,j,ii,jj,i1,j1;float p1,p2,p3;input();bw=bwidth();K=TwoArrayAlloc(nj2,bw);f=(double*)calloc(nj2,sizeof(double));if(f==NULL) exit(1);for(i=0;i<nj2;i++) f[i]=p[i+1];for(ie=1;ie<=ne;ie++){stiff(ie);ekzk(ie);}for(i=1;i<=nz;i++)/*约束处理（后处理的0、1置换法）*/{j=nzc[i]-1;for(i1=0;i1<bw-1;i1++){K[j][i1]=00.00;/*将有约束处的行置0*/ii=j+i1+1;if(ii<nj2) K[ii][bw-i1-2]=00.00;/*将有约束处的列置0*/}K[j][bw-1]=1.00;/*将对角线元素置1*/f[j]=0;}/*fprintf(outfile,"========================zk======================\n");for(i=0;i<nj2;i++){for(j=0;j<bw;j++){fprintf(outfile,"%6.2f",K[i][j]);}fprintf(outfile,"\n");}*/fprintf(outfile,"============jdhz=================\n");for(i=0;i<nj2;i++){fprintf(outfile,"%d %6.2f\n",i+1,f[i]);/*输出节点力*/}if(!BAND2(nj2,bw-1,K,f)){printf("Failed!\n"); exit(0);}fprintf(outfile,"=========jdwy==========\n");fprintf(outfile,"xy\n");for(ii=0;ii<nj;ii++){p1=f[2*ii];p2=f[2*ii+1];fprintf(outfile,"%2d %14.5f %14.5f\n",ii+1,p1,p2);/*输出节点位移*/ }force();fclose(outfile);}平面刚架静力分析有限元程序#include<io.h>#include"string.h"#include"stdlib.h"#include"stdio.h"#include"math.h"void hbw();void sncs(int nel);void fix(int np);void trmat();void fis(int nel);void fpj();void force();void stiff();void addsm();void restr();void matmul();void soleq();void outdis();float sm[6][6],tsm[300][30],res[60][2],pj[300],t[6][6],d[6][6],r[300], fo[6],foj[6],pf[200][4],x[100],y[100],ae[10][3],sl,sn,cs,eal,eil;int nj,ne,nt,nr,npj,npf,nn,mbw,jel[100][2],nae[100],is[6];FILE *infile,*outfile;/************主函数**************/void main(){char name1[30],name2[30];int i,j,nel,np;printf("please enter data-filename\n");scanf("%s",name1);printf("please enter result-filename\n");scanf("%s",name2);if((infile=fopen(name1,"r"))!=NULL){fscanf(infile,"%d%d%d%d%d%d",&nj,&ne,&nt,&nr,&npj,&npf);for(i=0;i<ne;i++) fscanf(infile,"%d%d%d",&jel[i][0],&jel[i][1],&nae[i]); for(i=0;i<nj;i++) fscanf(infile,"%f%f",&x[i],&y[i]);for(i=0;i<nt;i++) fscanf(infile,"%f%f%f",&ae[i][0],&ae[i][1],&ae[i][2]); for(i=0;i<nr;i++) fscanf(infile,"%f%f",&res[i][0],&res[i][1]);}else{printf("the data-file not exit!");exit(1);}nn=3*nj;outfile=fopen(name2,"w");if(outfile==NULL){printf("the result-file not exist!");exit(1);}fprintf(outfile,"statis analysis of plane frame\n");fprintf(outfile,"input data\n");fprintf(outfile,"************\n");fprintf(outfile,"control data\n");fprintf(outfile,"the num.of node:%3d\n",nj);fprintf(outfile,"the num.of mem:%3d\n",ne);fprintf(outfile,"the num.of type of section characteristic:%3d\n",nt); fprintf(outfile,"the num.of restricted degrees of freedom:%3d\n",nr); fprintf(outfile,"the num.of nodal loads:%3d\n",npj);fprintf(outfile,"the num.of non-nodal loads:%3d\n",npf);fprintf(outfile,"the num.of nodal degrees of freedom:%3d\n",nn); fprintf(outfile,"information of mem.\n");fprintf(outfile,"mem. start node end node type\n");for(i=0;i<ne;i++)fprintf(outfile,"%5d%10d%10d%10d\n",i+1,jel[i][0],jel[i][1],nae[i]); fprintf(outfile,"coordinates x and y of nodes\n");fprintf(outfile,"node x y\n");for(i=0;i<nj;i++)/*第二次输入*/fprintf(outfile,"%10d%10.2f%10.2f\n",i+1,x[i],y[i]);fprintf(outfile,"information of cross section each mem.\n");fprintf(outfile,"typeaeramoment of interia elastic modulus\n");for(i=0;i<nt;i++)fprintf(outfile,"%8d%15.5f%15.5f%20.2f\n",i+1,ae[i][0],ae[i][1],ae[i][2]); fprintf(outfile,"information restriction\n");fprintf(outfile,"rester.-norestr.-disp.-norestr.-disp.\n");for(i=0;i<nr;i++)fprintf(outfile,"%10d%19.3f%19.3f\n",i+1,res[i][0],res[i][1]);hbw();for(i=0;i<nn;i++) pj[i]=0.0;if(npj==0) goto aa;fprintf(outfile,"nodal loads\n");fprintf(outfile,"nodepxpyzm\n");for(i=0;i<npj;i++){ int no;float px,py,zm;fscanf(infile,"%d%f%f%f",&no,&px,&py,&zm);fprintf(outfile,"%10d%10.2f%10.2f%10.2f\n",no,px,py,zm);pj[3*no-3]=px; pj[3*no-2]=py; pj[3*no-1]=zm;}aa:for(i=0;i<nn;i++) r[i]=0.0;for(i=0;i<nr;i++){int ni;ni=floor(res[i][0]+0.1)-1;if(pj[ni]!=0.0)r[ni]=r[ni]-pj[ni];}if(npf==0) goto bb;fprintf(outfile,"non-nodal loads\n");fprintf(outfile,"loads leng.supp.to load mem type\n");for(i=0;i<npf;i++){fscanf(infile,"%f%f%f%f%f",&pf[i][0],&pf[i][1],&pf[i][2],&pf[i][3]);fprintf(outfile,"%15.2f%15.2f%15.2f%15.2f\n",pf[i][0],pf[i][1],pf[i][2],pf[i][3]); }for(np=0;np<npf;np++){nel=floor(pf[np][2]+0.1);sncs(nel-1);fix(np);trmat();fis(nel-1);fpj();}bb:for(i=0;i<nn;i++)for(j=0;j<mbw;j++) tsm[i][j]=0.0; for(nel=0;nel<ne;nel++){sncs(nel);trmat();stiff();matmul();fis(nel);addsm();}restr();soleq();outdis();force();}/*===================求最大半带宽========================*/ void hbw(){int nel;mbw=0;for(nel=0;nel<ne;nel++){int ma=abs(jel[nel][0]-jel[nel][1]);if(mbw<ma) mbw=ma;}mbw=3*(mbw+1);fprintf(outfile,"Half_Bandwidth Mbw:%5d\n",mbw);}/*矩阵相乘*/void matmul(){int i,j,k;for(i=0;i<6;i++)for(j=0;j<6;j++) d[i][j]=0.0;for(i=0;i<6;i++)for(j=0;j<6;j++)for(k=0;k<6;k++) d[i][j]=d[i][j]+t[k][i]*sm[k][j];for(i=0;i<6;i++)for(j=0;j<6;j++) sm[i][j]=0.0;for(i=0;i<6;i++)for(j=0;j<6;j++)for(k=0;k<6;k++)sm[i][j]=sm[i][j]+d[i][k]*t[k][j];}/*求单元常数*/void sncs(int nel){int ii,jj,k;float xi,yi,xj,yj;ii=jel[nel][0];jj=jel[nel][1];xi=x[ii-1];xj=x[jj-1];yi=y[ii-1];yj=y[jj-1];sl=sqrt((xj-xi)*(xj-xi)+(yj-yi)*(yj-yi)); sn=(yj-yi)/sl;cs=(xj-xi)/sl;k=nae[nel];eal=ae[k-1][0]*ae[k-1][2]/sl;eil=ae[k-1][1]*ae[k-1][2]/sl;}/*求单元坐标转换矩阵*/void trmat(){int i,j;for(i=0;i<6;i++)for(j=0;j<6;j++)t[i][j]=0.0;t[0][0]=cs;t[0][1]=sn;t[1][0]=-sn;t[1][1]=cs;t[2][2]=1.0;for(i=0;i<3;i++)for(j=0;j<3;j++) t[i+3][j+3]=t[i][j];}/*求单元固端力*/void fix(int np){float w,c,c1,c2,cc,cc1,cc2;int i,im;w=pf[np][0];c=pf[np][1];c1=c/sl; c2=c1*c1; cc=sl-c;cc1=cc/sl; cc2=cc1*cc1;for(i=0;i<6;i++) fo[i]=0.0;im=floor(pf[np][3]+0.1);if(im==1){fo[1]=-w*c*(1.0-c2+c2*c1/2.0);fo[2]=-w*c*c*(6.0-8.0*c1+3*c2)/12.0;fo[4]=-w*c-fo[1];fo[5]=w*c1*c*c*(4.0-3.0*c1)/12.0;}else if(im==2){fo[1]=-w*(1+2*c1)*cc2;fo[2]=-w*c*c2;fo[4]=-w*(1+2*cc1)*c2;fo[5]=w*cc*c2;}else if(im==3){fo[1]=6*w*c1*cc1/sl;fo[2]=w*cc1*(2-3*cc1);fo[4]=-fo[2];fo[5]=w*c1*(2-3*c1);}else if(im==4){fo[1]=0.25*w*c*(2.-3.*c2+1.6*c2*c1); fo[2]=-w*c*c*(2.-3.*c1+1.2*c2)/6.0; fo[4]=-w*c/2.0-fo[2];fo[5]=0.25*w*c*c*c1*(1.-0.8*c1);}else if(im==5){fo[0]=-w*cc1;fo[3]=-w*c1;}else if(im==6){fo[0]=-w*c*(1.-0.5*c1);fo[3]=-0.5*w*c*c1;}}/*形成总结点荷载向量*/void fpj(){ int i,j,k,ii,jj;for(k=0;k<6;k++)foj[k]=0.0;for(i=0;i<6;i++)for(j=0;j<6;j++) foj[i]=foj[i]+t[j][i]*fo[j]; for(ii=0;ii<6;ii++) pj[is[ii]]=pj[is[ii]]-foj[ii]; }/*形成单元刚度矩阵*/void stiff(){ int i,j;float s1,s2;for(i=0;i<6;i++)for(j=0;j<6;j++)sm[i][j]=0.0; sm[0][0]=eal;sm[3][0]=-eal; sm[3][3]=eal;s1=12.*eil/(sl*sl);s2=6.*eil/sl;sm[1][1]=s1; sm[4][1]=-s1;sm[4][4]=s1; sm[2][1]=s2;sm[5][1]=s2; sm[4][2]=-s2;sm[5][4]=-s2; sm[2][2]=4*eil;sm[5][5]=4*eil; sm[5][2]=2*eil;for(i=0;i<6;i++)for(j=0;j<i;j++)sm[j][i]=sm[i][j];}/*由单元位移分量码L形成总刚位移分量码IS（L）*/ void fis(int nel){ int i,j;for(i=0;i<2;i++)for(j=0;j<3;j++) is[3*i+j]=3*(jel[nel][i]-1)+j; }/*形成结构原始总刚度矩阵*/void addsm(){ int i,j;for(i=0;i<6;i++)for(j=0;j<6;j++){ int kc;kc=is[j]-is[i];if(kc>=0) tsm[is[i]][kc]=tsm[is[i]][kc]+sm[i][j];}}/*约束处理*/void restr(){ int i;for(i=0;i<nr;i++){ int ni;ni=floor(res[i][0]+0.1); tsm[ni-1][0]=1.0e25;pj[ni-1]=tsm[ni-1][0]*res[i][1]; }}/*解线性方程组*/void soleq(){ float c1;int k,ni,im,i,m,j,nm,jm;for(k=1;k<nn;k++){ if(fabs(tsm[k-1][0])<=0.000001){ fprintf(outfile,"*****singularity in row %4d*****\n",k+1); goto fin;}ni=k+mbw-1; im=ni;if(ni>nn) im=nn;for(i=k+1;i<im+1;i++){ m=i-k+1;c1=tsm[k-1][m-1]/tsm[k-1][0];for(j=1;j<mbw-m+2;j++)tsm[i-1][j-1]=tsm[i-1][j-1]-c1*tsm[k-1][j+i-k-1];pj[i-1]=pj[i-1]-c1*pj[k-1]; }}if(fabs(tsm[nn-1][0])<=0.000001){ fprintf(outfile,"******singularity in row %4d******\n",nn); goto fin; }pj[nn-1]=pj[nn-1]/tsm[nn-1][0];for(i=nn-1;i>0;i--){ nm=nn-i+1;jm=nm;if(nm>mbw)jm=mbw;for(j=2;j<jm+1;j++)pj[i-1]=pj[i-1]-tsm[i-1][j-1]*pj[j+i-2];pj[i-1]=pj[i-1]/tsm[i-1][0]; }fin:return;}/*输出位移向量*/void outdis(){ int i;fprintf(outfile,"output results\n");fprintf(outfile,"******\n");fprintf(outfile,"nodal displacements\n");fprintf(outfile," node u v ceta\n");for(i=0;i<nj;i++) fprintf(outfile,"%10d%15.6f%15.6f%15.6f\n",i+1,pj[3*i],pj[3*i+1],pj[3*i+2]); }/*求单元杆端力、支座反力或结合点力*/void force (){ float dj[6],f[6],fj[6],dd[6];int nel,np,i,j,ip;fprintf(outfile,"Mem.N1 Q1 M1 N2 Q2 M2 \n");for(nel=0;nel<ne;nel++){sncs(nel);trmat();stiff();fis(nel);for(i=0;i<6;i++) dj[i]=pj[is[i]];for(i=0;i<6;i++) {dd[i]=0.0;f[i]=0.0;}for(i=0;i<6;i++)for(j=0;j<6;j++) dd[i]=dd[i]+t[i][j]*dj[j];for(i=0;i<6;i++)for(j=0;j<6;j++) f[i]=f[i]+sm[i][j]*dd[j];if(npf==0) goto a;for(np=0;np<npf;np++){ip=floor(pf[np][2]+0.1)-1;if(ip==nel){ fix(np);for(j=0;j<6;j++) f[j]=f[j]+fo[j];}}a:fprintf(outfile,"%5d%11.2f%11.2f%11.2f%11.2f%11.2f%11.2f\n",nel+1,f[0],f[1],f[2],f[3],f[ 4],f[5]);for(i=0;i<6;i++)fj[i]=0.0;for(i=0;i<6;i++)fj[i]=0.0;for(i=0;i<6;i++)for(j=0;j<6;j++)fj[i]=fj[i]+t[j][i]*f[j];for(i=0;i<6;i++) r[is[i]]=r[is[i]]+fj[i]; }fprintf(outfile,"Nodal Reaction\n");fprintf(outfile,"Node RX RY RM\n");for(i=0;i<nj;i++)fprintf(outfile,"%10d%15.4f%15.4f%15.4f\n",i+1,r[3*i],r[3*i+1],r[3*i+2]);}。
有限元方法编程
有限元方法编程【实用版1篇】目录（篇1）1.有限元方法概述2.有限元方法的编程步骤3.有限元方法的应用实例4.总结正文（篇1）一、有限元方法概述有限元方法是一种数值分析方法，广泛应用于固体力学、流体力学、热传导等领域。
它的基本思想是将待求解的连续体划分为有限个小的、简单的子区域，即单元，然后用有限个简单的方程组来代替原来的连续方程，通过求解这些方程组得到近似解。
这种方法既能降低问题的复杂度，又能保证解的精度，因此在工程界有着广泛的应用。
二、有限元方法的编程步骤1.几何建模：根据实际问题，创建待求解的几何模型。
这通常包括划分单元、计算节点坐标等步骤。
2.选择单元类型：根据问题类型和求解需求，选择合适的单元类型，如有限元、无限元、矩形单元、六面体单元等。
3.编写有限元方程：根据单元类型和几何模型，编写有限元方程。
这包括计算单元的刚度矩阵、质量矩阵、载荷矩阵等。
4.组装总方程：将所有单元的有限元方程组装成总方程，通常是一个大型的线性或非线性方程组。
5.求解方程组：使用数值方法（如有限元法、直接解法、迭代法等）求解总方程组，得到近似解。
6.后处理：对求解结果进行分析和处理，如计算应力、应变、位移等。
三、有限元方法的应用实例以一个简单的二维拉伸问题为例，假设有一个长方形板，在左右两端施加均匀拉力，求解板上各个点的应力和应变。
1.几何建模：将长方形板划分为矩形单元，计算节点坐标。
2.选择单元类型：此处采用矩形单元。
3.编写有限元方程：计算单元的刚度矩阵、质量矩阵、载荷矩阵，组装总方程。
4.求解方程组：使用有限元法求解总方程组，得到应力和应变。
5.后处理：分析应力和应变分布，验证解的正确性。
四、总结有限元方法作为一种数值分析方法，通过将连续体划分为有限个小的、简单的子区域，然后用有限个方程组来代替原来的连续方程，降低了问题的复杂度，同时保证了解的精度。
在实际应用中，有限元方法需要经历几何建模、单元选择、编写有限元方程、组装总方程、求解方程组和后处理等步骤。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
有限元编程的c++实现算例 1. #include<> 2. #include<> 3. 4. 5. #define ne 3 #define nj 4 #define nz 6 #define npj 0 #define npf 1 #define nj3 12 #define dd 6 #define e0 #define a0 #define i0 #define pi 16. 17. 18. int jm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}}; /*gghjghg*/ 19. double gc[ne+1]={,,,}; 20. double gj[ne+1]={,,,}; 21. double mj[ne+1]={,a0,a0,a0}; 22. double gx[ne+1]={,i0,i0,i0}; 23. int zc[nz+1]={0,1,2,3,10,11,12}; 24. double pj[npj+1][3]={{,,}}; 25. double pf[npf+1][5]={{0,0,0,0,0},{0,-20,,,}}; 26. double kz[nj3+1][dd+1],p[nj3+1]; 27. double pe[7],f[7],f0[7],t[7][7]; 28. double ke[7][7],kd[7][7]; 29. 30. 31. 36. void jdugd(int); 38. void zb(int); 39. void gdnl(int); 40. void dugd(int); 41. 42. 43. void main() 45. { 46. int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0; 47. double cl,wy[7]; 48. int im,in,jn; 49. 50. 54. if(npj>0) 55. { 56. for(i=1;i<=npj;i++) 57. { j=pj[i][2]; 59. p[j]=pj[i][1]; 60. } 61. } 62. if(npf>0) 63. { 64. for(i=1;i<=npf;i++) 65. { hz=i; 67. gdnl(hz); 68. e=(int)pf[hz][3]; 69. zb(e); for(j=1;j<=6;j++) { 72. pe[j]=; 73. for(k=1;k<=6;k++) { 75. pe[j]=pe[j]-t[k][j]*f0[k]; 76. } 77. } 78. al=jm[e][1]; 79. bl=jm[e][2]; 80. p[3*al-2]=p[3*al-2]+pe[1]; p[3*al-1]=p[3*al-1]+pe[2]; 82. p[3*al]=p[3*al]+pe[3]; 83. p[3*bl-2]=p[3*bl-2]+pe[4]; 84. p[3*bl-1]=p[3*bl-1]+pe[5]; 85. p[3*bl]=p[3*bl]+pe[6]; 86. } 87. } 88. 89. 90. for(e=1;e<=ne;e++) { 94. dugd(e); for(i=1;i<=2;i++) { 97. for(ii=1;ii<=3;ii++) 98. { 99. h=3*(i-1)+ii; dh=3*(jm[e][i]-1)+ii; for(j=1;j<=2;j++) 102. { 103. for(jj=1;jj<=3;jj++) { 105. l=3*(j-1)+jj; zl=3*(jm[e][j]-1)+jj; dl=zl-dh+1; if(dl>0) 109. kz[dh][dl]=kz[dh][dl]+ke[h][l]; } 111. } 112. } 113. } 114. } 115. 116. for(i=1;i<=nz;i++) { 119. z=zc[i]; kz[z][l]=; for(j=2;j<=dd;j++) 122. { 123. kz[z][j]=; } 125. if((z!=1)) 126. { 127. if(z>dd) 128. j0=dd; 129. else if(z<=dd) 130. j0=z; for(j=2;j<=j0;j++) 132. kz[z-j+1][j]=; 133. } 134. p[z]=; } 136. 137. 138. 139. 140. for(k=1;k<=nj3-1;k++) 141. { 142. if(nj3>k+dd-1) im=k+dd-1; 144. else if(nj3<=k+dd-1) 145. im=nj3; 146. in=k+1; 147. for(i=in;i<=im;i++) 148. { 149. l=i-k+1; 150. cl=kz[k][l]/kz[k][1]; jn=dd-l+1; 152. for(j=1;j<=jn;j++) 153. { 154. m=j+i-k; 155. kz[i][j]=kz[i][j]-cl*kz[k][m]; 156. } 157. p[i]=p[i]-cl*p[k]; } 159. } 160. 161. 162. 163. 164. p[nj3]=p[nj3]/kz[nj3][1]; for(i=nj3-1;i>=1;i--) 166. { 167. if(dd>nj3-i+1) 168. j0=nj3-i+1; 169. else j0=dd; for(j=2;j<=j0;j++) 171. { 172. h=j+i-1; 173. p[i]=p[i]-kz[i][j]*p[h]; 174. } 175. p[i]=p[i]/kz[i][1]; } 177. printf("\n"); 178. printf("_____________________________________________________________\n"); 179. printf("NJ U V CETA \n"); for(i=1;i<=nj;i++) 181. { 182. printf(" %-5d % % %\n",i,p[3*i-2],p[3*i-1],p[3*i]); 183. } 184. printf("_____________________________________________________________\n"); 185. printf("E N Q M \n"); 187. for(e=1;e<=ne;e++) { 190. jdugd(e); zb(e); for(i=1;i<=2;i++) 193. { 194. for(ii=1;ii<=3;ii++) 195. { 196. h=3*(i-1)+ii; 197. dh=3*(jm[e][i]-1)+ii; wy[h]=p[dh]; 199. } 200. } 201. for(i=1;i<=6;i++) 202. { 203. f[i]=; 204. for(j=1;j<=6;j++) 205. { 206. for(k=1;k<=6;k++) { 208. f[i]=f[i]+kd[i][j]*t[j][k]*wy[k]; 209. } 210. } 211. } 212. if(npf>0) 213. { 214. for(i=1;i<=npf;i++) if(pf[i][3]==e) { 217. hz=i; 218. gdnl(hz); for(j=1;j<=6;j++) { 221. f[j]=f[j]+f0[j]; 222. } 223. } 224. } 225. printf("%-3d(A) % % %\n",e,f[1],f[2],f[3]); printf(" (B) % % %\n",f[4],f[5],f[6]); } 228. return; 229. } 230. 232. 236. void gdnl(int hz) 237. { 238. int ind,e; 239. double g,c,l0,d; 240. 241. 242. g=pf[hz][1]; c=pf[hz][2]; e=(int)pf[hz][3]; ind=(int)pf[hz][4]; l0=gc[e]; d=l0-c; 248. if(ind==1) 249. { 250. f0[1]=; 251. f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)))/2; f0[3]=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0))/12; 253. f0[4]=; 254. f0[5]=-g*c-f0[2];