有限元教材-第十章有限元程序设计

更新时间：2023-05-05 01:21:02 阅读量：实用文档文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

有限元教材 pdf推荐度：
相关推荐

第十章有限元程序设计

有限元方法作为一门系统的技术，仅学会了它的基本理论是远远不够的，只有形成完整的计算程序，问题才最终得到了解决。完成这样的有限元程序设计是一项工作量很大的工程。本章就是要结合简单的有限元教学程序FEMED，简要介绍有限元程序设计技术。FEMED是专为有限元程序设计教学编制的程序，它不包含复杂的前后处理功能，可进行平面问题及平面桁架的线弹性静力分析，在程序结构上与大型程序类似，具有计算单元的任意扩充功能，在方程的组集和求解上也采用了较为流行的变带宽存储方式。

有限元程序大致可分为两类，第一类是专用程序，主要用于研究或教学，一般这类程序规模较小，前后处理功能较弱。用于研究的程序能够解一些特殊的问题，满足研究工作的需要。而教学程序则是为了学生了解有限元的主要结构和设计方法设计的，程序比较简单，FEMED就属于这类程序。第二类是大型通用程序，是大型结构分析的得力工具，目前国际上流行的大约有2000多种。常用的有NASTRAN、MARC、ANSYS、ADINA和ABAQUS等。这类程序一般前后处理功能比较强，有友好的界面，能进行大型计算，但往往无法完成具有特殊要求的计算。通过本章的学习，使读者初步掌握有限元编程的基本方法，具有开发特殊功能的专用程序或为通用程序开发具有特殊功能的计算模块的能力。

§10.1有限元程序的基本结构

有限元程序一般包括三项基本内容：前处理、结构分析和后处理。早期有限元分析软件的研究重点在于推导新的高效率求解方法和高精度的单元，随着数值分析方法的逐步完善，尤其是计算机内存和运算速度的飞速发展，整个计算系统用于求解运算的时间越来越少，加之求解问题的日益大型化和复杂化，使得数据准备和运算结果的表现问题日益突出。因此目前几乎所有的商业化有限元程序系统都有功能很强的前后处理模块，这直接关系到分析软件的可推广性。它是商用有限元软件不可或缺的部分，但它不是有限元的中心部分，在本书中不作详细介绍。

10.1.1 前处理

在结构分析之前，必须完成的工作就是前处理的任务，它主要包括以下内容：结构造型、选取单元类型，、网格剖分，从而形成节点数、节点编码、节点坐标、单元数、单元节点编码等。另外根据所解题目的不同，还需读入不同的材料参数、边界约束条件及载荷工况等。有限元前处理程序的功能就是为用户提供一种工具，使其能尽可能方便地完成上述工作。而这一问题解决的好坏直接关系到程序使用者需付出的工作量，而网格剖分质量还直接影响计算精度。很多程序都建立了对用户非常友好的图形用户界面，使用户能以可视图形方式直观快速地进行结构造型和网格自动划分，生成有限元分析所需数据。

在上述工作中工作量最大的就是结构造型和网格剖分部分，目前越来越倾向于由专用软件完成上述工作，所以当今所有的商业化有限元软件开发商都开发了和著名的CAD软件的接口，从而极大地提高了设计水平和效率。今天，工程师可以在集成的CAD和FEA软件环境中快捷地解决一个在以前无法应付的复杂工程分析问题。

在FEMED程序中，只包含了很简单的节点、网格及给定位移的自动生成。所有前处理功能是

由INPUT及相关子程序完成的。

10.1.2 结构分析

结构分析部分是有限元程序的主体部分，主要需计算插值函数矩阵[N]，插值函数的导数矩阵[B]，进而进行数值积分得到单元刚度矩阵，组集总刚度矩阵及总载荷向量，解线性代数方程组，得到节点位移。对于动力问题则需计算质量矩阵，求解特征值问题。而对于动力响应问题及非线性问题，则需进行多增量步的计算。在得到节点位移的基础上还需进行应力的计算。

在本章中重点介绍的就是结构分析中的静力分析部分，在FEMED程序中主要包含的也就是这部分。通过这部分的学习，希望读者能掌握有限元程序的基本编程技巧。

10.1.3 后处理

在大型分析软件中，程序的后处理功能也是非常重要的，在有限元结构分析完成后，通过强有力的后处理图形功能，可以给出各应力分量、等效应力等的分布云图，结构的变形图，振型图和实

图10.1

时动画等，使得使用者对各物理量在整个结构中的变化情况有一个全面的认识，也可检查网格剖分和所加载荷及约束是否有误。否则面对输出成千上万个计算结果数据，往往使人如坠雾中，不得要领，需花费大量的时间进行数据的整理和解释。

由于是简单的教学程序，在FEMED程序中不包含后处理功能，仅仅是将计算结果输出。

图10.1给出了一般求解静力学问题的程序框图，FEMED程序的结构基本如图所示，但没有包含带宽优化部分。

§10.2 数组的半动态分配

在有限元程序中需开辟许多不同功能的数组，而这些数组的大小则与所解题目的大小及类型有关，如节点坐标数组、节点位移数组等与节点数相关，而单元节点编号数组等则与单元数相关。此外，上述数组的大小还和自由度维数等相关。对于一个有限元程序而言，应具有求解各类问题及各

种规模问题的机动性，因此求解不同的问题，应定义不同大小的数组。FORTRAN 语言规定数组的大小必须预先给定，在子程序中可以定义大小可变化的可调数组，但其内存分配必须由该子程序的上一级模块确定。

最简单的办法就是每个数组都定义得足够大，使得要求解的问题都能被满足，但这显然会带来内存的巨大浪费，因此，一般采用数组半动态分配的方法来解决上述问题，FEMED 也采用了这一方法。

在主控程序PCONTR 中开辟了两个大数组A(*)和M(*)，分别存储实型量和整型量，将各数组中的数顺序存放在大数组中，在读入控制信息NJ （节点数）、NE （单元数）、NUMMAT （材料数）、NDF （节点自由度数）、NDM （空间维数）、NEN （单元节点数）等之后，就可进行存储单元的分配。在图10.2中以大数组A 的存储为例，先存放节点坐标，占用NJ ×NDM 个存贮单元，再存放材料参数，占用NUMMA T ×20个存贮单元……

图10.2

在上述过程中，最重要的是计算出各数组存放的起始位置又称作地址，如NXC 为XC 数组在A 中的地址，ND 为D 数组在A 中的地址……，这项工作由POINT 子程序完成。需要分别计算出大数组A 和M 中所存数组的所有地址，最后计算出A 和M 数组需要的总长度，与预设的数组长度作比较，若没有超界计算继续进行，若超界则计算停止，显示错误信息。

需要指出的是，在上述数组的半动态分配时，将一系列存储容量不定的一维、二维数组均用地址对应关系分配到事先容量固定的大数组A 和M 中去，对于大数组A 和M 来说存储量是不变的，但各数组在大数组中的存储分配却是动态的。我们把这种关系称为数组的半动态分配，这一过程是通过哑实结合完成的。这些数组称为可调数组，即其数组长度是可调的。在子程序的数组说明中，一维数组写成下列形式是等价的，对数组长度没有限制：

