- Email: [email protected]
An iterative algorithm for adaptive element splitting using transition elements
Computers & Srrucmm Vol. 37. No. 3, pp. 283-294, 1990 Printed in Great Britain.
0045-J949po $3.00 * 0.00 Pergamon Press plc
AN ITERATIVE ALGORITHM FOR ADAPTIVE ELEMENT SPLITTING USING TRANSITION ELEMENTS R. GOPALAKRISHNAN, S. VUAYAKAR and H. BUSBY~ Department of Mechanical Engineering, The Ohio State University, Columbus, OH 43210, U.S.A. (Received 27 September 1989) Abstract-A simple iterative algorithm is described to implement the ‘h-adaptive’ convergence technique to two~jmensiona~ elements with quadratic insulation using transition elements. Examples of a thick
cylinder, a thick walled sphere and a Boussinesq problem using piecewise quadratic transition elements are presented.
Adaptive finite element methods that automatically carry out mesh refinement are classified into two categories, h-adaptive and p-adaptive methods. In p-adaptive methods, the mesh is refined by adding higher order terms, usually polynomials to the finite element shape function basis [1,2]. In the h-adaptive methods, the mesh is refined either by splitting elements [3,4] or by rebuilding the mesh from scratch using an automatic mesh generator. In the p-adaptive meshes the bandwidth of the associated system of equations tends to increase drastically with each iteration. Here it is noted that solving a banded symmetric linear system of equations using Gauss elimination is proportional to the number of unknowns and to the square of the bandwidth. Also, the size of the element matrices increases with each iteration and the CPU time required for its formulation is proportional to the square of the element matrix size. Besides, the p-convergence technique has an intrinsic limitation in terms of the attainable accuracy, because the number of ~lynomial mode shapes must be hard coded in the software. Thus, when all the mode shapes have been exhausted, no further improvement in the accuracy of the solution can be achieved if required. On the other hand, with the h-convergence technique, the amount of mesh refinement is limited only by the numerical stability and the capability of the hardware to handle large size finite element analysis. Adaptive finite analyses that use the h-adaptive scheme with element splitting offer computational advantages as compared to both the p-version and the h-version, with complete remeshing of the models. Element matrices do not need to be recomputed for the elements that have no split. The prevalent approach is to identify elements with error estimates greater than the specified value and to split them and enforce compatibility with their unsplit neighbors through constraint conditions applied I’To whom correspondence should be addressed.
during the formulation of the element stiffness matrices [3,4]. The use of such constraints, however, requires information about each element’s neighbors while formulating its stiffness matrix. A simple and efficient element representation scheme is presented in this paper which eliminates the need to store the information about the neighbors of each element. The new scheme also enables us to automati~~ly split any element without the need to modify the neighboring elements first to accommodate any new nodes that need to be added due to the splitting of the element. Roberti and Melkanoff [3] describe an hversion program using linear triangular elements. in their implementation, the critical triangles are split first, and then compatibility with the neighbors of the split elements is ensured by splitting the neighbors appropriately. The conventional quadrilateral and triangular elements when used for mesh gradation, however, very often lead to refined meshes with elements which are very much distorted. To avoid this element distortion due to automatic mesh reflnement, a set of transition elements was used in the algorithm. However, not much work has been done in the h-adaptive convergence technique as applied to higher order elements like elements with quadratic interpolation. The reason is the difficulty in maintaining mesh compatibility in the automatically refined mesh without much element distortion. However, an h-adaptive convergence technique, when applied to higher order elements, is expected to give better solution convergence rates in comparison with the technique as applied to elements with linear interpolation. A simple iterative procedure is presented here for implementing the h-adaptive convergence technique to two-dimensional quadrilateral elements with quadratic inte~olation using piecewise quadratic transition elements. ELEMENT REPRESENTATION The algorithm makes use of four different types of two-dimensional elements, a conventional eight-
283
R. GOPALAKRISHNANet al.
284 Side 4 4
Side 6 14 7
4
Side 5 13
v--h--Q3
3
0
tl 0
t-v Gj 0
0
5
:
0;
00
0
I
_:
2
s
-+t
1
Side 2
1 o -_) 0 -_)
Gauss Point Linear (corner)
t] ---)
Quadratic (midside) node
node
9 Side 1
0
---_)
0
----)
0
d
5
2
10 Side 2
Gauss Point Linear (corner) node Quadratic (mid&de) node
Fig. 1. Ten-noded piecewise quadratic two-dimensional quadrilateral transition element.
Fig. 3. ‘Fourteen-noded piecewise quadratic two-dimensional quadrilateral transition element.
noded two-dimensional quadrilateral element, a lonoded piecewise quadratic two-dimensional quadrilateral transition element (Fig. 1), a 1Znoded piecewise quadratic two-dimensional quad~lateral transition element (Fig. 2) and a 14-noded piecewise quadratic two-dimensional quadrilateral transition element (Fig. 3). A brief description of these elements along with the derivation of their shape functions has been included in the Appendix. A new and efficient element representation scheme is used in the algorithm. The proposed scheme requires storage of the coordinates of the various nodes of the finite element mesh along with their appropriate global node number in an array and also a form of two-dimensional label for the four corner nodes of each element in a different array. The two-dimensional labels of the four comer nodes of each element are used to represent each element. For example, the two-dimensional labels of the comer nodes of an element with curved sides are the coordinates of the
corner nodes in cylindrical coordinates which are later converted to the Cartesian coordinates for solving the system. This could be done because there always existed a one-to-one mapping between the two-dimensional labels and the Cartesian coordinate system, As for establishing the element connectivity, an algorithm, using the two-dimensional labels of the four corner nodes, evaluated for each element the coordinates of points along each side of the element where a node would exist if it were one of the elements used. A search algorithm searches through the array of nodes to check whether a node exists at the evaluated points. Thus, the element connectivity is established. The advantage of this type of element representation is that the splitting of any element could be carried out independently of its neighboring elements. The addition of new nodes to the interelement boundary due to the splitting of an element can be done without notifying the neighboring element. A compatible automatically refined mesh is ensured by a mesh compatibility algorithm which is executed over each element after the splitting of all the necessary elements is completed. The element connectivity needs to be updated only once after establishing mesh compatibility.
Side 5 7
1
9 Side 1
5
3
10 Side 2
2 0 +
0
0
q
---_)
--) +
Gauss Point Linear (corner) node Quadratic (midside) node
Fig. 2. Twelve-noded piecewise quadratic two-dimensional quadrilateral transition element.
Quadratic
(mIdsIde)
node
Linear (corner) nods New Quadratic(midside) node8 formd o-,
element
splitting
are ----
due to the II , 12 , 13 , 14 and IS
Fig. 4. Splitting of the IO-noded transition element in the C-direction.
285
Iterative algorithm for adaptive element splitting
from the splitting of the parent element conform with one of the elements used in the algorithm. For example, when a lo-noded transition element is split along the r-direction, the following steps need to be executed (see Fig. 4). 1. A new interelement boundary is formed with nodes 8 and 10 (of the parent element) as its comer nodes. 2. Nodes 8 and 10 of the parent elements are converted from quadratic midside nodes to linear comer nodes in the new elements. 3. Five additional quadratic midside nodes (nodes 11, 12, 13, 14 and 15) are formed. 4. The new element formed at the lower half of the parent element is assigned an element number the same as that of the parent element, whereas the new element formed at the top half of the parent element is assigned an element number equal to the number of elements (before the parent element was split) plus one. 5. The two new elements formed are assigned their appropriate four new comer nodes (node numbers correspond to the parent element). For the parent eldment, corner nodes are 1, 2, 3 and 4.
Fig. 5. Splitting of the IO-noded transition element in the
~-direction. ELEMENT SPLITTING
Each element has a (r, q) local coordinate system, and can be split into two new elements, either in the r-direction or in the t]-direction. Depending on the type of element and the direction in which it is being split some additional nodes need to be added to the new elements and some quadratic midside nodes need to be converted to linear corner nodes or vice versa. This is done to ensure that the new elements formed
13
‘k’
12
15
10
10
8-
8 ‘i’
14
11
H 0 --)
Quadratic (midside) node
O+
Linear (corner) node New Quadratic (midstde) nodes formed due to the element splitting are ----
13
, 14 and 15
Fig. 6. Splitting of the 12-noded transition element in the t-direction.
12
12
10 8
-10
11
0 --)
O+
Quadratic (midside) node Linear (corner) node
Three new Quadratic (midside) nodes formed due to element splitting in the 1 dhction are -- 13 , 14 and 15 Node 9 was converted from a Quadratic (mldside) node to a Linear (comer) node .
Fig. 7. Splitting of the 12-noded transition element in the q-direction.
8 11
R. GOPALAKRISHNAN et al.
286
m Y
1
6
5
I
2
I
6
5
7
2
q + Quadratic (midside) node O-W Linear (corner) node New Quadratic (midside) nodes formed due to the element splitting are ---- 15 , 16 and 17 Fig. 8. Splitting of the 1Cnoded transition element in the l-direction. 14
4
9
13
3
14
4
9
13
12
12
rl 10
8
t 5 ‘i’
3
8
-10
11
11
IIII 1
6
5
I
2
1
6
5
7
2
q --) Quadratic (midside) node 0 --) Linear (corner) node One new Quadratic (midside) node formed due to element splitting in the ? direction is -- 15. Fig. 9. Splitting of the 1Cnoded transition element in the o-direction.
For the new element at the lower half, corner nodes are 1, 2, 8 and 10. For the new element at the top half, corner nodes are 10, 8, 3 and 4. 6. The number of elements is incremented by one.
A similar set of six steps with appropriate node numbers was implemented when the IO-noded element was split along the q-direction (Fig. 5), but this time the new element formed at the left half of the parent element is assigned an element number the same as that of the parent element, whereas the new
element formed at the right half of the parent element is assigned an element number equal to the number of elements (before splitting the parent element) plus one. Similar steps are taken while splitting the other elements used in the algorithm, as shown schematically in Figs 5-9. MESH COMPATIBILITY
Since an element’s boundary is also shared by the element’s neighbor, addition of new nodes along the interelement boundary when the element is being split
Side 1 9 .-
Additional on side
node 1.
Fig. IO. Check for additional nodes on any side of an element.
Quadrhc node of element ‘I’ Spltt along c overlnppln9 alth llnenr aode+Two new element8 -- ‘I’ and ‘k’ Additional quadratic nodes-13,14 01 nei@bourln9 element ‘j’
Fig. 11. Check for quadratic node at quarter point of a side overlapping with a linear node.
287
Iterative algorithm for adaptive element splitting
g)yqqz~/ 1
5
Split along e Quadratic node of element ‘i’ OVertappingwith ttnur aode *Tao new elements -It’ and ‘b’ of neighbowing element ‘j’
Node 6 at the mtdpoint of side 2 is Quadratic . Hence , node 9 + should not have been there .
14. If the
5 To accamodate for the extra-node element ‘I’ 41 split along the direction t forming mo smaller
elements ‘I’ and ‘t’
at
.
a side
Fig. 12. Quadratic node at midpoint of a side of the quadrilateral overlapping with a linear node.
1
s
9
21
Missing Quadratic node --) pahvs linear nodes
105
2
9
p 300.0 psi a = 1.0 in. b = 2.0 in 0
Addttion~l Qmdratte node ‘10’ndded in between linear nodes 1 and 5 .
Fig. 13. Check for quadratic nodes on either side of a linear
node.
might lead to the following interelement boundary.
situations
along the
1. A node on the common boundary may be a quadratic midside node for the element being split, whereas for the neighboring element it may be a linear corner node. This will lead to discontinuity in displacement along the interelement boundary. 2. The additional nodes added along the common boundary might lead to a situation wherein the neighboring element has more nodes along one side than are needed to make it piecewise quadratic. A systematic and iterative algorithm is developed to check for the above undesired situations arising due to element splitting and to implement corrective measures to ensure a compatible finite element mesh. Various element mismatch cases, along with their
a b r
Fig. 15. A thick cylinder subjected to uniform internal pressure.
measures considered appropriate corrective in the algorithm, are shown schematically in Figs 10-14.
NUMERICAL
EXAMPLES
In the numerical examples that follow, an adaptive finite element program used iterative element splitting
to refine the meshes. Three example cases have been chosen, a thick cylinder problem, a thick walled sphere problem and a Boussinesq (half-space) problem. These examples have been chosen to demonstrate the working of an iterative element splitting scheme rather than the pros and cons of the various error estimation techniques that may be used. Detailed discussion of the various error estimators may be found elsewhere [6, 101.
Table 1. Error norms for the thick cylinder problem (plane strain approach) A0 Aa No. of elements, no. of nodes 1 2 3 4
l/WE-02 3.953E - 03 1.316E-03 9.857E - 04
1.212E-02 2.516E - 03 9.058E - 04 6.576E - 04
3.154E-02 8602E - 03 3.196E-03 2.458E - 03
4, 21 16, 65 40, 149 76, 277
R. GOFALAKRISHNAN et al. NUMBER OF -=4
NUMBRROFELEhtENTS=40 EROFNODALPOINTS=149
Fig. 16. Initial mesh for the thick cylinder problem.
Fig. 18. Mesh 3 (thick cylinder). NUMBER OFELEMENTs= OFNODALPOlNTS=277
Fig. 19. Mesh 4 (thick cylinder).
Fig. 17. Mesh 2 (thick cylinder).
solution of the finite element problem and u is the unknown exact stress field. As the exact stresses o are usually not available, cr is approximated by d which may be obtained from d in many ways. For the following examples, 8 has been obtained by applying the nodal averaging technique to the discontinuous stress field 6. Hence, e, may be approximately computed to be
Let
“2
( ) Ileu$,
da,=
1
where Q, is the area of element ‘i’,
lie,IIL,=
lo, W%,) dQ
e,=d
-8.
e, being the error in the computed stresses Then Aa, is a measure of the error in the element 7’. Let
e,=b-a. Here B is the discontinuous stress field corresponding to the C,, continuous displacement field ti which is the
Table 2. Error norms for the sphere problem (axisyrnmetric approach)
M,“D”.h 1 2 :
(--e)), 2.596E - 02 8.282E- 03 1407E - 03 1.939E
(smoothed str!c.s) 1.029E- 02 4.331E - 03 7.624E 1.221E-03 - 04
s$L) 5.349E - 02 lS07E-02 2.443E- 03 5.174E
No. of elements, no. of nodes 4, 21 16, 65 124, 40, 441 149
289
Iterative algorithm for adaptive element splitting NUMBEROF ELEMENTS=40 NUMBER OF NODAL POINTS = 149
NUMBER OF ELEhmNTs=4
I
I
I
Fig. 20. Initial mesh for the thick walled sphere problem.
Then Aa is a measure of the error in the whole mesh. All the integrations in the above error estimators were computed by Gaussian quadrature. In the following examples the initial mesh used was made as coarse as possible so that it could adequately define the topology of the region being modeled. The adaptive procedure simply solved the finite element problem for the given mesh, computed the error estimate Aa, for each element ‘i’ and split each element 7’ for which Aa.,
I'
(maxi(Aut)) K
Fig. 21. Mesh 2 (sphere).
NUMBEROFEIEMENTS=124 NUMBER OF NODAL POINTS = 441
(1)
into four new elements, by splitting the element once along the 5 direction and once along the 1 direction. The four new elements formed were assigned a Aa value of Auj/C. The splitting was repeated till none of the elements satisfied eqn (1). For the following examples the constant K was chosen to be 1.5 and the constant C was chosen to be 2.0. The finite element problem was then solved for the new mesh, the new error estimates were computed and the procedure was repeated till the overall mesh error Au was less than the prescribed value. NUMBEROF ELEMENTS=16 R OF NODAL POINTS = 65
Fig. 22. Mesh 3 (sphere).
Fig. 23. Mesh 4 (sphere).
(a) Thick cylinder subjected to uniform internalpressure A thick cylinder subjected to a uniform internal pressure is shown in Fig. 15. Figure 16 shows the crude initial mesh, consisting of four eight-noded quadrilateral elements, that was used to start off the adaptive mesh refinement process. With adaptive mesh refinement a significant correspondence was seen between the area of high stress gradients (near the inner radius) and high mesh density. In mesh 4 a
Fig. 24. Boussinesq half-space problem.
R. GOPALAKRISHNAN et al.
290
0 + Linear (comer) node 0 + Quadratic(midside)node
Fig. 25. Initial mesh for the Boussinesq problem.
(b) Thick walled sphere subjected to uniform internal pressure
refined mesh away from the inner radius is seen as a consequence of the sudden change in mesh density in that region in the mesh of the previous iteration, i.e. mesh 3. Table 1 shows the non-dimensionalized values of the error estimators for the first four meshes. Both max(Aa,) and Aa, based on the exact stresses and the smoothed stress, were found to decrease with h-adaptive mesh refinement.
A thick walled sphere subjected to an uniform internal pressure considered as an axisymmetric problem is shown in Fig. 15. Figure 20 shows the initial coarse finite element mesh, made of the eightnoded quadrilateral elements, which was used to
Table 3. Error norms for the Boussinesq problem with window at (1, 0, l,OH2,0,2,0)
Y
(mae))i 1 2 3 4 5 6
3.573E - 02 0.1552 0.6210 2.484 9.9368 39.7475
(smoothed st r Ies) 3.573E - 02 0.1552 0.3076 0.1697 4.134E-02 4.00 E-02
No. of elements, no. of nodes
s!ZZ) 0.1935 0.7816 1.1892 0.3371 0.3565 0.3534
6, 29 12, 53 18, 77 24, 101 30, 125 36, 149
Table 4. Error norms for the Boussinesq problem [window (2.0, l,OH3,0,3,0)1
No. of elements, no. of nodes : 3 4 5 6
0.1552 3.573E - 02 0.621 2.484 9.936 39.747
7.693E 3.573E - 02 4.351E -02 1.074E - 02 l.O81E-02 1.062E - 02
0.1936 0.2921 8.182E 8.697E 8.615E 8.619E -
02 02 02 02
12, 6, 18, 24, 30, 36,
29 53 77 101 125 149
Table 5. Error norms for the Boussinesq problem with window at (5,0,5,0)<6,0,6,0)
Y
(-ax(?))i 1 2 3 4 5 6
3.573E - 02 0.1552 0.6120 2.484 9.9368 39.7475
(smoothed str!es) 1.8351E 1.0796E 2.3209E 2.3502E 2.2963E 2.2938E -
02 02 03 03 03 03
sZZs.) 7.3029E - 02 2.1229E - 02 2.441 E-02 2.224 E - 02 2.225 E - 02 2.225 E-02
No. of elements, no. of nodes 6, 29 12, 53 18, 77 24, 101 30, 125 36. 149
Iterative algorithm for adaptive element splitting NUMBEROFBLBMBNTS= 12 NUMBER OF NODAL POINTS = 53
291
NUMBEROF ELEMENTS= 24 NUMBER OF NODAL POINTS = 101
t---H I Fig. 26. Mesh 2 (Boussinesq problem).
Fig. 28. Mesh 4 (Boussinesq problem).
begin the h-adaptive
mesh refinement process. Here too significant a correspondence was seen between areas of high stress gradient and high mesh density. Table 2 shows the values of the non-dimensionalized error estimators for the different refined meshes. All the error estimators were found to decrease, indicating an improvement in the finite element solution with h-adaptive mesh refinement.
(c) Boussinesq (half -space) problem Figure 24 shows the Boussinesq problem, which is nothing but a concentrated force acting on an elastic half-space. The problem is modeled as an axisymmetric problem about the z-axis. A stress singularity occurs at the point of application of the force. Figure 25 shows the initial mesh consisting of six eightnoded quadrilateral elements used to start the refinement process. Figures 26-30 show the subsequent meshes. Tables 3-5 show the non-dimensionalized error estimators for the first six meshes. The maximum Aai was found to increase with an increase in mesh density in the neighborhood of stress singularity. This is because the order of the stress singularity is too high for the quadratic elements that have been used. The overall mesh error estimators, Aa, which were evaluated over a selected window of mesh area, however, were found to decrease with h-adaptive mesh refinement. The closeness of the window location to the stress singularity region determined the attainable accuracy in the solution with mesh refinement. Better accuracy in the solution was h’UMBBROF ELEMENTS= 18 NUMBER OF NODAL POINTS= 77
I
I
I
I
I
Fig. 27. Mesh 3 (Boussinesq problem).
obtained for windows away from the region of stress singularity. CONCLUSIONS
This paper has shown a new and efficient technique for implementing the h-adaptive mesh refinement process for higher order elements using transition elements. An advantage of using transition elements is that it keeps the element shapes from deteriorating with automatic mesh refinement. The simple element representation technique proposed in this paper brought about a drastic reduction in the amount of data stored and the number of steps involved in element splitting and maintaining element compatibility and element connectivity. The systematic and iterative algorithm developed eliminates the complicated bookkeeping involved in automatic element splitting and maintenance of mesh compatibility. It is hoped to extend the algorithm to three-dimensional analysis. NUMBEROF ELEMENTS= 30 NUMBER OF NODAL POINTS = 125
t-t-l I
I
I
I
I
I
u--i t-t-i
I
I
I I
Fig. 29. Mesh 5 (Boussinesq problem). NUMBEROF ELEMENTS= 36 NUMBER OF NODAL POINTS = 149
1
I
Fig. 30. Mesh 6 (Boussinesq problem).
-
I I
R.
292
GOPALAKRISHNAN et al.
REFERENCES
1. 2. 3. 4.
5. 6.
I. 8. 9. 10.
I.
Babuska, B. A. Szabo and I. N. Katz, The p-version finite element method. SIAM J. Numer. Anal. 18, 515-545 (1981). B. A. Szabo and A. K. Mehta. o-Convergent finite element approximations in fracture mechanics. Int. J. Numer. Engng 12, 551-560 (1978). P. Roberti and M. A. Melkanoff, Self-adaptive stress analysis based on stress convergence. Int. J. Numer. Meth. Engng 24, 1973-1992 (1987). K. K. Ang and S. Valliappan, Mesh grading technique using modified isoparametric shape functions and its application to wave propagation problems. Int. J. Numer. Meth. Engng 23, 331-348 (1986). S. M. Vijayakar, A recursive algorithm for adaptive element splitting using transition elements (personal communication). 0. C. Zienkiewicz and J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. Numer. Meth. Engng 24, 337-351 (1987). R. D. Cook, Concepts and Applications of Finite Element Analysis, 2nd Edn. John Wiley, New York (1974). K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). 0. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, London (1977). 0. C. Zienkiewicz et al., Adaptive Analysis Refinement
AL JJfL 1
f,(s)
0
-1
+l
0
8
1
f&8)
0
-1
0
+ls
4
4
rn) -1
0
+1s
and Shape Optimization-Some New Possibilities. The Optimum Shape-automated Structural Design, pp. 3-27.
Plenum Press, New York (1986). Theory of Elasticity, 11. S. Timoshenko, McGraw-Hill, New York (1934).
APPENDM:
THE PIECEWISE QUADRATIC ELEMENTS
1st Ekln.
TRANSmON
Figure Al shows the one-dimensional ‘hat’ functions f(s) and g(s) which are defined as follows:
dY+ 1
g,(s)
0
-1
0
+ls
Fig. Al. One-dimensional ‘hat’ functions. g,(s)=
g2(s) =
{
{
(1 - (2s + I)?) for s < 0 0 for s 20 for s g0 Fl - (2s - l)*) for s > 0
g3(s)= (1 - s2).
Side 4 7 m I+
4
3
P
P
Figure A2 shows the IO-noded piecewise quadratic twodimensional quadrilateral element. Its shape functions in terms of the one-dimensional ‘hat’ functionsf,(s) and gi(s) are given as: N,,=g,(O.fi(n) N,=g,(U.fi(tl) Ns=fi(S).g,(rl)
1
9 Side 1
. + 0 +
5
10
2
Side 2
Gauss Point
Linear(comer)node 0 -D Quadratic(midside)node Fig. A2. Ten-noded piecewise quadratic transition element.
293
Iterative algorithm for adaptive element splitting
The integration elements was: The integration ruie for the IO-noded transition element is:
rule used for the above transition
JJ 0+I +JJ-I0fG,vldtda tl) dt dtl +JJ+’0At, flfl.fG +JJ dt dv. 0 0 1) $1
-I
+I
-I
0
0
At, vl) dt drt = _t _,fK rt)dt drl JJ 0
Each subintegral in the above expression was computed using Gaussian quadrature. Figure A3 shows the It-noded piecewise quadratic twodimensional quadrilateral transition element used in the algorithm. Its shape functions were also expressed in terms of the ‘hat’ functions as follows:
Side 5
Side 1
Side 2
0 -+ Liw&u (comer) node n + Quadratic(midside) node
piecewise element.
quadratic
In the above expression all the subintegrals shown were evaluated using the conventional Gaussian quadrature method. The 14-noded piecewise quadratic transition element, shown in Fig. A4, was also used in the algorithm. Its shape function were also expressed in terms of the one-dimensional ‘hat’ functions as:
Side 6
Side 5
Side 1
Side 2
l + Gauss Point 0 -+ Linear (comer) node q -+ Quadratic (midside) node
* -+ Gauss Point
Fig. A3. Twelve-noded
-1
transition
Fig. A4, Fourteen-noded . piecewise quadratic element.
transition
R. GOPALAKRISHNAN et al.
294
For the above 14-noded transition element the integration rule used for integration over the element’s area was:
4 =hW *h(v)
-(s)_(g -(%)-(!$_
4 =h(t) *Ah)
.(qy-(!gq
-(s)_(q) The different subintegrals in the above expression were evaluated in the computer program using Gaussian quadrature.
0045-J949po $3.00 * 0.00 Pergamon Press plc
AN ITERATIVE ALGORITHM FOR ADAPTIVE ELEMENT SPLITTING USING TRANSITION ELEMENTS R. GOPALAKRISHNAN, S. VUAYAKAR and H. BUSBY~ Department of Mechanical Engineering, The Ohio State University, Columbus, OH 43210, U.S.A. (Received 27 September 1989) Abstract-A simple iterative algorithm is described to implement the ‘h-adaptive’ convergence technique to two~jmensiona~ elements with quadratic insulation using transition elements. Examples of a thick
cylinder, a thick walled sphere and a Boussinesq problem using piecewise quadratic transition elements are presented.
Adaptive finite element methods that automatically carry out mesh refinement are classified into two categories, h-adaptive and p-adaptive methods. In p-adaptive methods, the mesh is refined by adding higher order terms, usually polynomials to the finite element shape function basis [1,2]. In the h-adaptive methods, the mesh is refined either by splitting elements [3,4] or by rebuilding the mesh from scratch using an automatic mesh generator. In the p-adaptive meshes the bandwidth of the associated system of equations tends to increase drastically with each iteration. Here it is noted that solving a banded symmetric linear system of equations using Gauss elimination is proportional to the number of unknowns and to the square of the bandwidth. Also, the size of the element matrices increases with each iteration and the CPU time required for its formulation is proportional to the square of the element matrix size. Besides, the p-convergence technique has an intrinsic limitation in terms of the attainable accuracy, because the number of ~lynomial mode shapes must be hard coded in the software. Thus, when all the mode shapes have been exhausted, no further improvement in the accuracy of the solution can be achieved if required. On the other hand, with the h-convergence technique, the amount of mesh refinement is limited only by the numerical stability and the capability of the hardware to handle large size finite element analysis. Adaptive finite analyses that use the h-adaptive scheme with element splitting offer computational advantages as compared to both the p-version and the h-version, with complete remeshing of the models. Element matrices do not need to be recomputed for the elements that have no split. The prevalent approach is to identify elements with error estimates greater than the specified value and to split them and enforce compatibility with their unsplit neighbors through constraint conditions applied I’To whom correspondence should be addressed.
during the formulation of the element stiffness matrices [3,4]. The use of such constraints, however, requires information about each element’s neighbors while formulating its stiffness matrix. A simple and efficient element representation scheme is presented in this paper which eliminates the need to store the information about the neighbors of each element. The new scheme also enables us to automati~~ly split any element without the need to modify the neighboring elements first to accommodate any new nodes that need to be added due to the splitting of the element. Roberti and Melkanoff [3] describe an hversion program using linear triangular elements. in their implementation, the critical triangles are split first, and then compatibility with the neighbors of the split elements is ensured by splitting the neighbors appropriately. The conventional quadrilateral and triangular elements when used for mesh gradation, however, very often lead to refined meshes with elements which are very much distorted. To avoid this element distortion due to automatic mesh reflnement, a set of transition elements was used in the algorithm. However, not much work has been done in the h-adaptive convergence technique as applied to higher order elements like elements with quadratic interpolation. The reason is the difficulty in maintaining mesh compatibility in the automatically refined mesh without much element distortion. However, an h-adaptive convergence technique, when applied to higher order elements, is expected to give better solution convergence rates in comparison with the technique as applied to elements with linear interpolation. A simple iterative procedure is presented here for implementing the h-adaptive convergence technique to two-dimensional quadrilateral elements with quadratic inte~olation using piecewise quadratic transition elements. ELEMENT REPRESENTATION The algorithm makes use of four different types of two-dimensional elements, a conventional eight-
283
R. GOPALAKRISHNANet al.
284 Side 4 4
Side 6 14 7
4
Side 5 13
v--h--Q3
3
0
tl 0
t-v Gj 0
0
5
:
0;
00
0
I
_:
2
s
-+t
1
Side 2
1 o -_) 0 -_)
Gauss Point Linear (corner)
t] ---)
Quadratic (midside) node
node
9 Side 1
0
---_)
0
----)
0
d
5
2
10 Side 2
Gauss Point Linear (corner) node Quadratic (mid&de) node
Fig. 1. Ten-noded piecewise quadratic two-dimensional quadrilateral transition element.
Fig. 3. ‘Fourteen-noded piecewise quadratic two-dimensional quadrilateral transition element.
noded two-dimensional quadrilateral element, a lonoded piecewise quadratic two-dimensional quadrilateral transition element (Fig. 1), a 1Znoded piecewise quadratic two-dimensional quad~lateral transition element (Fig. 2) and a 14-noded piecewise quadratic two-dimensional quadrilateral transition element (Fig. 3). A brief description of these elements along with the derivation of their shape functions has been included in the Appendix. A new and efficient element representation scheme is used in the algorithm. The proposed scheme requires storage of the coordinates of the various nodes of the finite element mesh along with their appropriate global node number in an array and also a form of two-dimensional label for the four corner nodes of each element in a different array. The two-dimensional labels of the four comer nodes of each element are used to represent each element. For example, the two-dimensional labels of the comer nodes of an element with curved sides are the coordinates of the
corner nodes in cylindrical coordinates which are later converted to the Cartesian coordinates for solving the system. This could be done because there always existed a one-to-one mapping between the two-dimensional labels and the Cartesian coordinate system, As for establishing the element connectivity, an algorithm, using the two-dimensional labels of the four corner nodes, evaluated for each element the coordinates of points along each side of the element where a node would exist if it were one of the elements used. A search algorithm searches through the array of nodes to check whether a node exists at the evaluated points. Thus, the element connectivity is established. The advantage of this type of element representation is that the splitting of any element could be carried out independently of its neighboring elements. The addition of new nodes to the interelement boundary due to the splitting of an element can be done without notifying the neighboring element. A compatible automatically refined mesh is ensured by a mesh compatibility algorithm which is executed over each element after the splitting of all the necessary elements is completed. The element connectivity needs to be updated only once after establishing mesh compatibility.
Side 5 7
1
9 Side 1
5
3
10 Side 2
2 0 +
0
0
q
---_)
--) +
Gauss Point Linear (corner) node Quadratic (midside) node
Fig. 2. Twelve-noded piecewise quadratic two-dimensional quadrilateral transition element.
Quadratic
(mIdsIde)
node
Linear (corner) nods New Quadratic(midside) node8 formd o-,
element
splitting
are ----
due to the II , 12 , 13 , 14 and IS
Fig. 4. Splitting of the IO-noded transition element in the C-direction.
285
Iterative algorithm for adaptive element splitting
from the splitting of the parent element conform with one of the elements used in the algorithm. For example, when a lo-noded transition element is split along the r-direction, the following steps need to be executed (see Fig. 4). 1. A new interelement boundary is formed with nodes 8 and 10 (of the parent element) as its comer nodes. 2. Nodes 8 and 10 of the parent elements are converted from quadratic midside nodes to linear comer nodes in the new elements. 3. Five additional quadratic midside nodes (nodes 11, 12, 13, 14 and 15) are formed. 4. The new element formed at the lower half of the parent element is assigned an element number the same as that of the parent element, whereas the new element formed at the top half of the parent element is assigned an element number equal to the number of elements (before the parent element was split) plus one. 5. The two new elements formed are assigned their appropriate four new comer nodes (node numbers correspond to the parent element). For the parent eldment, corner nodes are 1, 2, 3 and 4.
Fig. 5. Splitting of the IO-noded transition element in the
~-direction. ELEMENT SPLITTING
Each element has a (r, q) local coordinate system, and can be split into two new elements, either in the r-direction or in the t]-direction. Depending on the type of element and the direction in which it is being split some additional nodes need to be added to the new elements and some quadratic midside nodes need to be converted to linear corner nodes or vice versa. This is done to ensure that the new elements formed
13
‘k’
12
15
10
10
8-
8 ‘i’
14
11
H 0 --)
Quadratic (midside) node
O+
Linear (corner) node New Quadratic (midstde) nodes formed due to the element splitting are ----
13
, 14 and 15
Fig. 6. Splitting of the 12-noded transition element in the t-direction.
12
12
10 8
-10
11
0 --)
O+
Quadratic (midside) node Linear (corner) node
Three new Quadratic (midside) nodes formed due to element splitting in the 1 dhction are -- 13 , 14 and 15 Node 9 was converted from a Quadratic (mldside) node to a Linear (comer) node .
Fig. 7. Splitting of the 12-noded transition element in the q-direction.
8 11
R. GOPALAKRISHNAN et al.
286
m Y
1
6
5
I
2
I
6
5
7
2
q + Quadratic (midside) node O-W Linear (corner) node New Quadratic (midside) nodes formed due to the element splitting are ---- 15 , 16 and 17 Fig. 8. Splitting of the 1Cnoded transition element in the l-direction. 14
4
9
13
3
14
4
9
13
12
12
rl 10
8
t 5 ‘i’
3
8
-10
11
11
IIII 1
6
5
I
2
1
6
5
7
2
q --) Quadratic (midside) node 0 --) Linear (corner) node One new Quadratic (midside) node formed due to element splitting in the ? direction is -- 15. Fig. 9. Splitting of the 1Cnoded transition element in the o-direction.
For the new element at the lower half, corner nodes are 1, 2, 8 and 10. For the new element at the top half, corner nodes are 10, 8, 3 and 4. 6. The number of elements is incremented by one.
A similar set of six steps with appropriate node numbers was implemented when the IO-noded element was split along the q-direction (Fig. 5), but this time the new element formed at the left half of the parent element is assigned an element number the same as that of the parent element, whereas the new
element formed at the right half of the parent element is assigned an element number equal to the number of elements (before splitting the parent element) plus one. Similar steps are taken while splitting the other elements used in the algorithm, as shown schematically in Figs 5-9. MESH COMPATIBILITY
Since an element’s boundary is also shared by the element’s neighbor, addition of new nodes along the interelement boundary when the element is being split
Side 1 9 .-
Additional on side
node 1.
Fig. IO. Check for additional nodes on any side of an element.
Quadrhc node of element ‘I’ Spltt along c overlnppln9 alth llnenr aode+Two new element8 -- ‘I’ and ‘k’ Additional quadratic nodes-13,14 01 nei@bourln9 element ‘j’
Fig. 11. Check for quadratic node at quarter point of a side overlapping with a linear node.
287
Iterative algorithm for adaptive element splitting
g)yqqz~/ 1
5
Split along e Quadratic node of element ‘i’ OVertappingwith ttnur aode *Tao new elements -It’ and ‘b’ of neighbowing element ‘j’
Node 6 at the mtdpoint of side 2 is Quadratic . Hence , node 9 + should not have been there .
14. If the
5 To accamodate for the extra-node element ‘I’ 41 split along the direction t forming mo smaller
elements ‘I’ and ‘t’
at
.
a side
Fig. 12. Quadratic node at midpoint of a side of the quadrilateral overlapping with a linear node.
1
s
9
21
Missing Quadratic node --) pahvs linear nodes
105
2
9
p 300.0 psi a = 1.0 in. b = 2.0 in 0
Addttion~l Qmdratte node ‘10’ndded in between linear nodes 1 and 5 .
Fig. 13. Check for quadratic nodes on either side of a linear
node.
might lead to the following interelement boundary.
situations
along the
1. A node on the common boundary may be a quadratic midside node for the element being split, whereas for the neighboring element it may be a linear corner node. This will lead to discontinuity in displacement along the interelement boundary. 2. The additional nodes added along the common boundary might lead to a situation wherein the neighboring element has more nodes along one side than are needed to make it piecewise quadratic. A systematic and iterative algorithm is developed to check for the above undesired situations arising due to element splitting and to implement corrective measures to ensure a compatible finite element mesh. Various element mismatch cases, along with their
a b r
Fig. 15. A thick cylinder subjected to uniform internal pressure.
measures considered appropriate corrective in the algorithm, are shown schematically in Figs 10-14.
NUMERICAL
EXAMPLES
In the numerical examples that follow, an adaptive finite element program used iterative element splitting
to refine the meshes. Three example cases have been chosen, a thick cylinder problem, a thick walled sphere problem and a Boussinesq (half-space) problem. These examples have been chosen to demonstrate the working of an iterative element splitting scheme rather than the pros and cons of the various error estimation techniques that may be used. Detailed discussion of the various error estimators may be found elsewhere [6, 101.
Table 1. Error norms for the thick cylinder problem (plane strain approach) A0 Aa No. of elements, no. of nodes 1 2 3 4
l/WE-02 3.953E - 03 1.316E-03 9.857E - 04
1.212E-02 2.516E - 03 9.058E - 04 6.576E - 04
3.154E-02 8602E - 03 3.196E-03 2.458E - 03
4, 21 16, 65 40, 149 76, 277
R. GOFALAKRISHNAN et al. NUMBER OF -=4
NUMBRROFELEhtENTS=40 EROFNODALPOINTS=149
Fig. 16. Initial mesh for the thick cylinder problem.
Fig. 18. Mesh 3 (thick cylinder). NUMBER OFELEMENTs= OFNODALPOlNTS=277
Fig. 19. Mesh 4 (thick cylinder).
Fig. 17. Mesh 2 (thick cylinder).
solution of the finite element problem and u is the unknown exact stress field. As the exact stresses o are usually not available, cr is approximated by d which may be obtained from d in many ways. For the following examples, 8 has been obtained by applying the nodal averaging technique to the discontinuous stress field 6. Hence, e, may be approximately computed to be
Let
“2
( ) Ileu$,
da,=
1
where Q, is the area of element ‘i’,
lie,IIL,=
lo, W%,) dQ
e,=d
-8.
e, being the error in the computed stresses Then Aa, is a measure of the error in the element 7’. Let
e,=b-a. Here B is the discontinuous stress field corresponding to the C,, continuous displacement field ti which is the
Table 2. Error norms for the sphere problem (axisyrnmetric approach)
M,“D”.h 1 2 :
(--e)), 2.596E - 02 8.282E- 03 1407E - 03 1.939E
(smoothed str!c.s) 1.029E- 02 4.331E - 03 7.624E 1.221E-03 - 04
s$L) 5.349E - 02 lS07E-02 2.443E- 03 5.174E
No. of elements, no. of nodes 4, 21 16, 65 124, 40, 441 149
289
Iterative algorithm for adaptive element splitting NUMBEROF ELEMENTS=40 NUMBER OF NODAL POINTS = 149
NUMBER OF ELEhmNTs=4
I
I
I
Fig. 20. Initial mesh for the thick walled sphere problem.
Then Aa is a measure of the error in the whole mesh. All the integrations in the above error estimators were computed by Gaussian quadrature. In the following examples the initial mesh used was made as coarse as possible so that it could adequately define the topology of the region being modeled. The adaptive procedure simply solved the finite element problem for the given mesh, computed the error estimate Aa, for each element ‘i’ and split each element 7’ for which Aa.,
I'
(maxi(Aut)) K
Fig. 21. Mesh 2 (sphere).
NUMBEROFEIEMENTS=124 NUMBER OF NODAL POINTS = 441
(1)
into four new elements, by splitting the element once along the 5 direction and once along the 1 direction. The four new elements formed were assigned a Aa value of Auj/C. The splitting was repeated till none of the elements satisfied eqn (1). For the following examples the constant K was chosen to be 1.5 and the constant C was chosen to be 2.0. The finite element problem was then solved for the new mesh, the new error estimates were computed and the procedure was repeated till the overall mesh error Au was less than the prescribed value. NUMBEROF ELEMENTS=16 R OF NODAL POINTS = 65
Fig. 22. Mesh 3 (sphere).
Fig. 23. Mesh 4 (sphere).
(a) Thick cylinder subjected to uniform internalpressure A thick cylinder subjected to a uniform internal pressure is shown in Fig. 15. Figure 16 shows the crude initial mesh, consisting of four eight-noded quadrilateral elements, that was used to start off the adaptive mesh refinement process. With adaptive mesh refinement a significant correspondence was seen between the area of high stress gradients (near the inner radius) and high mesh density. In mesh 4 a
Fig. 24. Boussinesq half-space problem.
R. GOPALAKRISHNAN et al.
290
0 + Linear (comer) node 0 + Quadratic(midside)node
Fig. 25. Initial mesh for the Boussinesq problem.
(b) Thick walled sphere subjected to uniform internal pressure
refined mesh away from the inner radius is seen as a consequence of the sudden change in mesh density in that region in the mesh of the previous iteration, i.e. mesh 3. Table 1 shows the non-dimensionalized values of the error estimators for the first four meshes. Both max(Aa,) and Aa, based on the exact stresses and the smoothed stress, were found to decrease with h-adaptive mesh refinement.
A thick walled sphere subjected to an uniform internal pressure considered as an axisymmetric problem is shown in Fig. 15. Figure 20 shows the initial coarse finite element mesh, made of the eightnoded quadrilateral elements, which was used to
Table 3. Error norms for the Boussinesq problem with window at (1, 0, l,OH2,0,2,0)
Y
(mae))i 1 2 3 4 5 6
3.573E - 02 0.1552 0.6210 2.484 9.9368 39.7475
(smoothed st r Ies) 3.573E - 02 0.1552 0.3076 0.1697 4.134E-02 4.00 E-02
No. of elements, no. of nodes
s!ZZ) 0.1935 0.7816 1.1892 0.3371 0.3565 0.3534
6, 29 12, 53 18, 77 24, 101 30, 125 36, 149
Table 4. Error norms for the Boussinesq problem [window (2.0, l,OH3,0,3,0)1
No. of elements, no. of nodes : 3 4 5 6
0.1552 3.573E - 02 0.621 2.484 9.936 39.747
7.693E 3.573E - 02 4.351E -02 1.074E - 02 l.O81E-02 1.062E - 02
0.1936 0.2921 8.182E 8.697E 8.615E 8.619E -
02 02 02 02
12, 6, 18, 24, 30, 36,
29 53 77 101 125 149
Table 5. Error norms for the Boussinesq problem with window at (5,0,5,0)<6,0,6,0)
Y
(-ax(?))i 1 2 3 4 5 6
3.573E - 02 0.1552 0.6120 2.484 9.9368 39.7475
(smoothed str!es) 1.8351E 1.0796E 2.3209E 2.3502E 2.2963E 2.2938E -
02 02 03 03 03 03
sZZs.) 7.3029E - 02 2.1229E - 02 2.441 E-02 2.224 E - 02 2.225 E - 02 2.225 E-02
No. of elements, no. of nodes 6, 29 12, 53 18, 77 24, 101 30, 125 36. 149
Iterative algorithm for adaptive element splitting NUMBEROFBLBMBNTS= 12 NUMBER OF NODAL POINTS = 53
291
NUMBEROF ELEMENTS= 24 NUMBER OF NODAL POINTS = 101
t---H I Fig. 26. Mesh 2 (Boussinesq problem).
Fig. 28. Mesh 4 (Boussinesq problem).
begin the h-adaptive
mesh refinement process. Here too significant a correspondence was seen between areas of high stress gradient and high mesh density. Table 2 shows the values of the non-dimensionalized error estimators for the different refined meshes. All the error estimators were found to decrease, indicating an improvement in the finite element solution with h-adaptive mesh refinement.
(c) Boussinesq (half -space) problem Figure 24 shows the Boussinesq problem, which is nothing but a concentrated force acting on an elastic half-space. The problem is modeled as an axisymmetric problem about the z-axis. A stress singularity occurs at the point of application of the force. Figure 25 shows the initial mesh consisting of six eightnoded quadrilateral elements used to start the refinement process. Figures 26-30 show the subsequent meshes. Tables 3-5 show the non-dimensionalized error estimators for the first six meshes. The maximum Aai was found to increase with an increase in mesh density in the neighborhood of stress singularity. This is because the order of the stress singularity is too high for the quadratic elements that have been used. The overall mesh error estimators, Aa, which were evaluated over a selected window of mesh area, however, were found to decrease with h-adaptive mesh refinement. The closeness of the window location to the stress singularity region determined the attainable accuracy in the solution with mesh refinement. Better accuracy in the solution was h’UMBBROF ELEMENTS= 18 NUMBER OF NODAL POINTS= 77
I
I
I
I
I
Fig. 27. Mesh 3 (Boussinesq problem).
obtained for windows away from the region of stress singularity. CONCLUSIONS
This paper has shown a new and efficient technique for implementing the h-adaptive mesh refinement process for higher order elements using transition elements. An advantage of using transition elements is that it keeps the element shapes from deteriorating with automatic mesh refinement. The simple element representation technique proposed in this paper brought about a drastic reduction in the amount of data stored and the number of steps involved in element splitting and maintaining element compatibility and element connectivity. The systematic and iterative algorithm developed eliminates the complicated bookkeeping involved in automatic element splitting and maintenance of mesh compatibility. It is hoped to extend the algorithm to three-dimensional analysis. NUMBEROF ELEMENTS= 30 NUMBER OF NODAL POINTS = 125
t-t-l I
I
I
I
I
I
u--i t-t-i
I
I
I I
Fig. 29. Mesh 5 (Boussinesq problem). NUMBEROF ELEMENTS= 36 NUMBER OF NODAL POINTS = 149
1
I
Fig. 30. Mesh 6 (Boussinesq problem).
-
I I
R.
292
GOPALAKRISHNAN et al.
REFERENCES
1. 2. 3. 4.
5. 6.
I. 8. 9. 10.
I.
Babuska, B. A. Szabo and I. N. Katz, The p-version finite element method. SIAM J. Numer. Anal. 18, 515-545 (1981). B. A. Szabo and A. K. Mehta. o-Convergent finite element approximations in fracture mechanics. Int. J. Numer. Engng 12, 551-560 (1978). P. Roberti and M. A. Melkanoff, Self-adaptive stress analysis based on stress convergence. Int. J. Numer. Meth. Engng 24, 1973-1992 (1987). K. K. Ang and S. Valliappan, Mesh grading technique using modified isoparametric shape functions and its application to wave propagation problems. Int. J. Numer. Meth. Engng 23, 331-348 (1986). S. M. Vijayakar, A recursive algorithm for adaptive element splitting using transition elements (personal communication). 0. C. Zienkiewicz and J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. Numer. Meth. Engng 24, 337-351 (1987). R. D. Cook, Concepts and Applications of Finite Element Analysis, 2nd Edn. John Wiley, New York (1974). K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). 0. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, London (1977). 0. C. Zienkiewicz et al., Adaptive Analysis Refinement
AL JJfL 1
f,(s)
0
-1
+l
0
8
1
f&8)
0
-1
0
+ls
4
4
rn) -1
0
+1s
and Shape Optimization-Some New Possibilities. The Optimum Shape-automated Structural Design, pp. 3-27.
Plenum Press, New York (1986). Theory of Elasticity, 11. S. Timoshenko, McGraw-Hill, New York (1934).
APPENDM:
THE PIECEWISE QUADRATIC ELEMENTS
1st Ekln.
TRANSmON
Figure Al shows the one-dimensional ‘hat’ functions f(s) and g(s) which are defined as follows:
dY+ 1
g,(s)
0
-1
0
+ls
Fig. Al. One-dimensional ‘hat’ functions. g,(s)=
g2(s) =
{
{
(1 - (2s + I)?) for s < 0 0 for s 20 for s g0 Fl - (2s - l)*) for s > 0
g3(s)= (1 - s2).
Side 4 7 m I+
4
3
P
P
Figure A2 shows the IO-noded piecewise quadratic twodimensional quadrilateral element. Its shape functions in terms of the one-dimensional ‘hat’ functionsf,(s) and gi(s) are given as: N,,=g,(O.fi(n) N,=g,(U.fi(tl) Ns=fi(S).g,(rl)
1
9 Side 1
. + 0 +
5
10
2
Side 2
Gauss Point
Linear(comer)node 0 -D Quadratic(midside)node Fig. A2. Ten-noded piecewise quadratic transition element.
293
Iterative algorithm for adaptive element splitting
The integration elements was: The integration ruie for the IO-noded transition element is:
rule used for the above transition
JJ 0+I +JJ-I0fG,vldtda tl) dt dtl +JJ+’0At, flfl.fG +JJ dt dv. 0 0 1) $1
-I
+I
-I
0
0
At, vl) dt drt = _t _,fK rt)dt drl JJ 0
Each subintegral in the above expression was computed using Gaussian quadrature. Figure A3 shows the It-noded piecewise quadratic twodimensional quadrilateral transition element used in the algorithm. Its shape functions were also expressed in terms of the ‘hat’ functions as follows:
Side 5
Side 1
Side 2
0 -+ Liw&u (comer) node n + Quadratic(midside) node
piecewise element.
quadratic
In the above expression all the subintegrals shown were evaluated using the conventional Gaussian quadrature method. The 14-noded piecewise quadratic transition element, shown in Fig. A4, was also used in the algorithm. Its shape function were also expressed in terms of the one-dimensional ‘hat’ functions as:
Side 6
Side 5
Side 1
Side 2
l + Gauss Point 0 -+ Linear (comer) node q -+ Quadratic (midside) node
* -+ Gauss Point
Fig. A3. Twelve-noded
-1
transition
Fig. A4, Fourteen-noded . piecewise quadratic element.
transition
R. GOPALAKRISHNAN et al.
294
For the above 14-noded transition element the integration rule used for integration over the element’s area was:
4 =hW *h(v)
-(s)_(g -(%)-(!$_
4 =h(t) *Ah)
.(qy-(!gq
-(s)_(q) The different subintegrals in the above expression were evaluated in the computer program using Gaussian quadrature.