关 键 词：

钢结构 有限元分析

资源描述：

钢结构有限元分析,钢结构,有限元分析

内容简介：

破碎站钢结构有限元分析,主要设计工作,1建立破碎站钢结构有限元模型2有限元模型的模态分析3有限元模型静力学分析（应变和应力）,1在ANSYS9.0软件中建立有限元模型,破碎站的主视图,受料仓与给料机,受料仓和给料机共686个单元，其中梁单元598个，杆单元88个，节点总数为597个,破碎平台与控制室,破碎平台和控制室共1055个单元，其中梁单元999个，杆单元56个，节点总数1008个，,模态分析,表4.1受料仓与给料机各阶固有频率(单位：HZ),表4.1破碎平台与控制室各阶固有频率(单位：HZ),摘 要本文主要对某煤矿地面生产系统，一次破碎站钢结构进行有限元分析。破碎站由受料仓与给料机和破碎平台与控制室两部分组成。对两部分的钢结构分别进行有限元分析。在结果中找到危险的部位进行具体的分析。首先，建立受料仓与给料机的有限元实体模型。计算等效的载荷，计算出钢结构在载荷下的应力和变形并分析它们的分布情况。其次，破碎平台与控制室求解过程和上边的一样，但是破碎平台和控制室的连接是铰接，所以在建模的过程中采用耦合的方法进行处理。最后，对两个有限元实体模型进行模态分析，分别求解出固有频率和模态振型图。关键词 有限元；钢结构；模态分析ABSTRACTThis dissertation mainly to an open coalmine ground production system, one broken to stand steel construction finite element analysis. Store -give material machine and broken platform- control room two parts make up the crush station. Finite element analysis to the steel construction of two parts comparatively. Find the dangerous part to carry on concrete analysis of the result.First of all, set up the finite element of Store -give material machines entity model. Calculate the equivalent load; solve out the stress and strain of the steel construction under the load and analysis their distribution situation.The next place, the course of solving is the same as above. But the connections of the broken platform and control room are the hinged joint, so deal with by coupling in the course of modeling. Finally, carry on mode analysis to two finite element entity models; it is solve the intrinsic frequencies and mode picture of shaking, respectively.Keyword finite element；steel construction；mode analysis目 录中文摘要英文摘要1 前言11.1有限元分析方法介绍11.2大型有限元分析软件ANSYS介绍21.3主要工作32 受料仓与给料机的钢结构有限元分析42.1建立有限元模型42.2载荷等效计算62.2.1主要结构截面几何参数62.2.2实际载荷情况72.2.3实际等效计算结果72.3有限元分析结果102.3.1受料仓与给料机整体位移102.3.2分析部位图122.3.3支撑立柱结果132.2.4两根纵梁结果173 破碎平台与控制室的钢结构有限元分析193.1建立有限元模型193.2载荷等效计算223.2.1主要结构截面几何参数223.2.2破碎平台实际载荷情况233.2.3破碎平台实际等效计算结果243.3有限元分析结果263.3.1破碎平台与控制室整体位移263.3.2顶层横梁结果273.2.3破碎机支撑梁结果263.2.4破碎机立柱结果294 破碎站钢结构模态分析314.1受料仓与给料机的固有频率和振型图314.2破碎平台与控制室的固有频率和振型图32参考文献35致 谢36英文资料原文英文资料翻译IV2 受料仓与给料机的钢结构有限元分析2.1建立有限元模型如图2.1破碎站主视图和图2.2破碎机布置图，它的工作过程是：卸料卡车间歇把最大入料粒度为1500mm的煤块倒入受料仓，受料仓存储大粒度煤块。刮板给料机把受料仓的大粒度的煤块连续的刮给破碎平台的破碎机。破碎机把最大入料粒度为1500mm的煤块破碎成最大排料粒度为300mm的煤块，煤块由底部的传送带传出。图2.1 破碎站主视图图2.2 破碎机布置图破碎站钢结构的弹性模量E200000MPa，泊松比m=0.3，质量密度r=7.810-3kg/cm3。破碎站由支撑件H型钢和斜支撑(角钢)组成。在结构离散化时，由于角钢和其它部位铰接，铰接是具有相同的线位移，而其角位移不同。承受轴向力，不承受在其它方向的弯矩，相当于二力杆，所以H型钢用梁单元模拟，角钢用杆单元模拟。破碎站是由受料仓与给料机和破碎平台与控制室两部分组成，故计算时是分别对这两部分进行的。离散后，受料仓和给料机共686个单元，其中梁单元598 个，杆单元88个，节点总数为597个，有限元模型如图2.3和图2.4所示。图2.3 受料仓与给料机有限元模型图2.4 受料仓与给料机有限元模型俯视图2.2载荷等效计算2.2.1主要结构截面几何参数破碎站主要结构采用H型钢梁，截面尺寸如图2.5所示，各截面横截面积A，截面惯性矩Iy，Iz和极惯性矩I如下。图2.5 截面尺寸料仓及给料机支撑结构料仓及给料机六根支撑立柱(H5004001220)A= 215.2mm2，Iy101947104mm4，Iz21340104mm4，I240104mm4料仓BB面横梁和给料机EE、FF面横梁(H4003001220)A=16320mm2，Iy48026104mm4，Iz9005104mm4，I181104mm4料仓CC面和DD面横梁(H4004001220)A=20320mm2，Iy 62479104mm4，Iz21339104mm4，I234104mm4给料机两根纵梁(H5504001220)A=22120mm2，Iy125678104mm4，Iz21341104mm4，I243104mm4给料机六根横梁(H4004001220)A=20320mm2，Iy 62479104mm4，Iz21339104mm4，I234104mm4其它横梁(H4003001220)A=16320mm2，Iy 48026104mm4，Iz 9005104mm4，I181104mm4斜支撑的横截面积12512：A2856mm275 6：A864mm22.2.1实际载荷情况给料机自重载荷：65000kg最大驱动扭矩：2*150kN-m给料机动力载荷：垂直载荷系数:1.3；水平载荷系数:0.3受料仓支撑柱实际支撑物料载荷：150000kg给料机尾部受仓压载荷：50000kg给料机前部物料载荷：4800kg/m物料冲击载荷（或偏心弯矩载荷）：12000kg-m(作用于每个仓柱上)走道平台载荷:沿设备两侧考虑走道宽度为1米，按300kg/m2考虑结构风荷载：按0.6kN/m2考虑地震载荷：地震烈度8度料仓自重载荷：85000kg清扫溜槽自重：6000kg导料挡板自重载荷：30000kg支撑结构（除滑橇自重）:30000kg2.2.2实际等效计算结果走道平台载荷数值：q13.6(N/mm)位置：垂直作用在给料机两根纵梁上。风载数值：q22.295(N/mm)位置：-z方向的迎风梁上。 驱动扭矩数值：M195106(Nmm)位置：纵梁的前端，由破碎平台传递过来。料机尾部受仓压载荷数值：q369.02899(N/mm)位置：纵梁的尾部，方向为重力方向。给料机前部物料载荷位置：纵梁的E_E截面到给料机前端。数值：q467.200(N/mm)导料挡板自重数值：q512.85714(N/mm)位置：垂直作用在给料机两根纵梁上。清扫滑溜槽自重数值：q62.57143(N/mm)位置：垂直作用在给料机两根纵梁上。图2.6 作用在受料仓与给料机上的压力载荷模型受料仓支撑柱支撑物料载荷数值：F300000(N)位置：受料仓每根支撑柱的顶端，方向为重力方向。物料冲击载荷数值：M28.000(Nmm)位置：受料仓每根支撑柱的顶端，方向由右手定则判断。给料机自重载荷Fy153833(N)位置：受料仓每根支撑柱的顶端，方向为重力方向。Fz45500(N)位置：受料仓每根支撑柱的顶端，其水平载荷方向与风载相同。图2.7 作用在受料仓与给料机上的集中载荷模型图2.8 受料仓与给料机的载荷模型2.3有限元分析结果2.3.1受料仓与给料机整体位移 最大z方向的位移出现在受料仓连结BCD面的水平横梁上，其值为-4.272mm，如图2.9所示，主要是由物料冲击载荷和给料机自重载荷的水平分量引起的。图2.9 受料仓与给料机的z方向变形最大y方向位移位于给料机纵梁最前端，如图2.10所示，其中y方向位移为8.531mm，主要是由刮板给料机的驱动扭矩和给料机前部的物料载荷共同引起的。图2.10 受料仓与给料机的y方向变形如图2.11所示，x方向位移为4.492mm，最大等效位移在给料机纵梁最前端，位移为4.492mm，这里属于悬臂梁，虽然有大截面的斜支撑，但是却是要重点分析的部位。图2.11 受料仓与给料机的x方向变形如图2.12受料仓与给料机的整体变形，其中y方向的分量占的比重最大，它的方向和主要载荷在同一个方向。其它方向的分量所占的比重比较小。整体变形的最大数值为9.649mm。图2.12 受料仓与给料机的整体的变形2.3.1分析部位图为了便于对出现危险部位位置的描述，参考图2.13危险部位侧视图1 2 3 4 5 6 7 . 8 9 10 11图2.13 危险部位侧视图1BB截面2纵梁尾部3CC截面4斜支撑5DD截面6纵梁中部7EE截面8FF截面9斜支撑10斜支撑和纵梁的铰接处11纵梁前端如图2.13在y方向为实际物体的重力方向，也是立柱的方向。X方向为纵梁的方向。纵梁与水平面上倾10度。在垂直纸面的z方向有横梁连接立柱和纵梁。2.3.2支撑立柱结果 轴向力如图2.14所示。由图可见，支撑立柱受到轴向压力作用，EE面的两根立柱轴向力较小，FF面两根立柱轴向力大，最大轴向力作用FF面立柱底部至第一根水平横梁之间，其值为-497590N。图2.14 立柱轴向力图轴向应力如图2.15所示。对应的最大轴向应力为-26.034MPa。因为这几根立柱的截面几何参数一样，所以出现的位置与轴向力相同。图2.15 立柱轴向应力图相对于立柱梁单元局部坐标y轴的弯矩如图2.16所示，EE面两根立柱弯矩较大，最大弯矩位于EE面立柱顶端，最大值My0.993108Nmm，最小弯矩位于CC面立柱顶端，最小值My-0.995108Nmm，图2.16 局部坐标y轴的弯矩图对应的应力如图2.17所示，最大应力为33.031MPa。最小应力为-33.118MPa。图2.17 局部坐标y轴的弯曲应力图相对应立柱梁单元局部坐标z轴的弯矩如图2.18所示，最大弯矩位于BB面中风载作用面的立柱，底部最大弯矩0.126108Nmm，在BB面拉筋连结点处弯矩数值为-0.196108Nmm，DD面第一根水平横梁处弯矩为0.25563108Nmm，图2.18 局部坐标z轴的弯矩图对应的应力如图2.19所示，BB面中风载作用面立柱底部应力为19.625MPa，底部连结点处应力为30.605MPa，DD面第一根水平横梁处应力为24.024MPa。图2.19 局部坐标z轴的弯曲应力图2.2.3两根纵梁结果 两根纵梁轴向力如图2.20，可见两根纵梁轴向应力很小，最大轴向力192337N，位于F面和纵梁与斜支撑的接触之间。这里的变形也是最大的。轴的方向与大部分的载荷方向近似垂直。在斜支撑与纵梁连接到纵梁的前端只承受弯矩，不承受轴向力。图2.20 两根纵梁轴向力图两根纵梁轴向应力如图2.21所示，可见两根纵梁轴向应力很小，最大轴向应力13.066MPa，位于F面和纵梁前端之间。图2.21 两根纵梁轴向应力图相对于纵梁梁单元局部坐标y轴的弯矩如图2.22所示，其最小弯矩My-0.176109Nmm，位于斜支撑与纵梁连接处，这里的弯矩最大驱动扭矩作用在悬臂梁顶端。图2.22 局部坐标y轴的弯矩图对应的应力如图2.23所示，最小应力-97.755MPa，位于FF面处，。图2.23 局部坐标y轴的弯曲应力图相对于纵梁梁单元局部坐标Z轴的弯矩如图2.24所示，其最大弯矩Mz0.174108Nmm，这里是承受给料机尾部受仓压载荷，位置是纵梁和C_C截面相连接。图2.24 局部坐标Z轴的弯矩图对应的应力如图2.25所示，最小应力-47.911MPa。位于CC面处。图2.25 局部坐标Z轴的弯曲应力图193 破碎平台的钢结构有限元分析3.1建立有限元模型破碎站图纸是二维图纸，在建立有限元模型的过程中充分应用了并行过程，料仓与给料机和破碎平台与控制室两部分是基本上同时建立起来的，破碎平台和控制室两部分是属于同一个模型，应该分别建立。由于承载核心部件破碎机的是破碎平台，所以先从破碎平台开始建立模型，再定义材料、单元类型、截面形状还要兼顾后续进程。 图3.1 为破碎站右视图，一共有四层的结构：第一层，放置传送带，传送带把破碎机破碎下来的物料运走，第二层，放置核心部件破碎机器，第三层和第四层，主要是放置破碎机。最顶层有控制室和起重机。图3.1 破碎站右视图破碎站钢结构的弹性模量E200000MPa，泊松比m=0.3，质量密度r=7.810-3kg/cm3。破碎站由支撑件H型钢和斜支撑(角钢)组成。在结构离散化时，由于角钢和其它部位铰接，铰接是具有相同的线位移，而其角位移不同。只承受轴向力，而且无其它方向的弯矩，所以H型钢用梁单元模拟，角钢用杆单元模拟。破碎站是由料仓与给料机和破碎平台与控制室两部分组成，故分别对这两部分进行计算。离散后，破碎平台和控制室共1055个单元，其中梁单元999个，杆单元56个，节点总数1008个，整体限元模型如图3.2所示。 图3.2 破碎平台和控制室有限元模型这个模型的破碎平台是主要承载，破碎机自重载荷、最大驱动扭矩、物料冲击载荷、上罩和下溜槽自重载荷和破碎机动力载荷等一系列重要载荷的重要部位，采用的支撑件H型钢和斜支撑(角钢)的面积都是最大的，所以在建立模型的过程中，控制单元长度大致是控制室的一半，这样保证重要部位的精度，同时也符合网格划分的粗精结合的基本原则。控制室采用的支撑件H型钢和斜支撑（角钢）的面积都是较小的，主要承载，起重设备及起重载荷、控制室载荷、走道平台载荷、结构风荷载等。破碎站钢结构立柱底部与基础相连，故立柱底部按固定端处理，即每根立柱底端的6个自由度都加以约束。破碎平台与上部控制室采用铰接连结，故破碎平台与上部控制室在模型中采用特殊的方法进行连接。破碎平台的重要部分如下图3.3 破碎平台的有限元模型，共408个单元，420个节点。在建立模型的过程中H型钢和斜支撑(角钢)的连接处的节点坐标，应该进行必要的记录和严格的控制。最忌讳到后期要建立新的连接的时候，发现没有合适的节点。因为当单元部分建立好的时候，后加节点，再改单元，修改的工作量将十分巨大。尤其在加载荷以后，加在单元的压力载荷会因为单元号的改变而移动位置，严重情况会导致加的载荷难以再修改，再进行重新加载荷非常麻烦。图3.3 破碎平台的有限元模型3.2载荷等效计算3.2.1主要结构截面几何参数破碎站主要结构采用H型钢梁，截面尺寸如图31所示，各截面横截面积A，截面惯性矩Iy，Iz和极惯性矩I如下。图3.4 截面尺寸破碎平台四根支撑立柱(H5004001220)A=21520mm2，Iy101947104mm4，Iz21340104mm4，I240104mm4破碎平台横梁(H7003001220)A=19920mm2，Iy16750104mm4，Iz9010104mm4，I198104mm41.2. 2控制室立柱(H300300812)A =9412mm2， Iy16340104mm4，Iz5401104mm4，I39104mm4水平横梁(H300200812)A =7010mm2， Iy11361104mm4，Iz1601104mm4，I28104mm4斜支撑的横截面积12512：A2856mm2 75 6：A864mm2 90 6：A1044mm23.2.1 破碎平台结构载荷主要结构断面形式立柱断面：H5004001220水平梁断面：H7004001220、H7003001220斜支撑： 212512510、290906载荷情况破碎机自重载荷：70000kg破碎机动力载荷：垂直载荷系数：电机侧1.8；非电机侧1.55水平载荷系数:0.25最大驱动扭矩：50kN-m；（过载扭矩）：300kN-m物料冲击载荷（或堵料载荷）：12000kg走道平台载荷：按500kg/m2考虑结构风荷载：按0.6kN/m2考虑地震载荷：地震烈度8度上罩和下溜槽自重载荷：15000kg支撑结构自重载荷：25000kg控制室及起重维护结构载荷主要结构断面形式立柱断面：H300300812水平梁断面：H300200812斜支撑：290906载荷情况起重设备及起重载荷：10000kg（其中设备5000kg；起重量5000kg）控制室载荷：8000kg走道平台载荷；按300kg/m2考虑结构风荷载：按0.6kN/m2考虑地震载荷：地震烈度8度支撑结构自重载荷：15000kg3.2.2破碎平台实际等效计算结果风载荷数值q1y0.789498(N/mm)位置：-z方向的迎风梁。数值q1x0.84872(N/mm)位置：-x方向迎风梁。破碎机自重载荷数值：q2z182.20588(N/mm)位置：非电机侧支撑梁。数值：q2z218.23529(N/mm)位置：电机侧支撑梁。数值：q2x36.02941(N/mm)位置：方向与风载荷相同，水平作用在支撑梁。物料冲击载荷数值：q324.70588(N/mm)位置：破碎机支撑梁。上罩和下溜槽自重载荷数值：q46.47059(N/mm)位置：破碎机支撑梁。起重设备及起重载荷数值：F35000(N)位置：控制室顶层横梁。控制室载荷控制室支撑梁，位置：沿重力方向，在控制室支撑梁与控制室结合处。 数值q524.000(N/mm)破碎平台的走道平台载荷。 数值：q64.96671(N/mm)位置：控制室顶层横梁。控制室走道平台载荷第一层走道平台长其载荷。 数值：F13539.6(N)位置：如图3.5。控制室第二层和第三层平台长相同。数值：F12981.6(N)位置：如图3.5。控制室顶层走道平台载荷。 数值：q19.440(N/mm)位置：如图3.5。图3.5 破碎平台的载荷图3.3有限元分析结果3.3.1破碎平台与控制室整体位移分析变形是材料力学中一种重要的应变分析（还包括转角等），位移计算结果如下。由于风载荷在y-z平面的作用面积很大，在顶层的z向梁和其它两个方向的梁相比较较长，最大位移位于结构顶层，位移值为7.233mm，在破碎平台顶层的两根梁，变形也较大，因为它承载起重设备及起重载荷、控制室载荷、走道平台载荷、-z向的结构风荷载等，如图3.6所示。 图3.6 破碎平台和控制室整体变形x方向最大位移为-7.19mm，出现位置和最大挠度出现的位置相同；由于y方向结构刚度较大，故y方向位移较小，其最大值为-4.574mm，位于破碎平台顶层的那两根上，主要由起重设备及起重载荷、控制室载荷、走道平台载荷引起；z方向最大位移为-1.426mm，位于破碎平台的横梁上，它是由风载荷引起的，所以最小。3.3.2顶层横梁结果通过观察整体的弯曲应力图发现由沿着工字钢中间平面弯矩引起的最大应力部位为位于破碎平台顶层的那两根上，即与y方向的挠度相同。它的弯矩如图3.7，它的左右两边的不对称性是由于风载荷引起的，最小值是-0.4108 N/mm。图3.7控制室顶层横梁局部坐标y轴的弯矩图相对应的应力为,最大值为39.697 MPa。最小的值是-52.812 MPa，出现的位置与对应的弯矩出现的位置相同。图3.8 控制室顶层横梁局部坐标y轴弯曲应力图3.2.3破碎机支撑梁结果对于其底部的两根承载破碎机器的那两根梁，承载找主要的载荷，虽然有同截面的若干根梁焊接上，并且有大截面的斜支撑，即T字梁进行连接，以保证大截面工字钢的四边型稳定性。它的弯矩是0.309108N/mm。图3.7控制室顶层横梁局部坐标z轴的弯矩图最危险的位置是在与垂直方向上相接触的那一点上，这里的要接触两个大型工字钢和四个斜支撑，这样更容易产生应力集中。由于结构的对称性，两根钢梁的最值相同。所以校核一根就可，它的应力也是出现在相同的位置，最小的应力是-51.525MPa，最大价值35.509MPa。 图3.8 控制室顶层横梁局部坐标z轴弯曲应力图3.2.4破碎机立柱计算结果工字钢的轴向拉伸也是应该考虑的，最大的拉力是在迎风面的控制室的铅直梁上，它的最大轴向力，如图3.11是84919N,对应的最大轴向应力， 图3.11 破碎平台立柱轴向力图如图3.12为，9.026 MPa。出现的位置为控制室的铅直梁与破碎平台铰接处。最大的应力为12.7025MPa。它的值不是很大，但是由于是铰接处，对于分析还是很有意义的。 图3.128 破碎平台立柱轴向应力图304 破碎站钢结构固有频率分析4.1受料仓与给料机的固有频率和振型图模态分析用于确定钢结构的固有频率。这使设计工程师们可以避开这些频率或最大限度地减小对这些频率上的激励，从而消除过度振动。振动模态分析用于求出结构自然频率和模态形状，也称固有频率和主振型。该分析的结果对于实际工程设计有关参数的选择（如激振频率的确定、共振现象的避免与利用等）及进一步的动力分析都很重要，因为结构的基本频率和模态信息能够反映动态响应特性，所以先求出受料仓与给料机各阶的固有频率，如表4.1，其中，不同阶数数值相同的是复根。模态分析主要用于决定结构的固有频率和振型。这是结构承受动态载荷的重要数据，同时，模态分析也是其他动力学分析的起点，例如瞬态动力学分析、谐响应分析和谱分析的起点。表4.1 受料仓与给料机各阶固有频率 (单位：HZ)阶 数频 率阶 数频 率115.6551125.797215.6811225.800316.5861325.832417.1231426.741521.1451527.211621.6601627.578721.6821727.604821.8791828.756923.6331928.7581025.7952028.763本节通过对受料仓与给料机的钢结构的有限元动力学模型的模态分析，算出了该受料仓与给料机的钢结构的前20阶固有频率和前3阶振型图。通过振型图和动画显示，可以很直观地分析受料仓与给料机的钢结构的动态性能，为受料仓与给料机的钢结构设计提供理论依据。前3阶振型图。图4.1 受料仓与给料机的第1阶振型图图4.2 受料仓与给料机的第2阶振型图图4.3受料仓与给料机的第3阶振型图4.2破碎平台与控制室的固有频率和振型图对钢结构的有限元模型的求解，一般不需要求出振动系统的全部固有频率和振型，由于低阶模态对钢结构振动系统的影响较大，因此本文仅计算了前20阶模态。振动模态分析用于求出结构自然频率和模态形状，也称固有频率和主振型。如表4.1。表4.1 破碎平台与控制室各阶固有频率 (单位：HZ)阶 数频 率阶 数频 率15.60011111.07327.10651211.10239.43841311.226410.0211412.181510.1961512.285610.8881612.287711.0391712.287811.0471812.287911.0641912.2881011.0662012.288图4.4破碎平台与控制室的第1阶振型图4.5破碎平台与控制室的第2阶振型图4.6破碎平台与控制室的第3阶振型34一种适应的有限元分析方法;面向一个整体计算环境摘要 这个工作介绍一个自适应数字程序的方法,它依靠变量构成整体,面向目标, 计算环境包括分析的前处理和后处理建模。一个为数据实验和后续发展准备的基本平台,它允许新的单元/误差估计和灵敏度分析工具。一个整体的高级收敛碎片恢复(SPR)和最近用平衡计划碎片中现场恢复(REP)的工具。SPR和REP进行比较并且应用误差分析和自适应引导从新建造网格程序。然而, SPR应用在第一计算灵敏数量和高级命令。网格(从新)生成过程是应用现代的集合的嵌块与图谱方法手段和Delaney三角技术。表面网格生成在任意范围自动进行(例如。在无使用权干涉)在不是三角单元就是四边单元自适应分析。这些思想在有限元自适应软件系统技术工具。效力和多功能的FESTA是表示用典型数字举例说明有限元的相互联络分析,恢复程序,误差分析/自适应和自动网格生成。检索词 有限元分析,误差估计, 自适应,高端精密, 灵敏度, 高级收敛碎片恢复(SPR), 平衡计划碎片中现场恢复(REP),原目标程序(OOP),相互电脑绘图。介绍这个工作出现在完整的(原目标) 自适应计算分析的计算环境一般性的二维(2D)问题。这个环境包括确保给一个准确水平根据一定标准,也有分析程序去去产生和修改分离有限单元。这个计算系统,命名FESTA(自适应有限元系统技术),包含五个组成部分(参见图表1箱体颜色): 一个图表前处理,为定义几何问题,最初的有限元网格(包括边界条件),和自适应分析中应用主参数。这些几何模型是与有限员模型分离的。 为解决当前的边界准则有限元模型问题。开发编码为了更好的模块化,可扩展性,用户方便性。这样,它能更适合维修和升级。然而,其它的用户/开发人员 应该有能力修改基本程序系统去适合他们自己专门应用。误差分析和灵敏度模块。区分误差是根据估计可得到的重新覆盖程序，例如，高级收敛碎片恢复(SPR)和平衡计划碎片中现场恢复(REP)。各种灵敏度命令（一级，二级或者更高级）计算处理手段类似于SPR。用户选择想要的误差估计和灵敏度命令。 一个网格（重新）生成（而不是网格细化）程序，基于集合嵌块与图谱 和Delaney三角技术。根据放大误差，计算在前期模型，新的有限元网格在自动生成（例如：无用户干涉）,不是用三角单元就是用四边单元（高精度）。 最后,后处理模型,所有这些分析结果(例如外型破坏,灵敏度和受压轮廓)可视操作。有限元求解图表预处理（几何拓扑）网格循环设计流程图最终分离是误差分析(ZZ, SPR, REP)是否收敛否重新生成网格形象化（后处理）图形-简易相互作用网格示意图实质上,FESTA是一个提供基本平台为数值分析和后续开发例如,新的误差分析工具,单元,或材料模型,估计计算实验室。导向目标程序和完整的预处理分析和用做FESTA软件进行后处理适合工程实践应用和后续研究开发。这片论文的其余部分如下组织。动力工作和文献摘要考察提供在第二段以后,第三段出现一些在自适应模拟理论背景和用在FESTA软件的总揽接触面图表。第四段引入SPR(用于底重量调整系统),REP,和敏感度模型的数学公式。一个关于自动网格生成系统技术应用在这个工作的讨论在第五段相关的信息关于FESTA工具出现在第六段,尤其在SPR和REP恢复关系方面。为了达到计划计算体系的效能,典型数据举例在第七段给出,在第八段,推论出结论和讨论未来研究的方向。2相关研究的动向正常解决工程实际问题用有限元模型(FEM)或这边界单元模型(BEM)包含增加计算范围内离散点的数量和求解结果系统等价相关的改变在数值求解。一般的,这个程序消耗时间,它决定于分析者的经验,和求解可能被误导进入一个无任何症状的求解范围。理想的,用一个正确的和可靠的自适应方案,对象将能够指定最初离散模型,它可以充分的描述主要的几何/拓扑和边界,在指定适合的公差值,根据一个适当的标准。这样,这个系统将能够自动优化模型直到误差下降到低于指定的公差。这个过程将是自动进行的并且没有任何用户的干预。这个主要的目标是开发FESTA的动力因素。这种手段自从它不在依赖分析元的经验,非常规经验,从而增加整体分析的可靠性。FEM需要发展更好的预处理技术,为执行自动分析,和为获得自适应求解(它以成为一种FEM商业软件趋势)促进开发一种自动网格生成算法,例如不用任何用户的干预这种算法能离散任何几何导入有限元网格。一些二维的几何算法可以进行开发(例如Baehmann在1987年; Blacker and Stephenson 在1991年; Zhu 在。 1991;Potyondy 在1995年; Borouchaki 和 Frey 在1998年),三维几何模型手段最近出现(例如 Cass 在 1996年; Escobar 和Montenegro 在1996年; Beall 在。 1997年; Lo 在1998年)。目前这些工作集中在二维网格与自适应求解联合生成。全四边型和全三角型网格有效技术在细节上进行考虑。尽管目前在此处算法能扩展最大的网格,例如四边型和三角型网格单元(可以看见,例如, Borouchaki 和 Frey 在1998年),论题不是在这个工作的范围。目前有大量的著作关于误差分析和自适应,读者可以直接去参见参考文献1。 Blebbier 和 Aliabad (1993) 和 Babushka 。 (1986)编辑大量的FEM和BEM检查自适应的技术。 Ladeveze 和 Oden (1998)编辑的书收集从工厂的论文自适应计算技术的新进展发表在Cachan期刊,法国,1997年9月17-19期,它收录最新的进展在自适应技术方法和在工程实际的应用。一些其他杂志也有题献出自适应,例如:计算工程1996年12卷第二页, 工程软件进展1992年15卷四分之三页数。和计算和数学应用1991年36卷第一页。 调查发现Moor 和Babushka (1987), Odem 和Demoniacs (1989), Strouboulis和 Hague (1992a, b), Babushka 和 Sari (1994), 和 Ainsworth 和 Olden (1997)。 Mackerel (1993, 1994)编辑大量关于FEM产生网格文献,优化,误差分析和FEM 和BEM的自适应技术那些出版从1990到1993年。 Babushka et al。 (1983) 编辑出版 FEM 和差分有限模型（FDM）自适应技术。关于FEM的书最近强调在自适应求解技术领域。举例，Zienkiewicz and Taylor (1989)写一本书,其中有误差分析和自适应一章（14章），它是对Zienkiewicz and Zhu (1992a, b, 1994)的论文进行补充。然而，Seaboard and Babushkas (1991)是最初对这个课题进行研究。第一篇关于自适应有限元的论文出现在17世纪早期。从那以后，大量关于这个课题的论文出现在技术文献。Babushkas and Rheingold (1978)出现最先的论文关于留数求值近似求解的误差分析并且应用获得局部的，更准确的答案。他们开发自适应数学基础技术。用后期误差的观念，它是开发对FEM的自适应战略这样就会有一部分的有限单元得到优化出现自适应方法体系。在上世纪八十年代早期，电脑制图技术开始被应用在网格生成技术的标准的工具。Sheppard(1986)发布了一篇关于自动网格生成技术应用在,自适应连接方法的论文。Zienkiewiczand Zhu (1987) 引入进步在求解可覆盖的过程梯度（压力）变化曲线的值的误差分析。在有限元上用原代码转化成工具实现非常简单，这种类型的技术，是基于平均求值和命名为L2 的工程，都应用可覆盖的过程梯度变化曲线技术，可行性分析也取得一定成绩。在1992年，这项技术被同一个作者矫正/改进引导出高级收敛碎片恢复(SPR) (Zienkiewicz and Zhu 1992 a, b, 1994)。这种方法是当一片受压光滑的单元并且最小的方型进行离散。以其适应局部更高级的命令高级收敛试样从有限元的多项式中进行收索。各种对改进覆盖过程的尝试可以在例如;Weinberg et al。 (1994), Weinberg and Abdul Ahab (1993), Blacker and Belytschko (1994), Labara et al。 (1994), and Lee et al。 (1997)。的文献中被找到。基本上这些改进技术合并于目前的平衡方程技术和覆盖过程边界条件。通过Babushkas et al。 (1994a, b)进行广泛的研究发现，通过数据举例，SPR对于过余数手段表现出更加的优越性能。最近，Boomed and Zienkiewicz (1997)发表满足在较弱结构平衡条件的新的高级收敛方法。它不要求任何收敛点的相关知识。这种新的覆盖技术被命名为平衡计划碎片中现场恢复(REP)。SPR 和 REP都是目前特别主张的工作方法。 以上的贡献，自适应领域一般是显著的和最近有发展成绩的。举例，Paullina et al。 (1997) 基于节点敏感性的观念建立了误差估计新类别。它被应用在一般意图的计算方法的连接比如FEM, BEM 和 FDM。 Raunchier and Steelier (1997)建立了在FEM中的误差反馈控制手段。Mahtomedi 和 Konkani (1998)用拉力能量平衡理论建立自适应程序。更进一步，Lakeview androgen (1998)最近编辑并出版了自适应计算机制进展的总结性文章。在一定数量的质量模型中相互作用，相关性，在初期的分析中非常重要。这个理由是优化的估算方法，例如FEM，或显现方法。单元自由方法。均衡的 Galleria BEM (Paullina和Gray 1999), or边界节点方法 (Mukherjee 和Mukherjee 1997)。 在现代的综合技术环境下相关误差分析的观念和自适应在FEM中是目前集中的工作。3计算理论方面任何时间计算模型应用在求解管理不同的边界值均衡求解问题。误差在离散问题被引入,它减少连续数学模型自由有限数据问题。这种离散的误差被定义为精确求解和近似数值的差值。定义局部误差是一种精确值和近似求解数值的相差的测量手段。定义局部误差是类似于响应数量(例如置换)是一种典型的求解过程。定义局部误差等于精确值和近似求解数值的差值。自适应方法是数值计划,它自动调节自己改善自己的求解过程以其求得特别精确的值。着两种基本的自适应方法的组成,是误差分析和自适应策略。这些组成将在下被讨论。总体上,有两种的典型离散误差分析;优先法和后期处理法。尽管优先法误差分析是非常准确的对于较差条件特别求解等级问题。但是他们经常不能提供给的模型的准确误差分析。后期处理法应用信息误差估计获得技术在求解过程。,增加一些求解的先期误差估计求解的估算。它经过这里的采用使得离散误差分析可以提供更加精确数值。在自适应战略的内容。扩展的方法已经被其他的手段提供(例如双精度和互补方法)是目前集中的工作。这些方法包括h-,p-和重复扩展。这个计算工具关于h-,p-和重复扩展是分别独立的。在重复扩展。这个网格是自动优化。6A methodology for adaptive finite element analysis:Towards an integrated computational environmentG. H. Paulino, I. F. M. Menezes, J. B. Cavalcante Neto, L. F. MarthaAbstract This work introduces a methodology for self-adaptive numerical procedures, which relies on the variouscomponents of an integrated, object-oriented, computa-tional environment involving pre-, analysis, andpost-processing modules. A basic platform for numericalexperiments and further development is provided, whichallows implementation of new elements/error estimatorsand sensitivity analysis. A general implementation of theSuperconvergent Patch Recovery (SPR) and the recentlyproposed Recovery by Equilibrium in Patches (REP) ispresented. Both SPR and REP are compared and used forerror estimation and for guiding the adaptive remeshingprocess. Moreover, the SPR is extended for calculatingsensitivity quantities of first and higher orders. The mesh(re-)generation process is accomplished by means ofmodern methods combining quadtree and Delaunay tri-angulation techniques. Surface mesh generation in arbi-trary domains is performed automatically (i.e. with nouser intervention) during the self-adaptive analysis usingeither quadrilateral or triangular elements. These ideas areimplemented in the Finite Element System Technology inAdaptivity (FESTA) software. The effectiveness and ver-satility of FESTA are demonstrated by representative nu-merical examples illustrating the interconnections amongfinite element analysis, recovery procedures, error esti-mation/adaptivity and automatic mesh generation.Key words finite element analysis, error estimation, ad-aptivity, h-refinement, sensitivity, superconvergent patchrecovery (SPR), recovery by equilibrium in patches (REP),object oriented programming (OOP), interactive computergraphics.1IntroductionThis work presents an integrated (object-oriented) com-putational environment for self-adaptive analyses of ge-neric two-dimensional (2D) problems. This environmentincludes analysis procedures to insure a given level of ac-curacy according to certain criteria, and also the proce-dures to generate and modify the finite elementdiscretization. This computational system, called FESTA(Finite Element System Technology in Adaptivity), involvesfive main components (see shaded boxes in Figure 1):? A graphical preprocessor, for defining the geometry ofthe problem, the initial finite element mesh (togetherwith boundary conditions), and the main parametersused in a self-adaptive analysis. Here the geometricalmodel is dissociated from the finite element model.? A finite element module for solving the current boun-dary value problem. The code has been developed sothat it is highly modular, expandable, and user-friendly.Thus, it can be properly maintained and continued.Moreover, other users/developers should be able tomodify the basic programming system to fit their spe-cific applications.? An error estimation and sensitivity module. Discreti-zation errors are estimated according to available re-covery procedures, e.g. Zienkiewicz and Zhu (ZZ),superconvergent patch recovery (SPR) and recovery byequilibrium in patches (REP). Sensitivities of variousorders (1st., 2nd. or higher) are calculated by means of aprocedure analogous to the SPR. The user chooses thedesired error estimator and sensitivity order.? A mesh (re-)generation (rather than mesh enrichment)procedure, based on the combination of quadtree andDelaunay triangulation techniques. According to themagnitude of the error, calculated in the previousmodule, a new finite element mesh is automaticallygenerated (i.e. with no user intervention), using eithertriangular or quadrilateral elements (h-refinement).Computational Mechanics 23 (1999) 361388 ? Springer-Verlag 1999361G.H. PaulinoDepartment of Civil and Environmental Engineering,University of Illinois at urbana-champaign2209 Newmark Laboratory,205 North Mathews Avenue,Urbana, IL 61801-2352, U.S.A.I.F.M. Menezes, J.B. Cavalcante Neto, L.F. MarthaTeCGraf (Computer Graphics Technology Group), PUC-Rio,Rio de Janeiro, R.J., 22453-900, BrazilJ.B. Cavalcante Neto, L.F. MarthaDepartment of Civil Engineering, PUC-Rio,Rio de Janeiro, R.J., 22453-900, BrazilCorrespondence to: G.H. PaulinoG.H. Paulino acknowledges the support from the United StatesNational Science Foundation (NSF) under Grant No. CMS-9713798. I.F.M. Menezes acknowledges the financial supportprovided by the FAPERJ, which is a Brazilian agency for researchand development in the state of Rio de Janeiro. G.H. Paulino andI.F.M. Menezes also acknowledge the Department of Civil andEnvironmental Engineering at UC-Davis for hospitality while partof this work was performed. J.B. Cavalcante Neto and L.F. Marthaacknowledge the financial support provided by the Brazilianagency CNPq. The authors also thank an anonymous reviewer forproviding relevant suggestions to this work.? Finally, a postprocessor module, where all the analysisresults (e.g. deformed shape, sensitivity and stresscontours) can be visualized.Essentially, FESTA is a computational laboratory whichoffers a basic platform for numerical analysis and furtherdevelopment, e.g. implementation of new error estimators,elements, or material models (Cavalcante Neto et al. 1998).Object-oriented programming and integration of pre-,analysis, and post-processing modules make FESTA asoftware well-suited for both practical engineering appli-cations and further research development.The remainder of this paper is organized as follows.A motivation to the work and a brief literature revieware provided in Sect. 2. Afterwards, Sect. 3 presents sometheoretical background on self-adaptive simulations andan overview of the graphical interface used in the FESTAsoftware. Section 4 introduces the mathematical formula-tion of the SPR (using weighted least square systems), theREP, and the sensitivity method. A discussion about theautomatic mesh generation techniques used in this work isgiven in Sect. 5. Relevant information regarding the im-plementation of FESTA is presented in Sect. 6, especiallyaspects related to the SPR and REP recoveries. In order toassess the effectiveness of the proposed computationalsystem, representative numerical examples are given inSect. 7. Finally, in Sect. 8, conclusions are inferred anddirections for future research are discussed.2Motivation and related workNormal practice to solve engineering problems by meansof the Finite Element Method (FEM) or the BoundaryElement Method (BEM) involves increasing the number ofdiscretization points in the computational domain andresolving the resulting system of equations to examine therelative change in the numerical solution. In general, thisprocedure is time consuming, it depends on the experienceof the analyst, and it can be misleading if the solution hasnot entered an asymptotic range.Ideally, with a robust and reliable self-adaptive scheme,one would be able to specify an initial discrete modelwhich is sufficient to describe the geometry/topology ofthe domain and the boundary conditions (BCs), and tospecify a desired error tolerance, according to an appro-priate criterion. Then, the system would automaticallyrefine the model until the error measure falls below theprescribed tolerance. The process should be fully auto-matic and without any user intervention. This is the maingoal which motivated the development of FESTA. Thisapproach increases the overall reliability of the analysisprocedure since it does not depend on the experience, orinexperience, of the analyst.The need for developing better pre-processing tech-niques for the FEM, for performing automated analysis,and for obtaining self-adaptive solutions (which is be-coming a trend for commercial FEM software) have driventhe development of automatic mesh generation algo-rithms, i.e. algorithms which are capable of discretizing anarbitrary geometry into a consistent finite element meshwithout any user intervention. Several algorithms for 2Dgeometries have been developed (e.g. Baehmannet al. 1987; Blacker and Stephenson 1991; Zhu et al. 1991;Potyondy et al. 1995b; Borouchaki and Frey 1998), andapproaches for three-dimensional (3D) geometries haveappeared more recently (e.g. Cass et al. 1996; Escobar andMontenegro 1996; Beall et al. 1997; Lo 1998). The presentwork focus on automatic 2D mesh generation in connec-tion with adaptive solutions. Efficient techniques for gen-erating all-quadrilateral and all-triangular meshes areconsidered in detail. Although the algorithms presentedherein could be extended to mixed meshes, i.e. mesheswith both triangular and quadrilateral elements (see, forexample, Borouchaki and Frey 1998), this topic is notwithin the scope of this work.There exist a vast literature on error estimation andadaptivity, and the reader is directed to the appropriatereferences1. The volumes edited by Brebbia and Alia-badi (1993) and Babus ka et al. (1986) review adaptivetechniques for the FEM and the BEM. The book edited byLadeve ze and Oden (1998) presents a compilation of pa-pers from the workshop of New Advances in AdaptiveComputational Mechanics, held at Cachan, France, 1719September 1997, which dealt with the latest advances inadaptive methods in mechanics and their impact onsolving engineering problems. Several issues of journalshave also been dedicated to adaptivity, e.g. volume 12(1996), number 2 of Engineering with Computers, vol-ume 15 (1992), numbers 3/4 of Advances in EngineeringSoftware, and volume 36 (1991), number 1 of the Journalof Computational and Applied Mathematics. Surveys ofthe literature in FEM include articles by Noor and Ba-bus ka (1987), Oden and Demkovicz (1989), Strouboulisand Haque (1992a, b), Babus ka and Suri (1994), andAinsworth and Oden (1997). Mackerle (1993, 1994) hascompiled a long list of references on mesh generation,refinement, error analysis and adaptive techniques forFEM and BEM that were published from 1990 to 1993. The(ZZ, SPR, REP , .)Visualization(Postprocessor)Final DiscretizationMeshRegenerationGraphicalPreprocessor(Geometry,Topology, BCs)FiniteElement SolverConver-gence?FESTA ITERATIVE MESH DESIGN CYCLENYError Estimator(ZZ, SPR, REP , .)Visualization(Postprocessor)Final DiscretizationMeshRegenerationGraphicalPreprocessor(Geometry,Topology, BCs)FiniteElement SolverConver-gence?FESTA ITERATIVE MESH DESIGN CYCLENYError EstimatorFig. 1. Simplified diagram of the FESTA interactive meshing1The list of papers referred here is just a small sampling of theliterature, considering articles of particular interest to the presentwork, and is not intended to be a representative survey of theliterature in the field.362volume edited by Babus ka et al. (1983) presents adaptivetechniques for the FEM and the Finite Difference Method(FDM). Relatively recent textbooks in the FEM emphasizethe field of adaptive solution techniques. For example, thebook by Zienkiewicz and Taylor (1989) includes a Chapteron Error Estimation and Adaptivity (Chapter 14), whichis supplemented by the papers by Zienkiewicz andZhu (1992a, b, 1994). Moreover, the book by Szabo andBabus ka (1991) is primarily dedicated to this subject.The first papers on adaptive finite elements appearedin the early seventies. Since then, an explosive numberof papers on the subject have been published in thetechnical literature. Babus ka and Rheinboldt (1978)presented a pioneering paper about error estimates byevaluating the residuals of the approximate solution andusing them to obtain local, more accurate answers. Theydeveloped the mathematical basis of self-adaptive tech-niques. With the concept of a posteriori error estimates,one can develop a self-adaptive strategy for the FEMsuch that only certain elements should be refined.Zienkiewicz et al. (1982) presented a hierarchical ap-proach for self-adaptive methods. In the early 1980s,computer graphics techniques started to be used asstandard tools by mesh generation programs. Shep-hard (1986) published a paper where geometric model-ing and automatic mesh generation techniques wereused in conjunction with self-adaptive methods. Zien-kiewicz and Zhu (1987) introduced an error estimatorbased on obtaining improved values of gradients(stresses) using some available recovery processes. Easyto be implemented in any finite element code, this typeof technique, based on averaging and on the so called L2projection, has been used to recover the gradients, andreasonable estimators were achieved. In 1992, thistechnique was corrected/improved by the same authors,leading to the so called Superconvergent Patch Recovery:SPR (Zienkiewicz and Zhu 1992a, b, 1994). This methodis a stress-smoothing procedure over local patches ofelements and is based on a discrete least-squares fit of ahigher-order local polynomial stress field to the stressesat the superconvergent sampling points obtained fromthe finite element calculation. Attempts to improve fur-ther the recovery process can be found in various ref-erences, e.g. Wiberg et al. (1994), Wiberg andAbdulwahab (1993), Blacker and Belytschko (1994),Tabbara et al. (1994), and Lee et al. (1997). Essentiallythese improved techniques incorporate equilibrium andboundary conditions on the recovery process. Anexhaustive study by Babus ka et al. (1994a, b) showed,through numerous examples, the excellent performanceand superiority of the SPR over residual-type approaches.Recently, Boroomand and Zienkiewicz (1997) have pre-sented a new super-convergent method satisfying theequilibrium condition in a weak form, which does notrequire any knowledge of superconvergent points. Thenew recovery technique has been called Recovery byEquilibrium in Patches: REP. Both SPR and REP are ofparticular interest to the present work.As indicated above, the general field of adaptivity isbroad and has advanced significantly in recent years. Forinstance, Paulino et al. (1997) have proposed a new classof error estimates based on the concept of nodal sensi-tivities, which can be used in conjunction with generalpurpose computational methods such as FEM, BEM orFDM. Rannacher and Suttmeier (1997) have suggested afeedback approach for error control in the FEM. Mahomedand Kekana (1998) have presented an adaptive procedurebased on strain energy equalisation. Moreover, a summaryof recent advances in adaptive computational mechanicscan be found in the book edited by Ladeve ze andOden (1998).Quantification of the quality of a model with respect toanother one, taken as the reference, is of primary impor-tance in numerical analysis. This is the case with well es-tablished methods, such as the FEM, or emerging methods,such as the element free Galerkin: EFG (Chung and Bel-ytschko 1998), the symmetric Galerkin BEM (Paulino andGray 1999), or the boundary node method (Mukherjee andMukherjee 1997). Integration of concepts regarding errorestimation and adaptivity in the FEM, within a moderncomputational environment, is the focus of the presentwork.3Theoretical and computational aspectsWhenever a numerical method is used to solve the gov-erning differential equations of a boundary value problem,error is introduced by the discretization process whichreduces the continuous mathematical model to one havinga finite number of degrees of freedom. The discretizationerrors are defined as the difference between the actualsolution and its numerical approximation. By definition,the local errore ? / ?/?1?is a measure of the difference between the exact (/) and anapproximate solution (/). Here, / is analogous to a re-sponse quantity (e.g. displacements) in a typical numericalsolution procedure.Self-adaptive methods are numerical schemes whichautomatically adjust themselves in order to improve thesolution of the problem to a specified accuracy. The twobasic components in adaptive methods are error estima-tion and adaptive strategy. These components are discus-sed below.In general, there are two types of discretization errorestimates: a priori and a posteriori. Although a priori es-timates are accurate for the worst case in a particular classof solutions of a problem, they usually do not provideinformation about the actual error for a given model.A posteriori estimates use information obtained during thesolution process, in addition to some a priori assumptionsabout the solution. A posteriori estimates, which canprovide quantitatively accurate measures of the discreti-zation error, have been adopted here.In the context of adaptive strategies, extension methodshave been preferred over others approaches (e.g. dual orcomplementary methods) and are the focus of this work.These methods include h-, p-, and r-extensions. Thecomputer implementation is referred to as the h-, p-, andr-versions, respectively. In the h-extension, the mesh isautomatically refined when the local error indicator ex-363ceeds a preassigned tolerance. The p-extension generallyemploys a fixed mesh. If the error in an element exceeds apreassigned tolerance, the local order of the approxima-tion is increased to reduce the error. The r-extension(node-redistribution) employs a fixed number of nodesand attempts to dynamically move the grid points to areasof high error in the mesh. Any of these extensions can alsobe combined in a special strategy, for example, h-p-, r-h-,among others.3.1Error estimation and adaptive refinementAs pointed out by several authors (e.g. Zienkiewicz andTaylor 1989), the specification of local error in the mannergiven in Eq. (1) is generally not convenient and occa-sionally misleading. Thus mathematical norms are intro-duced to measure the discretization error. The exactdiscretization error in the finite element solution is oftenquantified on the basis of the energy norm for the dis-placement error, jjejj, which can be expressed, in terms ofstresses, asjjejj2?ZXrex? r?TD?1rex? r?dX?2?where rexand r are the exact and the finite element stressfields, D is the constitutive matrix, and X is the problemdefinition domain.The basic idea of error estimators is to substitute thefield rex, which is generally unknown, by the field ? r, ob-tained by means of recovery procedures (e.g. ZZ, SPR orREP). Therefore, the expression for computing the ap-proximate (estimated) relative error distribution jjejjescanbe expressed asjjejj2es?ZX? r ? r?TD?1? r ? r?dX?3?Taking into account the finite element discretization andconsidering a specific finite element i, Eq. (3) can be re-written asjjejjies?2?ZXi? r ? r?TD?1? r ? r?jJjdXi?4?where standard isoparametric elements have been as-sumed; ? r denotes the recovered stress field, jJj is the de-terminant of the Jacobian transformation matrix, and Xiisthe element domain.The energy norm for the error can be evaluated over thewhole domain or part of it. The contribution of all theelements in the mesh is given byjjejj2?Xmi?1jjejji?2?5?where m is the total number of elements, i refers to theelement with domain Xi, and mi?1Xi? X.The relative percentage error in the energy norm (gex)for the whole domain or part of it can be obtained asgex?jjejjexjjUjjex?6?where jjUjjexis the square root of twice the strain energy,and it is given byjjUjjex?ZXrTD?1rdX?1=2?7?The adaptive refinement strategy (h-extension) is dis-cussed next. The error estimator will define how the dis-cretization model will be refined or coarsened. A simplecriterion to achieve a solution error with an acceptablelevel for the whole domain can be stated asges? gmax?8?where gmaxis the maximum permissible error, and gesisgiven byges?jjejjesjjUhjj2? jjejj2es?1=2?9?where jjUhjj is the energy norm obtained from the finiteelement solution.A sound criterion for an optimal mesh consists ofrequiring that the energy norm error be equidistributedamong elements because it leads to meshes with highconvergence rates. Thus, for each element i,jjejjies 1:0?12?A more efficient procedure, which is adopted here, con-sists of designing a completely new mesh (re-generation)which satisfies the requirementni? 1:0?13?in the limit of mesh refinement. By assuming a certain rateof convergence (Zienkiewicz and Taylor 1989), the valueof the error ratio nican be used to decide the new size ofthe element. Thush ?hin1=pi?14?where hiis the initial size of the element, h is the final size,and p is the polynomial order of the approximation.This procedure does not account for singularity effects,such as those generated by the sharp corners and cracks. Atreatment of singular behavior, with varying degree ofsophistication, can be found in Zienkiewicz and Tay-lor (1989), Szabo and Babus ka (1991), Babus kaet al. (1994), and Coorevits et al. (1994).The ratio error (Eq. (11) is implemented in the nu-merical analysis module, and it is exported to the self-adaptive module of FESTA, where new element sizes arecalculated to satisfy the error criterion. This is done364through the Element class of the self-adaptive module,which by means of OOP, permits various element types tobe treated in the same framework of analysis, i.e. themethod is not restricted to a specific element. The deter-mination of the new element size is used either to refine orto coarsen the mesh where necessary (mesh re-genera-tion).3.2FESTA computational interfaceIn order to obtain an efficient and robust computationalsystem for self-adaptive simulations, a generic treatmentfor dealing with different finite elements (therefore, dif-ferent shapes and different degrees of interpolation func-tions) has to be taken into account. For this purpose, theFESTA software has been developed in C+ language,using an object-oriented framework2. It allows an easiercode maintenance and expansion with reduced probabilityof introducing errors (Lages et al. 1999).The most important task in an object-oriented programis the definition of the class structure. The capability ofcode reuse and extension depends heavily on the quality ofthis organization. With respect to the implementationof the finite element code and error estimator algorithms,several classes have been defined here, such as FEM, Node,Element, Material, AnModel, Gauss, Shape, LoadElem,Error, and Sens (Martha et al. 1996). The communicationamong these classes is performed through instances (ob-jects) of the classes. These objects contain, in their localrecords, abstract data pointers to objects of other classes.Figure 2 shows the class organization of FESTA.The FEM class represents the numerical discretizationof the model into finite elements. It contains references toobjects of several classes in the organization: to a listof objects of Node class (the nodes of the mesh), to a listof objects of Element class (the elements of the mesh), to alist of objects of Material class (the groups of distinctmaterials used), and a reference to an object of the An-Model class. The latter class is responsible for the specificfeatures of the type of analysis being performed (e.g.plane-stress, plane-strain, plate-bending, solid, shell). TheGauss class holds information about the numerical inte-gration process. An object of this class is associated to anelement integration point. The Shape class holds the geo-metric and field interpolation aspects of an element. Asexplained later in this paper, the Shape class is also re-sponsible for the recovery techniques, once the terms ofthe polynomial expansion are specified. The LoadElemclass consists of fictitious elements that transfer naturalboundary conditions from the elements to the nodes. TheError class is responsible for computing the relative errordistribution jjejj, as shown in Eqs. (2) to (14). Finally, theSens class is responsible for computing sensitivities usingan extension of the SPR technique.The development of the FESTA interface is based on aportable user interface toolkit, called IUP/LED (Levyet al. 1996). The portability achieved by using thisgraphical package allows the computational system to beimplemented in various computational environments, e.g.SUN SparcStation, IMB-RS6000, PC/Linux, and MicrosoftWindows-95.In order to illustrate the actual user-interface, Fig. 3shows different features of the FESTA software availableduring a self-adaptive simulation. All the windows displaydifferent instances of a frame under uniform compressiveloading on the top. The upper left window shows the initialfinite element mesh (obtained by means of transfinitemapping) and boundary conditions. The upper right andcenter right windows show the refined finite elementmeshes, using quads and triangles, respectively, after onestep of self-adaptive analysis. The center left window dis-plays both the original and displaced shape of the struc-ture. Both lower windows show contour plots: the left oneillustrates the distribution of the ratio of the error (withrespect to an average error) over the problem domain,while the right one shows the distribution of sensitivity ofthe stress field rxx, i.e. orxx=ox.4Patch recoveries, revisitedInitially, this section describes the SPR based on weightedleast square systems. Afterwards, the REP is described andthe connections between the two recovery procedures arepointed out. Next, an extension of the SPR technique forsensitivity analysis using the FEM is presented.4.1Super-convergent Patch Recovery (SPR)A generic field (e.g. stress) can be approximated by thepolynomial expansion? r ? Pa?15?where P contains the appropriate polynomial terms, anda is a set of unknown parameters. Note that this expan-sion is used for each component of the stress tensor. Forexample, for 2D problems and quadratic (8 noded iso-parametric) finite elements (see Fig. 4), the followingGaussShapeAnModelNodeFEMNodeMaterialElementAnModelLoadElemErrorSensFig. 2. FESTA OOP class organization2For further details about object-oriented programming, see thebook by Stroustrup (1991).365approximation is recommended (Zienkiewicz andZhu 1992a)P ? ?1; x; y; x2; xy; y2?16?a ? ?a0; a1; a2; a3; a4; a5?T?17?The unknown coefficient a can be obtained through aweighted least square fit of the polynomial expansion (15)to the values of r obtained from the finite element solutionat the sampling points, i.e. r. Small patches of elementsare used to perform local least square fits, and a weightingparameter (wi) is considered here to emphasize the in-fluence of the sampling points which are closer to thepatch assembly node (see Fig. 4). Thus,wi? 1=qpi?18?where qiis the Euclidean distance between the samplingpoint i and the patch assembly node, and p is an integer.In practical applications p is generally in the range be-tween 0 and 4. The case p ? 0 corresponds to the originalSPR (Zienkiewicz and Zhu 1992a), i.e. with uniformweighting. It is important to mention that the weightingfunction is effective for solving problems with steep gra-Fig. 3. Overview of the FESTA Interface366dients. It also alleviates ill-conditioning problems associ-ated with the standard SPR.To illustrate the discrete superconvergent recovery,consider a patch of elements containing m samplingpoints, as shown in Fig. 4. For a generic sampling point iin this patch, let ?xi;yi? be the Cartesian coordinates in theglobal axes. Thus, the weighted least square problem isreduced to the minimization of the following functionalG ?Xmi?1w2i r?xi;yi? ? ? r?xi;yi?2?19?Substituting Eq. (15) in Eq. (19), one obtainsG ?Xmi?1w2i r?xi;yi? ? P?xi;yi?a?2?20?The minimization problem is solved by setting oG=oa ? 0,which leads to the following set of linear algebraic equa-tionsAa ? b?21?whereA ?Xmi?1w2iPT?xi;yi?P?xi;yi?22?b ?Xmi?1w2iPT?xi;yi? r?xi;yi?23?Finally, the set of simultaneous equations given by (21) canbe solved for the unknown vector a.4.2Recovery by Equilibrium in Patches (REP)In the displacement-based FEM, once the equilibriumequation is solved, the discrete equilibrium of all elementsis assured. Thus, every isolated patch XPwill also be inequilibrium. Using this argument, Boroomand and Zien-kiewicz (1997) have shown thatZXPBTrdX ?ZXPBT rdX?24?where B is the strain-displacement matrix and, as before, r denotes the finite element solution for stresses. For theelastic case, r ? DB u?25?where u represents nodal displacement values and D is theelasticity matrix. Substitution of Eq. (25) in the right-hand-side of Eq. (24) leads toZXPBTrdX ?ZXPBTDB udX?26?Now, consider a smooth representation of the stress fieldover the patch XP, as given by Eq. (15). Using this repre-sentation on the left-hand-side of Eq. (26), one obtainsZXPBTPdX?a ?ZXPBTDBdX? u?27?This set of equations can be rewritten in matrix form asHa ? FP?28?whereH ?ZXPBTPdX?29?andFP?ZXPBTDBdX? u?30?The computation of the polynomial coefficients can beobtained by means of a least square scheme. Let the error dbed ? Ha ? FP?31?Thus the function to be minimized isF ? dTd ? Ha ? FP?THa ? FP?32?The minimization problem is solved by settingoF=oa ? 0, which leads toa ? HTH?1HTFP?33?where H and FPare given by Eqs. (29) and (30), re-spectively.An important factor affecting the accuracy of themethod is the number of elements and configuration of thepatch. For the sake of simplicity, the same procedureadopted for constructing the patches in the SPR is alsoadopted in the REP process.4.3Sensitivity calculationsIn this Section, a natural extension of the SPR is presented,which consists of its use for calculation of sensitivityquantities. The procedure is similar to the one proposedby Liszka and Orkisz (1980) for the FDM. The technique isa simple and effective means of directly computing de-rivatives in finite element analysis. This offers the possi-bility of using this method in conjunction with moderntechniques for sensitivity analysis and optimization incomputational mechanics.Fig. 4. Patch Recovery Notation: D Sampling points; Nodalvalues determined by recovery procedure; ? Nodal Points; ?Patch assembly node; qidistance between the sampling point iand the patch assembly node367Again, consider the patch of Fig. 4, and develop Taylorexpansions of a generic (differentiable) function f (e.g.stress, strain, displacement, or related quantities) aroundeach sampling point of the patch. For instance, let ?xo;yo?denote the coordinates of the patch assembly node Po.Thusfi? fo? hiofoox? kiofooy?h2i2o2foox2?k2i2o2fooy2? ? ? O?q3i?34?wherefo? f?xo;yo?;fi? f?xi;yi?;hi? xi? xo;ki? yi? yo;andqi?h2i? k2iq:?35?Note that the expansion (34) must be performed for everysampling point i in the patch. Similarly to the previoussubsection, the weighted least square problem is reducedto the minimization of the following functionalF ?Xmi?1w2ifi? fo? hiofoox? kiofooy? ?2#?36?The minimization problem can be solved by settingoF=oc ? 0; which leads to a linear system of the formAc ? d?37?where A is given by Eq. (22) withP ? 1;hi;ki;.?:?38?The right-hand-side of Eq. (37) isd ?Xmi?1w2iPTfi?39?and the linear system of equations is solved with respect tothe unknownsc ? fo;ofoox;ofooy;?T?40?The advantage of this formulation is that not only thefunction, but also its sensitivities can be easily computed.Moreover, the order of sensitivity to be calculated dependson the number of terms used in the Taylor expansion (seeEq. 34).5Mesh generation techniquesThe required meshing characteristics of an integratedsystem for interactive, self-adaptive finite element analysisare robustness, versatility, and computational efficiency.For 2D simulations, current meshing technology possessesthese characteristics. This section summarizes the devel-oped methodology for 2D self-adaptive mesh generationaimed on these characteristics. The meshing strategy isbased on four major components:? Geometry-based mesh generation The simulationmodel is defined by its geometry and topology. Thesimulation attributes are linked to topological entities.A simple topological description is adopted: a list ofregions, defined by a list of boundaries, which in turnare defined by a list of curves. The geometrical prop-erties are attached to the curves.? Independent boundary and domain recursive decom-position Boundary curves are discretized a prioriusing a recursive binary decomposition algorithm,which is based on the error estimation analysis. Eachregion domain is recursively enumerated using aquadtree decomposition, which is based on the boun-dary discretization and on the error estimation analysis.? Mesh generation combining quadtree and Delaunaytechniques The mesh is generated in two stages. Ini-tially, the interior of the region is meshed using tem-plates based on quadtree decomposition. In a secondstage, the areas between the interior mesh and theboundary of the region are meshed using a boundarycontraction technique based on a Delaunay triangula-tion property.? Automatic analysis attribute redistribution Becauseattributes are attached to topological and geometricalentities, simulation attributes are automatically reap-plied to the generated elements, element boundaries,and nodes of all meshes generated in the self-adaptiveprocess.5.1Mesh generation using triangular elementsThe present strategy for adaptive meshing is described bymeans of the example shown in Fig. 5. This figure showsan L shape in-plane stress, with fixed vertical displace-ments at the upper horizontal border, fixed horizontaldisplacements at the right vertical border, a distributedload along the left vertical border, and free tractionselsewhere on the boundary. An initial finite element meshfor this model is shown in this figure. This mesh is used asthe starting mesh for the self-adaptive process: the firsterror estimation analysis is performed based on this mesh.5.1.1Recursive boundary decompositionOne of the most important characteristics of the presentself-adaptive meshing strategy is that the boundary refine-ment is enforced independently of the domain refinement.In fact, the domain discretization requires an a prioriFig. 5. Initial mesh and BCs368boundary discretization. This enforces a better mesh gra-dation along the generated mesh boundary. In thisboundary refinement, each boundary curve is discretizedin its own parametric space. Therefore, the algorithm isgeneral for all classes of geometric curves.The algorithm used to refine each boundary curve is aone-dimensional version of the algorithm that is used torefine the domain, which is based on a quadtree technique.Each curve is decomposed using a binary tree technique.The idea consists of recursively subdividing the curve intosegments whose sizes are computed based on the charac-teristic sizes of finite elements adjacent to the curve. Thesesizes are dictated by the error estimation analysis of theprevious step in the adaptive mesh refinement.The recursive binary refinement algorithm is simple.Each test point is located with respect to the current tree.The cell where the point is located is recursively divided bytwo until the size of the resulting cell is less than thecharacteristic size of the element associated with the testpoint. This process is repeated for each element edgesegment along the curve. Each segment mid point is lo-cated in the tree resulting from the refinement due to theprevious point.Figure 6(a) shows an example of the binary refinementof a boundary curve. This straight line boundary curvecorresponds to one border of the model of Fig. 5. Thecurrent curve discretization is shown in Fig. 6(a), wherethe number of test points (in the middle of predicted re-finement segments) are indicated for each element edge.The resulting curve discretization and binary tree areshown in Fig. 6(b). Note that the binary tree is also in-fluenced by the characteristic size (hi) of the element.Figure 7 shows the refinement of all boundary curves ofthe example model. As mentioned previously, in the pres-ent meshing strategy, simulation attributes are attached tothe curves themselves instead of the element sides or nodes.These attributes, after the boundary refinement, are redis-tributed to the new element sides and nodes.5.1.2Recursive domain decompositionSeveral algorithms with different strategies have beenpublished in the literature for automatic meshing of 2Dregions. Among them, the algorithms based on recursivespatial enumeration using quadtree (Yerry and Shep-hard 1984; Baehmann et al. 1987) and the algorithmsbased on Delaunay triangulation (Joe 1986; Chew 1989;Lo 1989; Florani and Puppo 1992; Potyondy et al. 1995b)are the most robust ones.The quadtree technique has been successfully used for2D finite element meshes for several years (Yerry andShephard 1984; Baehmann et al. 1987). Due to propertiesof the quaternary tree, the algorithm is fast, efficient, andproduces a good transition between regions of differentdegrees of mesh refinement. One of the problems with thistechnique is that the boundary refinement is dictated bythe domain cell decomposition. This produces an irregulargeneration of boundary nodes with little positioning con-trol. The algorithm, in its original version, cannot conformwith a given boundary refinement, which is a very usefulproperty for the combination of different meshing algo-rithms and is essential for local mesh modifications.The present adaptive meshing strategy combines theadvantages of the quadtree and boundary contraction4444333322001012433322abFig. 6. a Current curve discretization; bCorresponding binary treeFig. 7. Boundary refinementFig. 8. Initial quadtree369techniques. The present technique can be considered as ahybrid algorithm between the traditional Shephardsquadtree approach (Baehmann et al. 1987) and Potyondysboundary contraction technique with a quadtree interiorpoint generation (Potyondy et al. 1995b). The domaindiscretization conforms with the a priori boundary re-finement (described previously). No additional nodes arecreated on the boundary curves.Figures 8 to 14 illustrate the current meshing algorithmcorresponding to the adopted L shape example. Figure 8shows the initial quadtree, which is defined solely on thepredefined boundary refinement. This tree is further re-fined based on the error analysis, i.e. the quadtree cell sizesare defined considering the characteristic finite elementsizes predicted by the error estimator. This is shown inFig. 9.Following the traditional quadtree meshing ap-proach (Baehmann et al. 1987), finite elements are gener-ated in the interior cells using templates. These templatesrequire the maximum difference in depth level for twoadjacent cells to be one. For the model example, Fig. 10shows the cell decomposition after refining the interiorcells to conform with this requirement. Figure 11 showsgenerated elements in the interior cells using templates fortriangles, together with the associated quadtree, andFig. 12 shows the final interior mesh.A major difference between the present algorithm andthe traditional quadtree meshing is that only interior cellsare considered for generating elements based on thequadtree decomposition. The narrow areas between theinterior cells and the domain boundary are meshed in aunique process. In the present case, a boundary contrac-tion procedure generates the mesh in the remaining area.No new interior node is generated in this process. Aproperty of the Delaunay triangulation is used for thecreation of triangular elements. Given a segment of thecurrent boundary, the selection of a boundary node for thecreation of a triangle is based on the maximum includedangle. Because the boundary is not necessarily convex,additional checks are needed to avoid triangle overlap-ping (Shaw and Pitchen 1978). Figure 13 shows the gen-erated triangular elements along the model boundary.A problem with the boundary contraction technique isthe quadratic complexity with respect of the number ofFig. 9. Quadtree considering error analysisFig. 10. Quadtree for one depth levelFig. 11. Quadtree and interior meshFig. 12. Interior meshFig. 13. Resulting mesh370mesh nodes. In the present case, the algorithm is enhancedby taking into account the existing quadtree domain de-composition. The selection of nodes for triangle or quad-rilateral creation exploits the quadtree data structure toavoid testing of nodes that are not in the vicinity of theboundary segment in consideration. This is certainly animportant efficiency gain of the present algorithm whencompared to previous boundary contraction procedures.The final step of the mesh generation is a node coor-dinate smoothing by averaging the coordinates of adjacentnodes. For the present example, the final mesh is shown inFig. 14.5.2Mesh generation using quadrilateral elementsFor the generation of quadrilateral elements, a combina-tion of the procedure devised by Potyondy et al. (1995b)and the templates of the quadtree technique is adoptedhere. An example, parallel to the one in the previoussection, is given in Figs. 15 to 26. The boundary contrac-tion algorithm generates quadrilaterals (preferentially)and triangles using pairs of boundary segments (cfFig. 17). To this effect, pairs of boundary segments anddouble interior cells (as mentioned previously) in thequadtree are considered. After the generation of quadri-laterals and triangles by the boundary contraction proce-dure, templates are used for the creation of quadrilateralfinite elements. This algorithm requires an even number ofboundary segments. This is achieved in the boundary re-finement by subdividing the smallest boundary segment ofeach curve into two if the final number of segments of theFig. 14. Final mesh (after smoothing)Fig. 15. Initial mesh and BCsFig. 16. Boundary refinementFig. 17. Initial quadtreeFig. 18. Quadtree considering error analysisFig. 19. Quadtree for one depth level371curve is odd. The subdivision of the smallest segmentguarantees the boundary gradation and provides a morerefined boundary where the error is larger.Consider the initial finite element mesh, made up oflinear quadrilateral elements (Q4), as illustrated in Fig. 15(cf Fig. 5). Figure 16 shows the boundary discretization, asa result of the recursive decomposition technique, fromthe error analysis obtained using the initial mesh.The use of pairs of segments leads to some modifica-tions in the treatment of the quadtree, with respect to therecursive domain decomposition. The generation of theinitial quadtree is performed considering these pairs ofsegments, which results in cells with double the size whencompared to the resulting sizes if the boundary segmentswere considered one by one. Figure 17 illustrates thispoint for the current example.The refinement of the quadtree due to the error anal-ysis calculated in the previous mesh also takes into ac-Fig. 20. Quadtree and interior meshFig. 21. Interior meshOriginal Polygon NodesAdded polygon NodesabOriginal Polygon NodesAdded polygon NodesabFig. 22 a, b. Transforming a collection of polygons into a col-lection of strictly 4-sided polygons. The original collection of 3and 4 sided polygons is shown in a, while the resulting collectionof 4-sided polygons is shown in bFig. 23. Mesh with quads and trianglesFig. 24. Resulting meshPi=e=1PMe1P+2e+2P3e= number of elementsadjacent to nodeiM14(M)123ieFig. 25. Quadrilateral element mesh smoothing procedure. Theposition of each internal node is adjusted to lie at the centroid ofthe polygon formed by its adjacent elementsFig. 26. Final mesh (after smoothing)372count (for the same reason) double the characteristic sizeof the element dictated by the error analysis, to be con-sistent with the generation of the initial quadtree. Fig-ure 18 shows the quadtree after consideration of the erroranalysis.Accordingly, the segment pairs are used to classify thecells as boundary, with respect to the distance to theboundary contour. The vertex cells consider only thevertices of the segment pairs, i.e. they do not consider thevertex between the segments, although this vertex belongsto the boundary. These particular vertices are consideredat the final stage of the algorithm. This classification is alsoused for mesh generation using triangles. However, in thiscase, each boundary segment is considered independently,while here, for quadrilateral mesh generation, pairs ofsegments are adopted.The treatment of the quadtree in order to ensure onlyone depth level between adjacent cells is the same one usedfor triangular elements, as explained in the previous Sec-tion. Figure 19 illustrates the quadtree with only one depthlevel.The templates for quadrilateral elements were devisedin such a way that the cells have twice the size of the cellsfor generating triangular elements. In the present case, thisis considered to refine the quadtree cells taking into ac-count the boundary discretization and the error estimationanalysis. These templates are then employed in the interiorcells, according to the observations described above. Thefinite elements generated by the templates, together withthe associated quadtree, are given in Fig. 20. The interiormesh is shown in Fig. 21.The boundary contraction algorithm also has to bemodified in order to consider pairs of boundary segments.It initially generates a mesh with quadrilateral (preferably)and triangular polygons that will be treated in order togenerate the final quadrilateral elements. The criterionfor the quadrilateral element generation (Potyondyet al. 1995b) is similar to the one adopted for the trian-gular meshing, i.e. it adopts a best angle criterion, ac-cording to the Delaunay triangulation. Each polygon of thecurrent mesh is subdivided into a set of quadrilateral el-ements by inserting a new vertex in its centroid. Figure 22illustrates this idea.Note that the new vertices obtained, when each edge ofthe polygon is subdivided, may already exist (it occurs, forexample, in a boundary edge). Thus the algorithm shouldverify if the vertex to be created already exists. An opti-mized algorithm, which uses information from the treeand avoids exhaustive searches, is adopted here. Figure 23shows the finite element mesh with both triangular andquadrilateral elements at the boundary region, and Fig. 24presents the resulting mesh with only quadrilateral ele-ments.Finally, a smoothing procedure is applied based on analgorithm proposed by Zhu et al. (1991). Each internalvertex is re-located towards the centroid of the polygon ofits vicinity. This procedure tends to generate quadrilateralsof optimum shape and is illustrated in Fig. 25. Fig. 26shows the final mesh after this smoothing procedure.6Numerical implementation aspectsThis Section presents some implementation aspects of theFESTA software, illustrating the interconnections amongdifferent modules of the whole system. Specific pseudo-codes for the SPR and REP techniques are provided, whichshow how to transform the governing formulae to matrixform using an OOP philosophy.The notation adopted for the pseudo-codes is explainednext.? The initial capital letters of each pseudo-code nameindicates the OOP class which is responsible for im-plementing the corresponding method, e.g. FEM_-StressRecoverySPR means that this method belongs tothe FEM class;? Each pseudo-code has a definition line on the top thatdefines the input and output explicitly, which are in-dicated as in or out, respectively;? Algorithmic key words such as if, begin, end, foreach,among others, are indicated in bold;? Vectors and matrices are denoted by f?g and ?, re-spectively, while components of vectors or matrices aredenoted by means of subscripts (e.g. f?g?i?, ?i;j?);? Comments are indicated in italic (after symbol /).Only the most important methods, which are necessary forunderstanding the SPR techniques, are presented here,namely FEM_StressRecoverySPR, NODE_ComputeStressSPR, ELEM_ComputeAbMatricesSPR, and ELEM_-RecoverStressSPR. These are presented in Pseudo-codes 1to 4, respectively. It is worth mentioning that, in order toimplement the REP technique, relatively small modifica-tions in the SPR pseudo-codes are needed. These modifi-cations are mostly related to the implementation ofEqs. (24) to (33), i.e. computation of the matrices H andF. Therefore, only the pseudo-code responsible for com-puting the matrices A and b (Eqs. (15) to (23) needs tobe changed in order to implement the REP technique.Thus, according to the notation adopted above, thispseudo-code is renamed ELEM_ComputeHFMatricesREP,and it is presented in Pseudo-code 5.Pseudo-code 1: FEM_StressRecoverySPRFEM_StressRecoverySPR(out: ?rec stress?; out: fcounterg)begin/ Get total number of nodes in the meshnum mesh node ? NODE_GetNumMeshNodes();/ Get number of stress componentsnum stress comp ? ANMODEL_GetNumStressComp();373/ Initialize recovered stress data structureforeach nodal point ?i?; i ? 1;.;num mesh nodebeginfcounterg?i? ? 0foreach stress component ?j?; j ? 1;.;num stress compbegin?rec stress?i;j? ? 0:0endend/ Loop over the nodes to recover nodal stresses by means of the SPR procedureforeach nodal point ?i?; i ? 1;.;num mesh nodebeginNODE_ComputeStresSPR (i; ?rec stress?; fcounterg)end/ Update nodal recovered stressesforeach nodal point ?i?; i ? 1;.;num mesh nodebeginforeach stress component ?j?; j ? 1;.;num stress compbegin?rec stress?i;j? ? ?rec stress?i;j?= fcounterg?i?;endendendPseudo-code 2: NODE_ComputeStressSPRNODE_ComputeStressSPR(in: node; out: ?rec stress?; out: fcounterg)begin/ Check to see if current node can be considered as a PATCH NODE:/ Get list of adjacent elementsnum adj elem ? ELEM_GetAdjElem(node; fadj elemg);if(num adj elem ? 2) return;/ Check to see if adjacent elements have different materialsmaterial ? ELEM_GetMaterial(fadj elemg?1?);foreach adjacent element ?i?; i ? 2;.;num adj elembegincurrent material ? ELEM_GetMaterial(fadj elemg?i?);if (material 6? current material ) return;end/ Check to see if current node is a corner nodeif (SHAPE_CornerNode(node ) ? false ) return;/ Get number of nodes per elementnum node ? SHAPE_GetNumElemNodes(fadj elemg?1?);/ Get number of polynomial termsnum poly term ? SHAPE_GetNumPolyTerms(fadj elemg?1?);/ Get number of stress componentsnum stress comp ? ANMODEL_GetNumStressComp();/ Get memory for vector of polynomial termsfpoly termg ? AllocVector(num poly term);/ Get memory for vector of recovered nodesfrec nodeg ? AllocVector(num node ? num adj elem);/ Get memory for matrices A, a, and b?A? ? AllocMatrix(num poly term; num poly term);?a? ? AllocMatrix(num poly term; num stress comp);?b? ? AllocMatrix(num poly term; num stress comp);/ Loop over the elements of the current PATCHforeach adjacent element ?i?; i ? 1;.;num adj elembegin/ Compute element contributions for matrices A and bELEM_ComputeAbMatricesSPR(node; fadj elemg?i?; ?A?; ?b?);end374/ Solve linear system of equations with multiple right hand sides: A a = b?a? ? ?A?1?b?/ Loop over the elements of the current PATCH for defining the nodes to be recoveredforeach adjacent element ?i?; i ? 1;.;num adj elembegin/ Select nodes to be recoveredELEM_GetRecoveredNodes(fadj elemg?i?; num rec node; frec nodeg);end/ Recover stress values for each selected nodeforeach recovered node ?i?; i ? 1;.;num rec nodebegin/ Get Cartesian coordinates of current recovered nodeCart coord ? NODE_GetCartCoord(frec nodeg?i?);/ Get polynomial terms for current recovered nodeSHAPE_GetPolyTerms(Cart coord; fpoly termg);/ Recover stress values for current nodeELEM_RecoverStressSPR( frec nodeg?i?; num poly term; fpoly termg, ?a?;?rec stress?; fcounterg);endendPseudo-code 3: ELEM_ComputeAbMatricesSPRELEM_ComputeAbMatricesSPR(in: node; in: elem; out: ?A?; out: ?b?)begin/ Get number of nodes per elementnum node ? SHAPE_GetNumElemNodes(elem);/ Get number of polynomial termsnum poly term ? SHAPE_GetNumPolyTerms(elem);/ Get number of stress componentsnum stress comp ? ANMODEL_GetNumStressComp();/ Get memory for vector of stressesfstressg ? AllocVector(num stress comp);/ Get memory for vector of polynomial termsfpoly termg ? AllocVector(num poly term);/ Get memory for vector of shape functionsfshapeg ? AllocVector(num node);/ Get memory for vector of nodal coordinatesfnodal coordg ? AllocVector(num node);/ Get Cartesian coordinates of element nodesSHAPE_GetCartCoord(elem; fnodal coordg);/ Get Cartesian coordinates of the patch nodepatch coord ? NODE_GetCartCoord(node);/ Get number of integration points of current elementnum int point ? GAUSS_GetNumIntPoints(elem);/ Loop over the integration pointsforeach integration point ?i?; i ? 1;.;num int pointbegin/ Get local coordinates of current integration pointlocal int point coord ? GAUSS_GetLocalCoord(i );/ Get shape functions for current integration pointSHAPE_GetShapeFunction(local int point coord; fshapeg);/ Get stress values at integration point (from previous finite element analysis)fstressg ? fFiniteElementStressDatag/ Compute Cartesian coordinates of the current integration pointforeach element node ?j?; j ? 1;.;num nodebeginCart int point coord ? Cart int point coord ? shape?j? nodal coord?j?end/ Compute Euclidean distance between integration point and patch nodedistance ? ComputeDistance(patch coord; Cart int point coord );375/ Compute weighting parameter: p is given by the user in the range from 0 to 4weight ? 1=distancep;/ Get polynomial terms for current integration pointSHAPE_GetPolyTerms(Cart int point coord; fpoly termg);/ Compute contribution for matrices A and bforeach polynomial term ?j?; j ? 1;.;num poly termbeginforeach polynomial term ?k?; k ? 1;.;num poly termbegin?A?j;k? ? ?A?j;k? weight2? poly term?j? ? fpoly termg?k?;endforeach stress component ?k?; k ? 1;.;num stress compbegin?b?j;k? ? ?b?j;k? weight2? poly term?j? ? fstressg?k?;endendendendPseudo-code 4: ELEM_RecoverStressSPRELEM_RecoverStressSPR(in: node; in: num poly term; in: fpoly termg; in: ?a?; out: ?rec stress?; out: fcounterg)begin/ Get number of stress componentsnum stress comp ? ANMODEL_GetNumStressComp();/ Update number of times that current node has been recoveredfcounterg?node? ? fcounterg?node? 1;/ Loop over the stresses componentsforeach stress component ?i?; i ? 1;.;num stress compbegin/ Loop over the polynomial termsforeach polynomial term ?j?; j ? 1;.;num poly termbegin?rec stress?node;i? ? ?rec stress?node;i? fpoly termg?j? ?a?j;i?;endendendPseudo-code 5: ELEM_ComputeHFMatricesREPELEM_ComputeHFMatricesREP(in: node; in: elem; out: ?H?; out: ?F?)begin/ Get number of nodes per elementnum node ? SHAPE_GetNumElemNodes(elem);/ Get number of degrees of freedom per nodenum dof node ? ANMODEL_GetNumDofNode(node);/ Compute number of degrees of freedom of current elementnum dof elem ? num node ? num dof node;/ Get number of polynomial termsnum poly term ? SHAPE_GetNumPolyTerms(elem);/ Get dimension of matrix B (strain-displacement matrix)dim B mat ? ANMODEL_GetDimBMatrix();/ Get total number of components to be recoverednum rec comp ? num poly term ? dim B mat;/ Get memory for element stiffness matrix?elem stiff? ? AllocMatrix(num dof elem; num dof elem);/ Get memory for element displacement vectorfelem dispg ? AllocVector(num dof elem);/ Get memory for vector of shape functions3767Computational experimentsThe effectiveness of the FESTA software is assessed bymeans of representative numerical examples. The exam-ples illustrate properties associated with the interconnec-tions among finite element analysis, error estimators/adaptive techniques and automatic mesh generationschemes using quadtree and Delaunay triangulation. Thefollowing examples are presented:(1) L-Shaped Plate(2) Plate with HolesThe first set of examples involves an L-shaped plate. Acomparative study involving both adaptive and uniformmeshes is conducted. Self-adaptive calculations are carriedout using both SPR and REP. In this study, Q8 finite ele-ments are considered. The second set of examples consistsof investigating self-adaptive solutions obtained for a platewith holes. Adaptive meshing is considered using the SPRand various element types, e.g. T3, T6, Q4, and Q8.Moreover, sensitivity calculations are also performed forthis example.7.1L-shaped plateA square plate of unit thickness with a square hole, underplane stress conditions, is considered here. Due to sym-metry, only one quarter of the domain is discretized, asshown in Fig. 5. The uniformly distributed load is q ? 1:0,the external dimensions are 100:0 ? 100:0, the dimensionsof the square hole are 50:0 ? 50:0, and the material hasE ? 105and m ? 0:3. Consistent units are used. The com-putational experiments are performed with Q8 elementsand using an error tolerance gmax? 5%. The exact errorsare computed in terms of the squared energy norm of theexact solution, jjUjj2ex? 0:311332 (see Baehmann andShephard 1989).Three adaptive steps are carried out using both the SPR(see Figs. 27 to 30) and REP (see Figs. 31 to 34). Moreover,a sequence of uniform meshes is provided in Figs. 35 to 40.Tables 1, 2 and 3 present some of the results obtainedduring the analysis, such as: number of the currentadaptive step (Step #, for the adaptive meshes) or theidentification of the mesh (Mesh #, for the uniformmeshes), total number of nodes (#Nodes), total number offshapeg ? AllocVector(num node);/ Get memory for vector of nodal coordinatesfnodal coordg ? AllocVector(num node);/ Get Cartesian coordinates of element nodesSHAPE_GetCartCoord(elem; fnodal coordg);/ Get memory for matrix B (strain-displacement matrix)?B? ? AllocMatrix(dim B mat; num dof elem);/ Get memory for matrix P (matrix with the polynomial terms)?P? ? AllocMatrix(dim B mat; num rec comp);/ Get number of integration points of current elementnum int point ? GAUSS_GetNumIntPoints(elem);/ Loop over the integration pointsforeach integration point ?i?; i ? 1;.;num int pointbegin/ Get local coordinates of current integration pointlocal int point coord ? GAUSS_GetLocalCoord(i);/ Get shape functions for current integration pointSHAPE_GetShapeFunction(local int point coord; fshapeg);/ Compute Cartesian coordinates of the current integration pointforeach element node ?j?; j ? 1;.;num nodebeginCart int point coord ? Cart int point coord ? shape?j? ? nodal coord?j?end/ Get B matrix at current integration pointELEM_GetBMatrix(local int point coord; ?B?);/ Get P matrix at current integration pointELEM_GetPMatrix(Cart int point coord; ?P?);/ Compute contribution of current element to matrix H?H? ? ?H? ? ?B?T? ?P?;end/ Compute current element stiffness matrix?elem stiff? ? ELEM_ComputeElemStiffMat(elem);/ Get element displacements (from previous finite element analysis)felem dispg ? fFiniteElementDisplacementDatag/ Compute contribution of current element to matrix F?F? ? ?F? ? ?elem stiff? ? felem dispg;end377Fig. 27. Initial meshFig. 28. 1st. Adaptive step (SPR)Fig. 29. 2nd. Adaptive step (SPR)Fig. 30. 3rd. Adaptive step (SPR)Fig. 31. Initial meshFig. 32. 1st. Adaptive step (REP)Fig. 33. 2nd. Adaptive step (REP)Fig. 34. 3rd. Adaptive step (REP)378elements (#Elem), total number of degrees of freedom(#DOF), central processing unit (CPU) time, and currentrelative error. As shown in Tables 1 and 2, the permissibleerror is achieved in the second iterative step. If the processis continued, the global error drops even more, as ex-pected. For comparable error levels, the REP leads to moreFig. 35. Uniform Mesh #1Fig. 36. Uniform Mesh #2Fig. 37. Uniform Mesh #3Fig. 38. Uniform Mesh #4Table 1. Results obtained with the adaptive meshes (using SPR technique)Step (#)#Node#Elem#DOFTime (sec)?Error (%)FEAMeshingInitial5312960.620.5812.6013441036708.601.625.282767242151635.633.482.8331194381236692.261.37?Using a SUN Sparc Station 20Fig. 39. Uniform Mesh #5Fig. 40. Uniform Mesh #6379DOFs than the SPR for this example. For the uniformmeshes, Table 3 shows that the permissible error isachieved after 3 refinement steps, i.e. for Mesh #4. Therefinement is continued a couple times leading to a finalerror of 3:49%. Although the #DOFs is comparable for Step#3 of Table 1 and Mesh #6 of Table 3, the error is smallerin the former case. Thus, selective refinement by means ofthe SPR leads to a better convergence rate than that ofuniform meshes.In Tables 4, 5 and 6, the effectivity of the error estimatesis studied. These Tables provide the Step # for adaptivemeshes, or Mesh # for uniform meshes, the energy norm ofthe finite element solution jjUhjj2, the energy norm of theapproximation of the error jjejj2esand jjejjes, the estimatederror ges(see Eq. (9), the energy norm of the exact errorjjejj2exand jjejjex(where jjejj2ex? jjUjj2ex? jjUhjj2), the exacterror gex(%), and the effectivity index h. This index isdefined byh ?jjejjesjjejjex:?41?If h 1:0, the error is overestimated, and if h 1:0, theerror is underestimated. The effectivity index for all themeshes using the SPR and REP is very close to 1:0. Asexpected, the effectivity index for uniform meshes is not asclose to 1:0 as the adaptive meshes. This can be seen fromthe last column of Tables 4 and 6.Adaptive meshes using the SPR and REP (Tables 4 and5) converge very fast to the actual value of jjUhjj2up to 3decimal digits accuracy. The uniform meshes (Table 6)show a much slower convergence. This can be verified bycomparing the 3rd. column of Tables 4, 5 and 6 with theanalytical solution.Finally,Fig. 41showsagraphoflog?jjejjes? ? log?#DOF?for the adaptive refinement using SPR and REP, and forthe uniform meshes. The adaptive meshes show betterTable 5. Parameters obtained with the adaptive meshes (using REP technique)Step#DOFkUhk2kek2eskekesges(%)kek2exkekexgex(%)hInitial960.3070740.0086630.09307516.560.0042580.06525311.691.42637115180.3096940.0013840.0372026.670.0016380.0404727.250.919203217560.3109200.0005730.0239374.290.0004120.0202973.641.179337339860.3112930.0000550.0074161.330.0000390.0062451.121.187510Table 4. Parameters obtained with the adaptive meshes (using SPR technique)Step#DOFkUhk2kek2eskekesges(%)kek2exkekexgex(%)hInitial960.3070740.0049540.07038512.600.0042580.06525311.691.07864216700.3100380.0008650.0294115.270.0012940.0359726.450.817603215160.3110400.0002490.0157792.830.0002920.0170883.060.923396323660.3112630.0000590.0076811.380.0000690.0083071.490.924684Table 3. Results obtained with the uniform meshes (using SPR technique)Mesh (#)#Node#Elem#DOFTime (sec)?(FEA)Error (%)15312960.6212.602177483363.106.9233731087209.605.334641192124824.934.465981300192061.103.89613934322736120.463.49?Using a SUN Sparc Station 20Table 2. Results obtained with the adaptive meshes (using REP technique)Step (#)#Node#Elem#DOFTime (sec)?Error (%)FEAMeshingInitial5312960.980.7216.561270775188.873.856.672883282175656.2012.434.29320106453986235.191.33?Using a SUN Sparc Station 20380convergence than the uniform ones. For this specific ex-ample, the behavior of both SPR and REP is in agreementwith observations by Boroomand and Zienkiewicz (1997).7.2Plate with holesThe numerical model here consists of a three-hole plateunder plane-strain, for which a complete self-adaptiveanalysis is performed. The geometry and boundary con-ditions, as well as some numerical parameters (in consis-tent units) used in this example, are shown in Fig. 42. Thisproblem has been studied by Baehmann (1989) and Ca-valcante Neto (1994). However, they have used triangular0.10log (llelles)0.011001000log (#D.O.F.)Adaptive Meshes (SPR)Adaptive Meshes (REP)Uniform Meshes (SPR)0.10log (llelles)0.011001000log (#D.O.F.)Adaptive Meshes (SPR)Adaptive Meshes (REP)Uniform Meshes (SPR)Fig. 41. log?kekes? ? log?#DOF?Table 7. Self-adaptive results for plate under plane-strainElementAdaptive Step #Node#Elem#DOFTime (sec)?Error (%)FEAMeshing(Initial Mesh)811181522.7841.5228.71T319781740192588.51427.0513.562291054155782737.497.45(Initial Mesh)2821185465.155.8515.91T61935421184219.559.267.0721394642275834.324.85(Initial Mesh)81591521.9843.1027.79Q41143813352846105.11422.479.72252044964103461310.534.66(Initial Mesh)223594283.605.8816.56Q811421439280838.5019.686.1022766872549496.933.83?Using a SUN Sparc Station 20= 0.35E = 10PermissibleError = 10%Thickness = 1.0q = 1.0150100200100250q30120751003075751159555= 0.35E = 10PermissibleError = 10%Thickness = 1.0q = 1.0150100200100250q30120751003075751159555Fig. 42. Plate with holesTable 6. Parameters obtained with the uniform meshes (using SPR technique)Mesh (#)#DOFkUhk2kek2eskekesges(%)kek2exkekexgex(%)h1960.3070740.0049540.07038512.600.0042580.06525311.691.07864223360.3095800.0014880.0385756.920.0017520.0418577.500.92159237200.3102410.0008860.0297665.330.0010910.0330305.920.901173412480.3105470.0006190.0248794.460.0007850.0280185.020.887970519200.3107220.0004720.0217263.890.0006100.0246984.430.879660627360.3108360.0003790.0194683.490.0004960.0222713.990.874139381elements, while in the present work, both linear andquadratic, triangular (T3, and T6) and quadrilateral (Q4,and Q8) elements are used. As indicated in Fig. 42, thepermissible relative error is gmax? 10% for all elements.Table 7 shows that, for the T6, Q4, and Q8 meshes, thepermissible error is achieved in just one iteration ofthe adaptive procedure. If the process is continued, thenthe (global) error drops even more however, as expected,there is a substantial increase in the number of DOFs. TheTable presents, for each iterative step, some results ob-tained during the self-adaptive analysis: number of thecurrent step (Adaptive Step), total number of nodes(#Node), total number of elements (#Elem), total numberof degrees-of-freedom (#DOF), FEA and meshing time,and the current relative error. The initial finite elementmeshes and the two refined meshes corresponding to theself-adaptive steps (for each element type) are shown inFigs. 43 to 45 for the T3 mesh, Figs. 46 to 48 for the T6Fig. 43. Initial finite element mesh (T3)Fig. 44. 1st. Self-adaptive step (T3)Fig. 45. 2nd. Self-adaptive step (T3)Fig. 46. Initial finite element mesh (T6)Fig. 47. 1st. Self-adaptive step (T6)Fig. 48. 2nd. Self-adaptive step (T6)382mesh, Figs. 49 to 51 for the Q4 mesh, and Figs. 52 to 54 forthe Q8 mesh. Note that the adapted mesh with T3 ele-ments, shown in Fig. 45, displays refinement out of theregion of stress concentration. This problem might beattributed to the actual error estimator and the type ofelement employed. In general, all the elements provide arelatively good mesh gradation, however, the aspect ratioobtained with triangles (Figs. 43 to 48) is better than thatobtained with quads (Figs. 49 to 54).7.2.1Sensitivity calculationThe performance of the FESTA software is evaluated bymeans of sensitivity calculations. For this purpose, theFig. 49. Initial finite element mesh (Q4)Fig. 50. 1st. Self-adaptive step (Q4)Fig. 51. 2nd. Self-adaptive step (Q4)Fig. 52. Initial finite element mesh (Q8)Fig. 53. 1st. Self-adaptive step (Q8)Fig. 54. 2nd. Self-adaptive step (Q8)383function corresponding to the rxxstress field (see the ex-ample given in the previous section) is recovered using thesensitivity formulation proposed in Sect. 3.3. The sensi-tivity results are comparable to the ones obtained with thefinite element analysis, as illustrated by Figs. 55 and 56.The sensitivities orxx=ox and orxx=oy are given in Figs. 57and 58.8Concluding remarks and extensionsA methodology for self-adaptive numerical analysis usingthe FEM has been presented. Based on five major modules namely preprocessor, finite element code, error/sensi-tivity analysis, mesh (re-)generator, and postprocessor the self-adaptive simulations are performed in a robust,versatile, flexible, and integrated (object oriented) com-putational environment called FESTA (Finite ElementSystem Technology in Adaptivity). It provides a basicplatform for numerical analyses and further development,e.g. implementation of new elements, analysis procedures,error estimators, and sensitivity methods. The effective-ness of the FESTA software is demonstrated by represen-tative numerical examples illustrating features ofautomatic mesh generators and the interconnectionsamong finite element analysis, recovery procedures, errorestimator/adaptivity and mesh generation.A general implementation of the SPR and the recentlyproposed REP is presented. Both SPR and REP are com-pared and used for error estimation and for guiding sev-eral adaptive remeshing processes. Moreover, the SPR isalso extended for calculating sensitivities of first, second,and higher orders. The mesh regeneration (rather thanFig. 55. Stress field rxxcalculated by means of sensitivity analysisFig. 56. Nodal smoothing finite element values of rxxstress field384enrichment) procedure is accomplished by combiningquadtree and Delaunay triangulation techniques. Surfacemesh generation in arbitrary domains is performed auto-matically during the self-adaptive analysis using eitherquadrilateral or triangular elements. Currently, this tech-nology is being used for 2D fracture propagation simula-tion by Potyondy et al. (1995a, b) and Arau jo et al. (1997).A natural subject for future research is the extension ofFESTA to three-dimensional (3D) problems. The mainingredient for 3D self-adaptive analysis is an automatic 3Dmesh generator. Extension of the 2D mesh generationtechniques presented here (e.g. quadtree-based) to 3D(octree-based) is currently beinginvestigated. The 3D meshgenerator provides the essential feature for 3D adaptivecalculations. It would be worth studying error estimatorsbased on SPR, REP, or the nodal sensitivity-based errorestimator recently proposed by Paulino et al. (1997), andtheir use as drivers for a 3D adaptive mesh refinementscheme. Another area of major interest is adaptive analysisin conjunction with multigrid techniques (e.g. Paulino1997; Bank 1998).Other relevant areas for future investigation includefracture mechanics and no linear problems. The solutionof general crack problems in a linear elastic fracture me-chanics (LEFM) setting involve issues such as developmentof special crack tip singular elements (e.g. Gray andPaulino 1998), corresponding meshing criteria using arosette of elements around the crack tip (e.g. Gerstle andFig. 57. Sensitivity orxx=oxFig. 58. Sensitivity orxx=oy385Abdalla Jr. 1990), and error estimation/adaptivity consid-ering singularity of the stress/strain fields (e.g. Lo andLee 1992; Jayaswal and Grosse 1992; Coorevits et al. 1994).This procedure would be especially relevant for an auto-mated crack propagation procedure, which could benefitfrom the basic capabilities of FESTA.The environment proposed herein can also be extendedfor nonlinear problems. Although adaptive procedures fornonlinear problems are not as well understood as those forlinear problems, several techniques have been presentedby, among others, Wiberg et al. (1996), Barthold et al.(1998), and Rannacher and Suttmeier (1998) for plasticityproblems, Wriggers and Scherf (1998) for contact prob-lems in plasticity, Yang et al. (1989), Zienkiewiczet al. (1990), and Owen et al. (1998) for forming processes,and Ortiz and Quigley (1991), Batra and Ko (1992), Peric et al. (1994), and Zienkiewicz et al. (1995) for strain lo-calization problems. The book edited by Ladeve ze andOden (1998) presents error estimators for nonlinear timedependent problems and corresponding adaptive compu-tational methods. Recently, Mosalam and Paulino (1999)have used mesh enrichment procedures for nonlinearproblems involving damage evolution in continuum solidmechanics, where the error is accounted for in an incre-ment-wise manner. Mesh regeneration (rather than meshenrichment) procedures would be specially advantageousfor this type of problems since areas of mesh refinement/derefinement can be effectively treated. Such problems arecurrently under investigation by the authors.As discussed above, the range of application of FESTAis potentially broad. Further use of the techniques pre-sented in this paper could contribute towards a reliableand automated environment in computational mechanics(e.g. see Tworzydlo and Oden 1993) employing finiteelement techniques.ReferencesAinsworth M, Oden JT (1997) A posteriori error estimation infinite element analysis. Comp. Meth. Appl. Mech. Eng. 142(12):188Arau jo TDP, Cavalcante Neto JB, Carvalho MTM, BittencourtTN, Martha LF (1997) Adaptive simulation of fracture pro-cesses based on spatial enumeration techniques. Int. J. RockMechanics and Mining Sciences & Geomechanics Abstracts34(34): 551Babus ka I, Chandra J, Flaherty JE (Eds.) (1983) Adaptive Com-putational Methods for Partial Differential Equations. SIAM(Society for Industrial and Applied Mathematics), Philadel-phiaBabus ka I, Von Petersdorff T, Andersson B (1994) Numericaltreatment of vertex singularities and intensity factors formixed boundary value problems for the Laplace equation inR3. SIAM Journal on Numerical Analysis 31(5):12651288Babus ka I, Rheinboldt WC (1978) A posteriori error estimates forthe finite element method. Int. J. Num. Meth. Eng.12(10):15971615Babus ka I, Strouboulis T, Upadhyay CS, Gangaraj SK (1994a) Amodel study of the quality of a posteriori estimators for linearelliptic problems error estimation in the interior of patchwiseuniform grids of triangles. Comp. Meth. Appl. Mech. Eng.114:307378Babus ka I, Strouboulis T, Upadhyay CS, Gangaraj SK, Copps K(1994b) Validation of a posteriori error estimators by a nu-merical approach. Int. J. Num. Meth. Eng. 37(7):10731123Babus ka I, Suri M (1994) The P and H-P versions of the finiteelement method, basic principles and properties. SIAM Re-view 36(4):578632Babus ka I, Zienkiewicz OC, Gago J, Oliveira ER de A (Eds.)(1986) Accuracy Estimates and Adaptive Refinements in Fi-nite Element Computations. John Wiley & Sons, ChichesterBaehmann PL (1989) Automated finite element modeling andsimulation. Ph.D. thesis, Rensselaer Polytechnic Institute,Troy, N.Y. U.S.A.Baehmann PL, Shephard MS (1989) Adaptive multiple-level h-refinement in automated finite element analysis. Engineeringwith Computers 5:235247Baehmann PL, Wittchen SL, Shephard MS, Grice KR, Yerry MA(1987) Robust geometrically based, automatic two-dimen-sional mesh generation. Int. J. Num. Meth. Eng. 24(6):10431078Bank RE (1988) PLTMG: A Software Package for Solving EllipticPartial Differential Equations: Users Guide 8.0. SIAM (Soci-ety for Industrial and Applied Mathematics), PABarthold F-J, Schmidt M, Stein E (1998) Error indicators andmesh refinements for finite-element-computations of elasto-plastic deformations. Comput. Mech. 22(3):225238Batra RC, Ko K-I (1992) An adaptive mesh refinement techniquefor the analysis of shear bands in plane strain compression ofa thermoviscoplastic solid. Comput. Mech. 10:369379Beall MW, Dey S, Klaas O, Shephard MS (1997) Adaptive anal-ysis within a geometry-based analysis framework. In: FourthU.S. National Congress on Computational Mechanics, p. 16,San Francisco, CABlacker TD, Belytschko T (1994): Superconvergent patch recov-ery with equilibrium and conjoint interpolant enhancements.Int. J. Num. Meth. Eng. 37(3):517536Blacker TD, Stephenson MB (1991): Paving: A new approach toautomated quadrilateral mesh generation. Int. J. Num. Meth.Eng. 32(4):811-847Boroomand B, Zienkiewicz OC (1997) Recovery by equilibriumin patches (REP). Int. J. Num. Meth. Eng. 40(1):137164Borouchaki H, Frey PJ (1998) Adaptive triangular-quadrilateralmesh generation. Int. J. Num. Meth. Eng. 41(5):915934Brebbia CA, Aliabadi MH (Eds.) (1993) Adaptive Finite andBoundary Element Methods. Computational MechanicsPublications, Southampton; Elsevier Applied Science, LondonCass RJ, Benzley SE, Meyers RJ, Blacker TD (1996) Generali