F(*)，F(1)，F(NJ)，

对于二维数组，则只可一个方向自由变动。由于FORTRAN 规定二维数组先变行后变列，因此行的长度必须是确定的。，在子程序的数组说明中，下列形式是等价的：

XC(2,*)，XC(2,1)，XC(2,NJ)，XC(NDM,*)，(NDM=2)

注意，只有对上一模块中已定义过的数组才可作如上的数组说明，否则，则必须准确确定数组的长度。

§10.3 构造单元刚度矩阵及总刚度矩阵的组集

构造单元刚度矩阵是由单元模块完成的，不同的单元调用不同的模块，在FEMED 中提供了平面桁架单元（TRUSS ）和平面单元（PLANE ）两个单元计算模块。读者可根据需要随意开发单元类型，通过ELEMLIB 子程序接入程序中。针对不同的单元，在不同的单元模块中构造不同的[B]、

[D]矩阵，单元刚度矩阵可表示为：

V dV =?e e T K B DB （10.1）

对于常应变单元，上述被积函数为常数，直接用被积函数值乘积分域的大小即可。而对于等参元等非常应变单元，则需对上式做GAUSS 积分，最终得到一个二维的单元刚度矩阵。

将单元刚度矩阵正确组集为总刚度矩阵，这一工作是由ADDSTF 子程序完成的。下面将举例说明总刚度矩阵的组集方法。

例1：平面桁架结构如图示，共有5个节点，6个单元。每个单元形成一个4×4的单元刚度矩阵，总刚度阵为10×10。

图10.3 以1号单元为例，其左节点为2，右节点为4，其4×4的单元刚度矩阵的各个参数应迭加在(10.2)式总刚度矩阵中的图示位置上。需要指出的是，总刚度矩阵的一个位置上往往有多个单元刚度矩阵做贡献，如与3号节点对应的(5,5)，(6,5)，(6,6)位置上就有2、3、4和6号单元都作了贡献。

全部单元组集完毕后，总刚度矩阵如式（10.2），形成一对称带状阵，即仅对角元附近的元素不为零，而远离对角线的元素皆为零。解题越大带状性质越明显。

????????????????????????????????????????????????????????????

???=00000000

444342410033323100

222111][称对K （10.2）

若储存的是总刚度矩阵的上三角元素，则带宽为每行从对角元，到最后一个非零元素的个数；若储存的是总刚度矩阵的下三角元素，则从第一个非零元素到对角元元素的个数称为这行的带宽。每行的带宽各不相同，且与节点编号的方式有关，在大型商业软件中，一般都包含带宽优化模块，就是为了减少总体带宽。

§10.4 总刚度矩阵的变带宽存储

总刚度阵占用的空间是有限元占用空间的绝大部分，一个软件的解题规模往往是受此限制。所以，一直以来人们在这方面想了很多的办法。很显然，将总刚度阵作为一个二维数组，完全储存起来是最直接和最自然的方法，但这也是占用空间最大、最不合理的方法，一般很少采用。

考虑到总刚度阵具有带状稀疏对称的特征，使用半带宽存储无疑是一个最自然的改进，它利用一个二维数组，只存储每行从对角元到这一行的最后一个非零元，如图10.4所示

????

?????

???

????????

???

?????

??????

???

??????

???????

????????????????????????0 000000000000000000000000000000000000称对

图10.4

由10.4的示意图知，原10×10的矩阵现在可用10×4的矩阵存储，存储空间减少了很多，但若其中有个别行的带宽特别大，则这种方法依然有存贮空间的大量浪费。为解决上述问题，一般多采用变带宽存储又称一维压缩存储的方法，这样刚度阵中大量的零既不需要存储也不参与计算，大大减少了存储空间和计算量，是目前使用较为广泛的方法。

在有限元的发展历史中波前法也是一种不可不提的计算方法，当计算机的内存非常有限时，即使采用变带宽存储也无法求解较大规模的问题，波前法曾经起到了不可忽视的作用。由于波前法无需组集总刚，所以对内存要求非常低，曾经被广泛地使用，但它也有非常明显的缺点，内外存交换非常频繁，导致计算速度受到很大的限制。随着计算机技术突飞猛进地发展，计算机的内存有了本质的提高，所以波前法已经很少被采用了。本书重点介绍变带宽存储的方法。

变带宽存储就是考虑到总刚度矩阵的对称性及带状稀疏的特点，用一维数组将总刚度矩阵下三角（或上三角）中的非零元素储存起来，如图10.5中虚线以上部分，注意一行中夹在中间的零必须储存。

图10.5

设A(i,j)为具有带状对称性质的二维数组，若将其存储在一维数组B(k)中，则必须开辟一个存放对角元地址的数组JP(i)，指明每一个对角元在一维数组中的位置。如上图中矩阵的JP(i)应为：1，3，6，9，11，15，18，21，24，27。若第i 行，第j 列的数据存在一维数组的第k 个位置上，则k 与i 、j 之间的对应关系为

k=JP(i)-i+j

(10.3)

第i 行的第一个非零元素在m i 列，则

m i =i-JP(i)+JP(i-1)+1 (i ≠1)

(10.4)

???????????????

???????????????????

???????????????????????

0000000000000000000000称

对

当i=1时，很明显m 1=1。有了（10.3、4）两式，我们完全能够对一维数组很方便地进行方程求解，对正定对称方程A(i,j)的求解方法都可以使用，只是对A(i,j)的操作都需先计算k ，然后对B(k)做同样的操作即可。另外，由于m i 为第i 行的第一个非零元素，所以对第i 行的操作只需从第m i 列开始即可。

对于上述10 10的矩阵，我们需要存储27个元素，节省有限，而对于大规模的问题则节省的内存和计算工作量都将是非常显著的。

另外需要指出的是，带宽的大小与编码方式有关，我们将中间有单元相连的节点称为相关节点，尽量减小相关节点间的编号差，则可减小带宽。减小带宽对于减小内存占用和减小计算工作量都是很有意义的。在图10.6的二维结构中，（a ）、（b ）两种编码方式的半带宽就相差很多，（a ）中的最大半带宽为8，而(b)中的最大半带宽则为16。这会使内存占用和计算工作量都相差很多。

例2：试分别给出图10.6中两种编码的带宽、对角元地址及总刚长度。

解：j 节点对应的是2j-1行和2j 行，两节点若包含在同一单元中则称相关节点，设与j 节点相关的节点中，编号最小的为i ，则

2j-1行的带宽d=2(j-i)+1。2j 行的带宽d=2(j-i)+2 对于（

对于（

（a ）编码形式的总纲长度为140，（b ）编码形式的总纲长度为220。

上述带宽及对角元地址的计算在PROFIL 子程序中完成，所不同的是PROFIL 子程序中还自动将约束自由度从总刚度矩阵中剔除，而本例中没考虑约束情况。

§10.5 给定位移约束的处理

经过总刚度矩阵的组集，并将给定载荷加在载荷向量的相应位置，我们即可得到线性代数方程

(a) (b)

图10.6

组

Ku F

(10.5)

每行对应一个自由度，需要指出的是，与每个自由度对应的有两个量——节点位移u i 和节点载荷f i ，在这两个量中有且仅有一个是已知的。在例1中前四个自由度已知节点位移而节点力未知，后面各自由度则已知节点力而节点位移未知。在求解线性代数方程组时，要求右端项全部已知，而未知量应在左端的u 中，因此需要根据边界约束条件对方程组进行改造。

若第i 个自由度上u i =?i 为已知，而f i 未知，则有三种处理方法：

105.1 加大数法

加大数法是程序处理最简单的方法，但也是一种近似的方法。作法是令对角元k ii =M i ，f i =?i M i ，其中M i 远远大于第i 行中其他各量，则其他变量的影响就可忽略，因此有δi ≈?i 。需要指出的是M i 的选择是技术性比较强的，若M i 选择得不够大则其他变量的影响不能较好地消除，则误差较大，若M i 选择得太大则影响方程的性态，引起较大误差甚至奇异，因此这是一种简单但不能尽如人意的方法。

10.5.2 划行划列法

划行划列法是约束处理的一种精确的方法，由于节点力f i 未知，第i 个方程引入了一个新的未知数，所以对求解a 是无意义的，因此放弃这一方程，换为u i =?i 。作法是将第i 行除对角元外的各元素都充零，对角元充1，令f i = ?i 。即

i i ij ii f i j n j k k ?=≠===)

,1(0

完成上述操作就可实现给定约束的处理，处理之后为

1211112122222212001

0i n i

n i i n n ni

nn n n k k k k u f k k k k u f u f k k k k u f ????????????

????????????=??????????

?????????????

?????

(10.6)

很明显上述操作破坏了方程的对称性，是不能接受的，因此需将第i 列也都充零（除对角元），以恢复其对称性。但这一操作不能影响其他各方程，因此需将?i 对其它各行的影响移至等号右边，即做如下操作

121111121222222212

0000100n i i n i i i i i

n n nn n n n ni i k k k u f f k k k k u f f k u f k k k u f f k ????-????????????????

-????????

??????

??==????????????????

????????????????-?

???????

(10.7)

完成上述操作则完成了约束处理，既保持了方程的对称性，又精确引入了位移边界条件，若还有其他约束则再做同样处理即可。

10.5.3 消自由度法

划行划列的方法无疑是一个切实可行的方法，但若受约束自由度很多则方程中多储存了许多的零，而且在解方程中也增加了许多计算量，因此就有了对上述算法的改进算法。FEMED 就采用了

这一算法。

对于需要划掉的行列一开始就不再组集入总刚度阵，若有l 个自由度被约束，则方程组就从n 元线性代数方程组降为n-l 元，存储空间和计算量都会有比较大的降低。具体作法是，开辟一个数组ID(2,*)，记录约束情况，即活化自由度与方程组的对应关系。对于例1，形成的ID 数组应为：

-1，-2；-3，-4；1，2；3，4；5，6；7，8；9，10。

负值代表该自由度被约束，如-3代表2节点的x 方向是第3个位移约束，而6则代表第5个节点的y 方向的自由度是与方程组的第6个未知数对应。在组集总刚度矩阵时，应先查ID 数组，例如原来应组集在k 85中的数现在则应该组集在k 41中，原来应组集在k 41中的数则无需组集，原来应组集在k 72中的数k ，对应的是第3自由度和第2个约束，现在则在总刚度矩阵中不体现，而应将f 3修改为f 3-k ?2。

在FEMED 中，ID 数组是在BOUN 和PROFIL 中自动形成的。

§10.6 单元应力及节点力的计算

在所有节点位移求出后，要进一步求出应力及反力等，这时需再次对单元循环，对各单元的积分点计算应力

e =σDBa

其中，a e 为单元的节点位移向量，B 、D 矩阵与前面各章相同分别为应变矩阵和弹性矩阵，在B 矩阵中代入不同的局部坐标（ξ, η）可以得到不同点的应力。但一般认为，积分点的应力精度最好，而节点应力是最差的，所以一般有限元程序得到的都是积分点应力而非节点应力，当需要节点应力时，往往再由积分点应力通过插值格式得到节点应力。

计算单元节点力可由下式得到

e e e T T e V V

dV dV ==??σF B B DB a 将所有各单元上计算出的节点力迭加起来，即可得到各节点的节点力。对于给定位移的自由度，如此求得的即是节点反力，而对于自由节点则会求出一个很小的节点力，这是计算误差使然。在线性分析中，上述计算主要是为了得到节点反力，而在非线性分析中，上述节点力则作为残差参与迭代计算。

§10.7 FEMED 程序

FEMED 程序是一个可用来求解线弹性静力学问题的有限元教学程序，它不包含复杂的前后处理，目的是让读者通过程序的学习初步掌握有限元主体程序的初步编程技巧，能够读懂和编制有限元程序，FEMED 程序采用动态数组存储技术。在线性代数方程组的存储上，采用了变带宽存储技术。在约束节点的处理上，采取了预先消除自由度的方法。所有这些都使其内存能够得到最有效的利用。本程序由于篇幅所限只包含了平面单元和平面桁架单元两种单元，但它采用了单元库的设计方式，使得使用者可以方便地添加自己所需要的单元，丰富自己的单元库。为了方便读者阅读，在10.9节给出了本程序的子程序索引、输入索引及主要变量索引。在程序之后给出了3个简单的算例。

PROGRAM FEMED !MAIN 1

C.....FINITE ELEMENT PROGRAM FOR EDUCATIONAL PURPOSE ONLY. !MAIN 2

C.....AVAILABLE FOR SOLVING MULTI-TYPES ELEMENTS STATIC PROBLEMS !MAIN 3

COMMON/PSIZE/MAXM,MAXA !MAIN 4

OPEN(5,FILE='EXP.DAT') !MAIN 5 OPEN(6,FILE='EXP.OUT') !MAIN 6 MAXM=6000 !MAIN 7 MAXA=90000 !MAIN 8 CALL PCONTR !MAIN 9 STOP !MAIN 10 END !MAIN 11 C............................................

SUBROUTINE PCONTR !PCON 1 C.....MAIN CONTROLING SUBROUTINE !PCON 2

COMMON/M/M(6000) !PCON 3 COMMON/A/A(90000) !PCON 4 COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PCON 5 COMMON/PDATA/MIX,MID,MLD,MJP,NXC,ND,NXL,NUL,NS,NP,NF,NU,NSK !PCON 6 COMMON/PSIZE/MAXM,MAXA !PCON 7 COMMON/MDATA/NUMMAT !PCON 8 C.....READ IN CONTROL DATA !PCON 9

READ(5,*) NJ,NE,NUMMAT,NDF,NDM,NEN !PCON 10 WRITE(6,2000) NJ,NE,NUMMAT,NDF,NDM,NEN !PCON 11 C.....CALCULATE THE ADDRESS OF ARRAY !PCON 12

CALL POINT(NNEQ) !PCON 13 C.....READ MOST OF THE PROBLEM DATA !PCON 14

CALL INPUT(A(NXC),M(MIX),M(MID),A(ND),A(NF)) !PCON 15 C.....DEFINE THE PROFIL FOR GLOBAL STIFFNESS MATRIX STORAGE !PCON 16

CALL PROFIL(M(MJP),M(MID),M(MIX),NK) !PCON 17 NEND=NSK+NK !PCON 18 IF(NEND.GT.MAXA) STOP '{A} ERR' !PCON 19 CALL PZERO(A(NSK),NK) !PCON 20 C.....FORM LOAD ARRAY !PCON 21

CALL PLOAD (M(MID),A(NF),A(NU),NNEQ) !PCON 22 C.....FORM ELEMENT STIFFNESS !PCON 23

CALL PFORM(A(NF),M(MLD),A(NUL),A(NXL),A(NS),A(NP),A(NU),A(NXC), !PCON 24

1 M(MIX),A(ND),M(MID),M(MJP),A(NSK),2) !PCON 25

C.....LDLT TRIANGULAR DECOMPOSITION !PCON 26

CALL LDLT(A(NSK),M(MJP),NEQ) !PCON 27 C.....FORWARD ELIMINATION AND BACKWARD SUBSTITUTION !PCON 28

CALL FORBACK(A(NSK),A(NU),M(MJP),NEQ) !PCON 29 C.....PRINT DISPLACEMENTS !PCON 30

CALL PRTDIS(M(MID),A(NU),A(NF)) !PCON 31 CALL PZERO(A(NF),NNEQ) !PCON 32 85abd76e7e21af45b307a8f6PUTE ELEMENT STRESSES AND THE REACTION !PCON 33

CALL PFORM(A(NF),M(MLD),A(NUL),A(NXL),A(NS),A(NP),A(NU),A(NXC), !PCON 34

1 M(MIX),A(ND),M(MID),M(MJP),A(NSK),3) !PCON 35

C.....PRINT THE REACTION !PCON 36

CALL PRTREAC(A(NF)) !PCON 37 2000 FORMAT(' NUMBER OF NODAL POINTS =',I5/ !PCON 38

1 ' NUMBER OF ELEMENTS =',I5/ !PCON 39

1 ' NUMBER OF MATERIALS =',I5/ !PCON 40

1 ' NUMBER OF DEGREES/NODE =',I5/ !PCON 41

1 ' NUMBER OF DIMENSION SPACE =',I5/ !PCON 42

1 ' NUMBER OF NODES/ELEMENT =',I5/) !PCON 43

RETURN !PCON 44

END !PCON 45 C............................................

SUBROUTINE POINT(NNEQ) !POIN 1 C.....CALCULATE THE ADDRESS OF ARRAY !POIN 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !POIN 3 COMMON/PSIZE/MAXM,MAXA !POIN 4 COMMON/PDATA/MIX,MID,MLD,MJP,NXC,ND,NXL,NUL,NS,NP,NF,NU,NSK !POIN 5 COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL,NST !POIN 6 COMMON/MDATA/NUMMAT !POIN 7 COMMON/A/A(90000) !POIN 8 NEN1=NEN+1 !POIN 9 NNEQ=NJ*NDF !POIN 10 NST=NDF*NEN !POIN 11 C.....CALCULATE THE ADDRESS OF INTEGRAL ARRAY {M} !POIN 12

MIX=1 !POIN 13 MID=MIX+NEN1*NE !POIN 14 MLD=MID+NNEQ !POIN 15 MJP=MLD+NST !POIN 16 MEND=MJP+NNEQ !POIN 17 IF(MEND.GT.MAXM) STOP '{M} ERR' !POIN 18 C.....CALCULATE THE ADDRESS OF REAL ARRAY {A} !POIN 19

NXC=1 !POIN 20 ND =NXC+NJ*NDM !POIN 21 NXL=ND +20*NUMMAT !POIN 22 NUL=NXL+NEN*NDM !POIN 23 NS =NUL+NST*2 !POIN 24 NP =NS +NST*NST !POIN 25 NF =NP +NST !POIN 26 NU =NF +NNEQ !POIN 27 NSK=NU +NNEQ !POIN 28 CALL PZERO(A(1),NSK) !POIN 29 RETURN !POIN 30 END !POIN 31 C............................................

SUBROUTINE INPUT(XC,IX,ID,D,F) !INPU 1 C.....READ MOST OF THE PROBLEM DATA !INPU 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !INPU 3 COMMON/PDATA/MIX,MID,MLD,MJP,NXC,ND,NXL,NUL,NS,NP,NF,NU,NSK !INPU 4 COMMON/PSIZE/MAXM,MAXA !INPU 5 COMMON/MDATA/NUMMAT !INPU 6 COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL,NST !INPU 7 DIMENSION XC(NDM,*),IX(NEN1,*),ID(NDF,*),D(20,*),F(NDM,*) !INPU 8 CHARACTER COOR*6,XX*18 !INPU 9 DATA XX/'NODAL COORDINATES'/COOR/'COOR'/ !INPU 10 CHARACTER FD*18,FORC*6 !INPU 11 DATA FORC/' FORC '/FD/' NODAL FORCE/DISPL'/ !INPU 12 DIMENSION P(1),S(1) !INPU 13 C.....INPUT NODAL COORDINATES DATA !INPU 14

CALL GENVEC(XC,XX,COOR) !INPU 15 C.....INPUT ELEMENT DATA !INPU 16

CALL ELDAT(IX) !INPU 17 C.....INPUT MATERIAL DATA !INPU 18

DO N=1,NUMMAT !INPU 19

READ(5,*) MA,IEL !INPU 20

WRITE(6,2000) MA,IEL !INPU 21

IF(MA.EQ.0) EXIT !INPU 22

D(20,MA)=IEL !INPU 23

CALL ELEMLIB(D(1,MA),P,XC,IX,S,P,NDF,NDM,1) !INPU 24 END DO !INPU 25 C.....INPUT BOUNDARY CONDITIONS !INPU 26

CALL BOUN(ID) !INPU 27 C.....INPUT FORCES AND DISPLACEMENTS !INPU 28

CALL GENVEC(F,FD,FORC) !INPU 29 2000 FORMAT(/' MATERIAL TYPE ',I4,' ELEMENT TYPE ',I4/) !INPU 30 RETURN !INPU 31 END !INPU 32 C.................................................

SUBROUTINE GENVEC(X,CD,FC) !GENV 1 C.....GENERATE REAL DATA ARRAYS BY LINEAR INTERPOLATION !GENV 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !GENV 3 CHARACTER CD*18,FC*6 !GENV 4 DIMENSION X(NDM,*),XL(6) !GENV 5 DO K=1,NJ+1 !GENV 6

READ(5,*) N,NG,(XL(I),I=1,NDM) !GENV 7

IF(N.EQ.0) EXIT !GENV 8

DO I=1,NDM !GENV 9

X(I,N)=XL(I) !GENV 10 END DO !GENV 11 IF(NG.NE.0) THEN !GENV 12

LI=(N-L)/NG !GENV 13

IF(LI.LT.1) STOP 'GENVEC ERR' !GENV 14

DO I=1,NDM !GENV 15

XL(I)=(X(I,N)-X(I,L))/LI !GENV 16 END DO !GENV 17

L1=L+NG !GENV 18

L2=N-NG !GENV 19 DO L=L1,L2 !GENV 20 DO I=1,NDM !GENV 21

LMLG=L-NG !GENV 22

X(I,L)=X(I,LMLG)+XL(I) !GENV 23 END DO !GENV 24 END DO !GENV 25 ENDIF !GENV 26 L=N !GENV 27 END DO !GENV 28 CALL OUTPUT(X,CD,FC) !GENV 29 RETURN !GENV 30 END !GENV 31 C.................................................

SUBROUTINE OUTPUT(X,CD,FC) !OUTP 1 COMMON /CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !OUTP 2 CHARACTER CD*25,FC*6 !OUTP 3 DIMENSION X(NDM,*) !OUTP 4 WRITE(6,2000) CD,(K,FC,K=1,NDM) !OUTP 5

DO J=1,NJ !OUTP 6

WRITE(6,3000)J,(X(L,J),L=1,NDM) !OUTP 7 END DO !OUTP 8 2000 FORMAT(/5X,A18/2X,5HNODE ,6(I6,A6)) !OUTP 9 3000 FORMAT(I6,6F12.4) !OUTP 10 RETURN !OUTP 11 END !OUTP 12 C...............................................................

SUBROUTINE ELDAT(IX) !ELDA 1 C.....ELEMENT DATA INPUT !ELDA 2

COMMON /CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !ELDA 3 DIMENSION IX(NEN1,*),IXL(9) !ELDA 4 WRITE(6,2000) (K,K=1,NEN) !ELDA 5 DO J=1,NE+1 !ELDA 6

READ(5,*) N,NN,NG,LK,(IXL(K),K=1,NEN) !ELDA 7

IF(N.EQ.0) EXIT !ELDA 8

N1=N+NN-1 !ELDA 9

IX(NEN1,N) =LK !ELDA 10 DO K=1,NEN !ELDA 11

IX(K,N) =IXL(K) !ELDA 12 END DO !ELDA 13 IF(NN.EQ.1) CYCLE !ELDA 14 DO I=N+1,N1 !ELDA 15

IX(NEN1,I) =LK !ELDA 16

I1=I-1 !ELDA 17

DO K=1,NEN !ELDA 18

IX(K,I)=IX(K,I1)+NG !ELDA 19

IF(IX(K,I1).EQ.0) IX(K,I)=0 !ELDA 20

END DO !ELDA 21 END DO !ELDA 22 END DO !ELDA 23 DO I=1,NE !ELDA 24 WRITE(6,3000) I,IX(NEN1,I),(IX(K,I),K=1,NEN) !ELDA 25 END DO !ELDA 26 2000 FORMAT(//5X,'ELEMENTS CONNECTION'//1X,'ELEM MATE', !ELDA 27

1 10(I2,' NODE')/(10X,10(I2,' NODE'))) !ELDA 28 3000 FORMAT(2I5,10I7/(10X,10I7)) !ELDA 29 RETURN !ELDA 30 END !ELDA 31 C........................................................

SUBROUTINE BOUN(ID) !BOUN 1 C.....RESTRAINT CONDITIONS INPUT !BOUN 2

COMMON /CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !BOUN 3 DIMENSION ID(NDF,*),IDL(6) !BOUN 4 CHARACTER*6 BC !BOUN 5 DATA BC/' B.C. '/ !BOUN 6 WRITE(6,2000) (I,BC,I=1,NDF) !BOUN 7 DO I=1,NJ !BOUN 8

DO J=1,NDF !BOUN 9

ID(J,I)=0 !BOUN 10 END DO !BOUN 11 END DO !BOUN 12

DO K=1,NJ !BOUN 13 READ(5,*) N,NN,NG,(IDL(I),I=1,NDF) !BOUN 14 IF(N.EQ.0) EXIT !BOUN 15 N1=N+(NN-1)*NG !BOUN 16 DO J=N,N1,NG !BOUN 17

DO I=1,NDF !BOUN 18

ID(I,J) =IDL(I) !BOUN 19

END DO !BOUN 20

WRITE(6,3000) J,(ID(I,J),I=1,NDF) !BOUN 21 END DO !BOUN 22 END DO !BOUN 23 2000 FORMAT(/5X,'NODAL BOUNDARY CONDITIONS'//2X,'NODE',7(I4,A6)/1X) !BOUN 24 3000 FORMAT(I6,7(I7,3X)) !BOUN 25 RETURN !BOUN 26 END !BOUN 27 C.........................................................

SUBROUTINE PROFIL(JP,ID,IX,NK) !PROF 1 C.....SET UP EQUATION NUMBERS !PROF 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PROF 3 DIMENSION JP(1),ID(NDF,1),IX(NEN1,*) !PROF 4 85abd76e7e21af45b307a8f6PUTE THE NUMBER OF EQUATION !PROF 5

NEQ=0 !PROF 6 NAD=0 !PROF 7 DO N=1,NJ !PROF 8

DO I=1,NDF !PROF 9

IF(ID(I,N).EQ.0) THEN !PROF 10

NEQ=NEQ+1 !PROF 11

ID(I,N)=NEQ !PROF 12

JP(NEQ)=0 !PROF 13

ELSE !PROF 14

NAD=NAD-1 !PROF 15

ID(I,N)=NAD !PROF 16

ENDIF !PROF 17 END DO !PROF 18 END DO !PROF 19 85abd76e7e21af45b307a8f6PUTE ROW BANDWIDTH !PROF 20

DO N=1,NE !PROF 21 LL=NEQ !PROF 22 DO I=1,NEN !PROF 23

II=IX(I,N) !PROF 24

IF(II.EQ.0) CYCLE !PROF 25

DO K=1,NDF !PROF 26

KK=ID(K,II) !PROF 27

IF(KK.LE.0) CYCLE !PROF 28

LL=MIN0(LL,KK) !PROF 29

END DO !PROF 30 END DO !PROF 31 DO J=1,NEN !PROF 32

JJ=IX(J,N) !PROF 33

IF(JJ.EQ.0) CYCLE !PROF 34

DO L=1,NDF !PROF 35

KK=ID(L,JJ) !PROF 36

IF(KK.LE.0) CYCLE !PROF 37

M=KK-LL !PROF 38

JP(KK)=MAX0(JP(KK),M) !PROF 39

END DO !PROF 40 END DO !PROF 41 END DO !PROF 42 85abd76e7e21af45b307a8f6PUTE DIAGONAL POINTERS FOR PROFIL !PROF 43

NK=1 !PROF 44 JP(1)=1 !PROF 45 DO N=2,NEQ !PROF 46 JP(N)=JP(N)+JP(N-1)+1 !PROF 47 END DO !PROF 48 NK=JP(NEQ) !PROF 49 WRITE(6,2000) NK !PROF 50 2000 FORMAT(' THE MAXIMUM STORAGE OF GLOBAL STIFFNESS =',I6/) !PROF 51 RETURN !PROF 52 END !PROF 53 C...........................................................

SUBROUTINE PLOAD(ID,F,U,NN) !PLOA 1 C.....FORM LOAD ARRAY !PLOA 2

DIMENSION F(1),ID(1),U(1) !PLOA 3 DO N=1,NN !PLOA 4 J=ID(N) !PLOA 5 IF(J.GT.0) U(J)=U(J)+F(N) !PLOA 6 END DO !PLOA 7 RETURN !PLOA 8 END !PLOA 9 C............................................................

SUBROUTINE PFORM(F,LD,UL,XL,S,P,U,X,IX,D,ID,JP,AA,ISW) !PFOR 1 C.....CONTROLING SUBROUTINE OF LOOPING OVER ELEMENT !PFOR 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PFOR 3 COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL,NST !PFOR 4 DIMENSION UL(NDF,*),U(NDF,*),F(NDF,*),S(NST,*),P(*),XL(NDM,*), !PFOR 5

1 X(NDM,*),D(20,*),IX(NEN1,*),ID(NDF,*),LD(NDF,*),JP(*), !PFOR 6

1 AA(*) !PFOR 7

IEL=0 !PFOR 8 MCT=0 !PFOR 9 DO N=1,NE !PFOR 10 CALL PZERO(UL,NST*(NST+3)) !PFOR 11 MA=IX(NEN1,N) !PFOR 12 C.....FORM ELEMENT DATA !PFOR 13

DO I=1,NEN !PFOR 14

II=IX(I,N) !PFOR 15

IF(II.EQ.0) THEN !PFOR 16

DO J=1,NDF !PFOR 17

LD(J,I)=0 !PFOR 18

END DO !PFOR 19

ELSE !PFOR 20

IID=II*NDF-NDF !PFOR 21

NEL=I !PFOR 22

DO J1=1,NDM !PFOR 23

XL(J1,I)=X(J1,II) !PFOR 24

END DO !PFOR 25

DO J=1,NDF !PFOR 26

IF(ISW.EQ.2) THEN !PFOR 27

K=ID(J,II) !PFOR 28

IF(K.LT.0) UL(J,I)=F(J,II) !PFOR 29

LD(J,I)=K !PFOR 30

ELSE !PFOR 31

UL(J,I)=U(J,II) !PFOR 32

LD(J,I)=IID+J !PFOR 33

END IF !PFOR 34

END DO !PFOR 35

END IF !PFOR 36 END DO !PFOR 37 IE20=D(20,MA) !PFOR 38 IF(IE20.NE.IEL) MCT=0 !PFOR 39 IEL=IE20 !PFOR 40 C.....CALL SUBROUTINE OF DIFFERENT ELEMENT TYPE !PFOR 41

CALL ELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,ISW) !PFOR 42 IF(ISW.EQ.2) THEN !PFOR 43 C.....STIFFNESS MATRIX ASSEMBLAGE !PFOR 44

CALL ADDSTF(LD,JP,S,AA,UL,U) !PFOR 45 ELSE !PFOR 46 C.....FORM FORCE OF NODE !PFOR 47

CALL BASBLY(F,P,LD,NST) !PFOR 48 ENDIF !PFOR 49 END DO !PFOR 50 RETURN !PFOR 51 END !PFOR 52 C.............................................................

SUBROUTINE ADDSTF(LD,JP,S,A,UL,U) !ADDS 1 C.....STIFFNESS MATRIX ASSEMBLAGE !ADDS 2

COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL,NST !ADDS 3 DIMENSION LD(*),JP(*),S(NST,NST),A(*),UL(*),U(*) !ADDS 4 DO J=1,NST !ADDS 5

K=LD(J) !ADDS 6

IF(K.GT.0) THEN !ADDS 7

L=JP(K)-K !ADDS 8

DO I=1,NST !ADDS 9

M=LD(I) !ADDS 10

IF(M.GT.K.OR.M.LE.0) CYCLE !ADDS 11

M=L+M !ADDS 12

A(M)=A(M)+S(J,I) !ADDS 13

END DO !ADDS 14 ELSE !ADDS 15

UU=UL(J) !ADDS 16

DO I=1,NST !ADDS 17

M=LD(I) !ADDS 18

IF(M.LE.0) CYCLE !ADDS 19

U(M)=U(M)-S(J,I)*UU !ADDS 20

END DO !ADDS 21 ENDIF !ADDS 22 END DO !ADDS 23

RETURN !ADDS 24 END !ADDS 25 C...........................................................

SUBROUTINE PZERO(P,N) !PZER 1 DIMENSION P(N) !PZER 2 DO I=1,N !PZER 3

P(I)=0.0 !PZER 4 END DO !PZER 5 RETURN !PZER 6 END !PZER 7 C...........................................................

SUBROUTINE BASBLY(B,P,LD,NS) !BASB 1 C.....ASSEMBLE THE LOAD VECTOR !BASB 2

DIMENSION B(*),P(*),LD(*) !BASB 3 DO I=1,NS !BASB 4

II=LD(I) !BASB 5

IF(II.EQ.0) CYCLE !BASB 6

B(II)=B(II)-P(I) !BASB 7 END DO !BASB 8 RETURN !BASB 9 END !BASB 10 C..........................................................

SUBROUTINE LDLT(A,JP,N) !LDLT 1 C.....LDLT TRIANGULAR DECOMPOSITION !LDLT 2

DIMENSION A(*),JP(*) !LDLT 3 NK=JP(N) !LDLT 4 A(1)=1.0/A(1) !LDLT 5 IJ=1 !LDLT 6 II=2 !LDLT 7 ID1=JP(1) !LDLT 8 DO I=2,N !LDLT 9

ID2=JP(I) !LDLT 10 MI=I-ID2+ID1+1 !LDLT 11 MJ=1 !LDLT 12 DO J=MI,I-1 !LDLT 13

IJ=IJ+1 !LDLT 14

IF(J.NE.1) MJ=J-JP(J)+JP(J-1)+1 !LDLT 15

AIJ=A(IJ) !LDLT 16

MIJ=MAX0(MI,MJ) !LDLT 17

IK=JP(I)-I+MIJ !LDLT 18

JK=JP(J)-J+MIJ !LDLT 19

DO K=MIJ,J-1 !LDLT 20

AIJ=AIJ-A(IK)*A(JK) !LDLT 21

IK=IK+1 !LDLT 22

JK=JK+1 !LDLT 23

END DO !LDLT 24

A(IJ)=AIJ !LDLT 25 END DO !LDLT 26 II=IJ+1 !LDLT 27 AII=A(II) !LDLT 28 IJ=JP(I-1)+1 !LDLT 29 DO J=MI,I-1 !LDLT 30

JJ=JP(J) !LDLT 31

AIJ=A(IJ) !LDLT 32

A(IJ)=A(IJ)*A(JJ) !LDLT 33

AII=AII-AIJ*A(IJ) !LDLT 34

IJ=IJ+1 !LDLT 35 END DO !LDLT 36 A(II)=1.0/AII !LDLT 37 ID1=ID2 !LDLT 38 END DO !LDLT 39 RETURN !LDLT 40 END !LDLT 41 C.................................................................

SUBROUTINE FORBACK(A,B,JP,N) !FORB 1 C.....FORWARD REDUCTION !FORB 2

DIMENSION A(*),B(*),JP(*) !FORB 3 IJ=2 !FORB 4 DO I=2,N !FORB 5

MI=I-JP(I)+JP(I-1)+1 !FORB 6

DO J=MI,I-1 !FORB 7

B(I)=B(I)-A(IJ)*B(J) !FORB 8

IJ=IJ+1 !FORB 9

END DO !FORB 10 IJ=IJ+1 !FORB 11 END DO !FORB 12 C.....REFORM THE RIGHT HAND SIDE !FORB 13

DO I=1,N !FORB 14 JPI=JP(I) !FORB 15 B(I)=B(I)*A(JPI) !FORB 16 END DO !FORB 17 C.....BACK SUBSTITUTION !FORB 18

IJ=IJ-1 !FORB 19 DO I=N,2,-1 !FORB 20 MI=I-JP(I)+JP(I-1)+1 !FORB 21 DO J=I-1,MI,-1 !FORB 22

IJ=IJ-1 !FORB 23

B(J)=B(J)-A(IJ)*B(I) !FORB 24 END DO !FORB 25 IJ=IJ-1 !FORB 26 END DO !FORB 27 RETURN !FORB 28 END FORB 29 C..........................................................

SUBROUTINE PRTDIS(ID,U,F) !PRTD 1 C.....PRINT DISPLACEMENTS !PRTD 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PRTD 3 DIMENSION ID(NDF,*),U(*),F(NDF,*) !PRTD 4 WRITE(6,2000)(K,K=1,NDF) !PRTD 5 NN=NJ*NDF !PRTD 6 NM=NEQ !PRTD 7 C.....FORM COMPLETE DISPLEMENT ARRAY !PRTD 8

DO N=NJ,1,-1 !PRTD 9

DO I=NDF,1,-1 !PRTD 10

K=ID(I,N) !PRTD 11

IF(K.LT.0) THEN !PRTD 12

U(NN)=F(I,N) !PRTD 13

ELSE !PRTD 14

U(NN)=U(NM) !PRTD 15

NM=NM-1 !PRTD 16

ENDIF !PRTD 17

NN=NN-1 !PRTD 18 END DO !PRTD 19 END DO !PRTD 20 C.....PRINT DISPLACEMENTS !PRTD 21

K=0 !PRTD 22 DO N=1,NJ !PRTD 23 WRITE(6,3000) N,(U(K+I),I=1,NDF) !PRTD 24 K=K+NDF !PRTD 25 END DO !PRTD 26 RETURN !PRTD 27 2000 FORMAT(//1X,'NODAL DISPLACEMENTS'/2X,'NODE',6(I8,'DISP')) !PRTD 28 3000 FORMAT(I6,6E12.4) !PRTD 29 END !PRTD 30 C.........................................................

SUBROUTINE PRTREAC(R) !PRTR 1 C.....PRINT NODAL REACTIONS !PRTR 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PRTR 3 REAL R(NDF,*),RSUM(6) !PRTR 4 NNEQ=NDF*NJ !PRTR 5 DO K=1,NDF !PRTR 6

RSUM(K)=0.0 !PRTR 7 END DO !PRTR 8 WRITE(6,2000) !PRTR 9 DO I=1,NJ !PRTR 10 DO K=1,NDF !PRTR 11

RSUM(K)=RSUM(K)+R(K,I) !PRTR 12

IF(ABS(R(K,I)).GT.1.0E-3) WRITE(6,3000) I,K,R(K,I) !PRTR 13 END DO !PRTR 14 END DO !PRTR 15 WRITE(6,4000)(RSUM(K),K=1,NDF) !PRTR 16 2000 FORMAT(/5X,'NODAL REACTIONS'/) !PRTR 17 3000 FORMAT(2X,'NODE',I3,I20,'DOF',F12.4) !PRTR 18 4000 FORMAT(3X,'SUM',6F12.4) !PRTR 19 RETURN !PRTR 20 END !PRTR 21 C..........................................................

SUBROUTINE ELEMLIB(D,U,X,IX,S,P,I,J,ISW) !ELEM 1 C.....ELEMENT LIBRARY !ELEM 2

REAL P(NST),S(NST,NST),U(1) !ELEM 3 DIMENSION IX(1),D(1),X(1) !ELEM 4 COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !ELEM 5 COMMON /ELDATA/ LM,N,MA,MCT,IEL,NEL,NST !ELEM 6 SELECT CASE(IEL) !ELEM 7 CASE(1) !ELEM 8 CALL TRUSS(D,U,X,IX,S,P,I,J,ISW) !ELEM 9

CASE(2) !ELEM 10 CALL PLANE(D,U,X,IX,S,P,I,J,ISW) !ELEM 11 CASE(3:) !ELEM 12 WRITE(6,4000) IEL !ELEM 13 STOP !ELEM 14 END SELECT !ELEM 15 RETURN !ELEM 16 4000 FORMAT(5X,' **FATAL ERROR** ELEMENT CLASS NUMBER',I5,' INPUT') !ELEM 17 END !ELEM 18 C............................................................

SUBROUTINE TRUSS(D,UL,XL,IX,S,P,NDF,NDM,ISW) !TRUS 1 C.... TRUSS LINEAR ELASTIC ELEMENT ROUTINE !TRUS 2

COMMON /ELDATA/ LM,N,MA,MCT,IEL,NEL,NST !TRUS 3 DIMENSION D(10),UL(NDF,1),XL(NDM,1),IX(1),S(NST,1),P(1),DU(2) !TRUS 4 C.... GO TO CORRECT ARRAY PROCESSOR !TRUS 5

SELECT CASE (ISW) !TRUS 6 C.... INPUT MATERIAL PROPERTIES !TRUS 7

CASE(1) !TRUS 8

READ(5,*) E,A !TRUS 9

D(1)=E*A !TRUS 10 WRITE(6,2000) E,A !TRUS 11 C.....CALCULATE ELEMENTAT STIFFNESS MATRIX !TRUS 12

CASE(2,3) !TRUS 13 DX=XL(1,2)-XL(1,1) !TRUS 14 DY=XL(2,2)-XL(2,1) !TRUS 15 DL=SQRT(DX*DX+DY*DY) !TRUS 16 EAL=D(1)/DL !TRUS 17 CC=DX/DL !TRUS 18 SS=DY/DL !TRUS 19 S(1,1)=CC*CC*EAL !TRUS 20 S(2,2)=SS*SS*EAL !TRUS 21 S(1,2)=CC*SS*EAL !TRUS 22 S(2,1)=S(1,2) !TRUS 23 IF(ISW.EQ.2) THEN !TRUS 24

DO I=1,NDF !TRUS 25

DO J=1,NDF !TRUS 26

S(I+2,J) = -S(I,J) !TRUS 27

S(I,J+2) = -S(I,J) !TRUS 28

S(I+2,J+2) = S(I,J) !TRUS 29

END DO !TRUS 30

END DO !TRUS 31 C.....CALCULATE ELEMENTAT NODAL FORCE !TRUS 32

ELSE !TRUS 33

DU(1)=UL(1,1)-UL(1,2) !TRUS 34

DU(2)=UL(2,1)-UL(2,2) !TRUS 35

DO I=1,NDF !TRUS 36

DO J=1,NDF !TRUS 37

P(I)=P(I)-S(I,J)*DU(J) !TRUS 38

END DO !TRUS 39

P(I+2)=-P(I) !TRUS 40

END DO !TRUS 41 C.....CALCULATE ELEMENTAT AXIAL FORCE !TRUS 42

FN=P(1)*CC+P(2)*SS !TRUS 43

WRITE(6,3000) N,IX(1),IX(2),FN !TRUS 44 END IF !TRUS 45 END SELECT !TRUS 46 2000 FORMAT(/5X,'PLANE TRUSS LINEAR ELASTIC ELEMENT'// !TRUS 47

1 10X,7HMODULUS,E18.5/10X,4HAREA,E21.5/) !TRUS 48 3000 FORMAT(/5X,'ELEMENT',I8,5X,'NODE1',I5,5X,'NODE2',I5,5X, !TRUS 49

1 'AXIAL FORCE',E15.6) !TRUS 50

RETURN !TRUS 51 END !TRUS 52 C............................................................

SUBROUTINE PLANE(D,UL,XL,IX,S,P,NDF,NDM,ISW) !PLAN 1 C.... PLANE LINEAR ELASTIC ELEMENT ROUTINE !PLAN 2

COMMON /ELDATA/ LM,N,MA,MCT,IEL,NEL,NST !PLAN 3 DIMENSION D(10),UL(NDF,1),XL(NDM,1),IX(1),S(NST,1),P(1), !PLAN 4

1 SHP(3,9),SG(16),TG(16),WG(16),SIG(6),EPS(3),WD(2) !PLAN 5

CHARACTER WD*4 !PLAN 6 DATA WD/'RESS','RAIN'/ !PLAN 7 C.... GO TO CORRECT ARRAY PROCESSOR !PLAN 8

SELECT CASE (ISW) !PLAN 9 C.... INPUT MATERIAL PROPERTIES !PLAN 10

CASE(1) !PLAN 11 READ(5,*)E,XNU,L,I !PLAN 12 D(4)=L !PLAN 13 IF(I.NE.0) I = 1 !PLAN 14 IF(I.EQ.0) I = 2 !PLAN 15 D(1) = E*(1.+(1-I)*XNU)/(1.+XNU)/(1.-I*XNU) !PLAN 16 D(2) = XNU*D(1)/(1.0E0+(1-I)*XNU) !PLAN 17 D(3) = E/2.0E0/(1.0E0+XNU) !PLAN 18 WRITE(6,2000) WD(I),E,XNU,L !PLAN 19 C.... STIFFNESS COMPUTATION !PLAN 20

CASE(2) !PLAN 21 L = D(4) !PLAN 22 IF(L*L.NE.LINT) CALL PGAUSS(L,LINT,SG,TG,WG) !PLAN 23 C.... FAST STIFFNESS COMPUTATION, COMPUTE INTEGRALS OF SHAPE FUNCTIONS !PLAN 24

DO L = 1,LINT !PLAN 25

CALL SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX) !PLAN 26

XSJ = XSJ*WG(L) !PLAN 27 C.... LOOP OVER ROWS !PLAN 28

J1 = 1 !PLAN 29

DO J=1,NEL !PLAN 30

W11 = SHP(1,J)*XSJ !PLAN 31

W12 = SHP(2,J)*XSJ !PLAN 32 C.... LOOP OVER COLUMNS (SYMMETRY NOTED) !PLAN 33

K1 = J1 !PLAN 34

DO K = J,NEL !PLAN 35

S(J1 ,K1 ) = S(J1 ,K1 ) + W11*SHP(1,K) !PLAN 36

S(J1 ,K1+1) = S(J1 ,K1+1) + W11*SHP(2,K) !PLAN 37

S(J1+1,K1 ) = S(J1+1,K1 ) + W12*SHP(1,K) !PLAN 38

S(J1+1,K1+1) = S(J1+1,K1+1) + W12*SHP(2,K) !PLAN 39

K1 = K1 + NDF !PLAN 40

END DO !PLAN 41

J1 = J1 + NDF !PLAN 42

END DO !PLAN 43 END DO !PLAN 44 C.... ASSEMBLE THE STIFFNESS MATRIX FROM INTEGRALS AND MATERIAL PROPS. !PLAN 45

NSL = NEL*NDF !PLAN 46 DO J = 1,NSL,NDF !PLAN 47

DO K = J,NSL,NDF !PLAN 48

W11 = S(J,K) !PLAN 49

W12 = S(J,K+1) !PLAN 50

W21 = S(J+1,K) !PLAN 51

W22 = S(J+1,K+1) !PLAN 52

S(J ,K ) = D(1)*W11 + D(3)*W22 !PLAN 53

S(J ,K+1) = D(2)*W12 + D(3)*W21 !PLAN 54

S(J+1,K ) = D(2)*W21 + D(3)*W12 !PLAN 55

S(J+1,K+1) = D(1)*W22 + D(3)*W11 !PLAN 56 C.... FORM LOWER PART BY SYMMETRY !PLAN 57

S(K,J) = S(J,K) !PLAN 58

S(K,J+1) = S(J+1,K) !PLAN 59

S(K+1,J) = S(J,K+1) !PLAN 60

S(K+1,J+1) = S(J+1,K+1) !PLAN 61

END DO !PLAN 62 END DO !PLAN 63 C.... COMPUTE ELEMENT STRESSES, STRAINS AND FORCES !PLAN 64

CASE(3) !PLAN 65 L = D(4) !PLAN 66 IF(L*L.NE.LINT) CALL PGAUSS(L,LINT,SG,TG,WG) !PLAN 67 DO L = 1,LINT !PLAN 68 C.... COMPUTE ELEMENT SHAPE FUNCTIONS !PLAN 69

CALL SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX) !PLAN 70 C.... COMPUTE STRAINS !PLAN 71

DO I = 1,3 !PLAN 72

EPS(I) = 0.0 !PLAN 73

END DO !PLAN 74

XX = 0.0 !PLAN 75

YY = 0.0 !PLAN 76

DO J = 1,NEL !PLAN 77

XX = XX + SHP(3,J)*XL(1,J) !PLAN 78

YY = YY + SHP(3,J)*XL(2,J) !PLAN 79

EPS(1) = EPS(1) + SHP(1,J)*UL(1,J) !PLAN 80

EPS(3) = EPS(3) + SHP(2,J)*UL(2,J) !PLAN 81

EPS(2) = EPS(2) + SHP(1,J)*UL(2,J) + SHP(2,J)*UL(1,J) !PLAN 82