- Email: [email protected]
Flexible mesh generator for triangular and quadrilateral areas
Flexible mesh generator for triangular and quadrilateral areas TAKEO TANIGUCHI
Department of Ovil Engineering, School of Engineering, Okayama University, Okayama, Japan
The object of this paper is to propose a new method to generate the topological properties among the input data neceesary for the user of the Finite Element Method. The method can generate the element-node relation for trlang-l-r and quadrilateral areas which are required to subdivide into different numbers of segments on their edges. This flexibility can largely remove the complexity of the mesh generation of two-dimensional area. The algorithm of the proposed method is so simple as to be used on a microcomputer. The main part of program listing written in BASIC is attached to this paper.
INTRODUCTION The term 'Finite Element Method' generally indicates the analysis, but its program consists of the following three parts: Pre-Processor, Analysis, and Post-Processor. Preprocessor works to decrease the tedious works of the FEM user at the preparation of all the input data of Analysis part. Analysis part works in accordance with the input data which are generated at Pre-processor, and the results are governed by the data generated at Pre-processor. This fact indicates that a better pre-processor can not only save the cost for the preparation of the input data but also give better results and improve the reliability of the results. Input data required at Analysis are classified into three kinds; the topological, geometrical, and physical properties. For example, the element-node relation, the co-ordinates of nodes, and Young's modulus of elements belong to the first, the second and the third property, respectively. A slight mistake of the input of the topological property may change the structure being treated as a whole, and, on the other hand, the influence of the mistake of the latter two properties may be restricted only in a part of the structure. This fact clarifies the reason of the requirement of the automatic mesh generator. We consider the mesh generation of a two-dimensional area with complicated boundary configuration which is consisted of several subareas with different physical properties. The first problem is how we can recognise the geometrical configurations of subareas, and the second problem is how they are subdivided into finite elements. Accepted Match 1987. Discussion closes September 1987.
142
Adv. Eng. Software, 1987, Vol. 9, No. 3
The Quad-Tree Method can resolve these problems at the same time, but it requires a number of operations to fit the boundary configurations.* The same difficulty is found at the application of the Mapping Function Method. 2 The Node Pattern Method has no difficulty on this aspect, but it requires long execution-time for meshing of each subarea.a The easiest way is the Blocking Method which subdivides the area into triangular and/or quadrilateral subareas and gives a mesh pattern for each subarea. The Blocking Method is easy to make a program of and to use, but it may become stiff and can not satisfy its user unless various types of mesh patterns are not prepared. The subdivision of three edges of a triangular area into the same number of segments and to subdivide two edges of a quadrilateral area into different numbers are systematic, and their programmings are easily obtained. The combination of these systematic mesh patterns may often give troubles to its user, because the user may fail into the case where all edges of a triangular or quadrailateral subarea must be subdivided into a different number of segments.
IINPUT DATAI number
of subareas
configuration mesh pattern
of subareas of subareas
number
of edges
number
of nodes
subarea-edge number
of segments
relation on each edge
x and y coordinates
[SUBDIVISION
OF AREA INTO SIMPLE
[GENERATION
[CONNECTION
Figure 1.
of nodes
SUBAREAS]
OF MESH FOR EACH SUBAREA[
OF SUBAREAS
INTO ORIGINAL
CONFIGURATION]
Flow-chart of the Blocking Method
0141-1195/871030142-08 $2.00 © 1987 Computational Mechanics Publications
N
J
M N
N
pattern for each of A(i), and the final stage joints them together to form the original configuration. The flow-chart of the Blocking Method is illustrated in Fig. 1. Input data necessary for the Blocking Method are (1) number of subareas, (2) number of nodes and edges, (3) block-node or block-edge relation, (4) mesh pattern of each subarea, (5) number of segments on each edge, (6) x and y co-ordinates of nodes, and some other information. Nodes and edges must be given so that the boundary configuration of all subareas can be sufficiently drawn by using them. Since the input-data of (3), (4), and (5) is concerned with the number of subareas, the amount of the input data is largely governed by the number of subareas. It suggests that the decreasing of the number of subareas is important for the effective use of the Blocking Method. Since each subareas is transformed into a meshed system by the introduction of mesh patterns prepared in the method, the flexibility of the Blocking Method is wholly governed by the kinds of mesh patterns prepared there. Conventional systematic meshing methods shown in Fig. 2 can be easily prepared, but the method with only this kind of mesh generator becomes very stiff. Its flexibility can be improved by preparing the mesh generators for transition region. At the same time the ones for the transition region can decrease the number of subareas required for appropriate mesh generation of the total area. We can easily find several systematic methods for this purpose, and some of them are shown in Fig. 3. Though they are useful for the increase of the flexibility, they still require some restriction of their use. The most preferable and desirable mesh genera-
L (=N+M)
N Figure 2.
Regular mesh patterns
The Blocking Method may become flexible and easy to use when this kind of mesh generator is prepared. The purpose of this study is to propose a method which can divide all edges of a triangular and quadrilateral area into different number of segments. Since the method is derived as an extension of conventional systematic mesh generation methods, it can sufficiently work on a microcomputer. Therefore, the method presented in this paper is developed on a Japanese microcomputer, PC-9801 of the Nippon Electric Company, and the program is written in BASIC.
M
N L
(=N.2M)
AkA/t/ /tA/t M
BLOCKING METHOD Let A be a two-dimensional area with complex boundary configuration. The Blocking Method consists of three stages: at the first stage, divide A into a gathering of simple triangular and/or quadrilateral subareas by connecting nodes on the boundary and also additional nodes placed in A. We note these subareas by A(/), where i indicates the ith subarea. The second stage generates an appropriate mesh
Figure 3.
N Mesh patterns for transition region
Adv. Eng. Software, 1987, Vcl. 9, No. 3
143
tor is the one which can subdivide all edges of triangular and quadrilateral areas into different number o f segments. The Blocking Method becomes really flexible by the introduction o f this kind of meshing method.
L
K
M
FLEXIBLE MESH GENERATOR In this section we propose a new flexible mesh generator applied to a triangular or quadrilateral area which is required to subdivide all edges into different number of segments. The most preferable method is a systematic one as shown in Figs. 2 and 3, because this type of meshing method requires relatively short execution-time and their programmings have simple structure. By this reason we aim to propose a new one as the extension of them. The method consists of two steps, i.e. the generation of topological properties and the modification of the geometrical properties.
(a)
N
L
M
M
(1) Generation of topological properties Firstly we consider the mesh generation for a triangular area. Let A be a triangular area whose three edges are to be subdivided into N, M, and L segments as shown in Fig. 4a, where we assume N > M > L . By the addition of an intersection line the triangle is divided into two triangles whose two edges have the same number of segments (see Fig. 4b).
K-M b)
N
L
M M
K-M N
(c) Figure 5.
M (b)
Further subdivision of these two triangles as shown in Fig. 4c clarifies that the flexible mesh generation method for a triangular area requires a mesh generator for a quadrilateral area a pair of whose opposite edges are to be divided into the same number of segments. The case of a quadrilateral area is explained in Fig. 5. We assume K > M, and N is the minimum among them. The original area is divided into a quadrilateral area, a pair of whose opposite edges have the same number o f segments, and a triangle area whose two edges have the same number of segments (see Fig. 5b). And the addition of an intersection line for the triangle clarifies that the flexible mesh generator for a quadrilateral area can be designed by using the same mesh generator as was required for the triangular
M-1
N-M
Figure 4.
I
(•) / N-M
N-I
M-I (c) Subdivision of triangular area
144 Adv. Eng. Software, 1987, Vol. 9, No. 3
Subdivision of quadrilateral area
case.
1
Now, consider the mesh generation method for a quadrilateral area a pair of whose edges are required to have the same number o f segments. The area may be considered as a gathering of M narrow subareas as shown in Fig. 6. The above consideration leads to the conclusion that if we find a method to give a mesh pattern for a narrow quadrilateral area as shown in Fig. 6, the flexible mesh generator for triangular and quadrilateral areas is completed.
L
M
M
N
11
l,
I1
e
11
I1 d
STEP 4: (We assume M = K.) Set N N = N, MM = M, and LL = L, and go to *MESH. Go to STEP 15. STEP 5: Find edges with the maximum and minimum segments. STEP 6: Divide the area into triangular subareas as shown in Fig. 4. STEP 7: Apply *MESH for these two subareas. STEP 8: Joint these subareas into the original form. Go to STEP 15. (Quadrilateral Cases) STEP 9: If a pair of opposite edges have the same number of segments, then go to STEP 10 else to STEP 11. STEP 10: (We assume N = L.) Find the minimum number among M and K. (We assume M < K . ) Set N N = M, LL = K , and M M = N . And apply *MESH. Go to STEP 15. STEP 11: Find the edge with the minimum segments. (we assume M is the minimum.)
-i-
11
d c
I1 N+2
N+3
N+4
N+5 2N+I
2N+2
-p
b
11 11
a
+ a
i1
1
2
3
4
N+5
N+~
N
11
N Figure 6
Subdivision o f quadrilateral area
N+2 N+3 N+4
I ,' A narrow quadrilateral area must be divided so that a pair of its opposite edges have a different number of segments. A procedure for this purpose is explained in Fig. 7. At first, the area is systematically subdivided and as many nodes as the number of the difference of segments on these two edges are generated on the upper edge. The generated nodes are systematically connected to nodes on the lower edge. For example, if the upper edge is required to have more two nodes than the lower edge, the third figure from the top in Fig. 7 corresponds to it. Mesh generation according to the different ordering of the additional node generation on the upper edge is explained in Fig. 8. The application of the proposed procedure to a different mesh pattern is also illustrated in Fig. 9. The repetition of this procedure to each narrow subarea in Fig. 6 finishes t h e mesh generation of a quadrilateral area, and its application to subareas in Figs. 4 and 5 completes the design of the flexible mesh generator for triangular and quadrilateral areas. The algorithm of this method is expressed as following: STEP 1 : Input anticlockwisely the number of segments of edges of the area, namely N, M, L, and K. Set N = 0 for a triangular case. STEP 2: I f N = 0 then go to STEP 3 else to STEP 9. (Triangular Cases) STEP 3: If two edges have the same number of segments, then go to STEP 4 else to STEP 5.
l" 1
', I 2
\
X,l
I 3
N+2 N+3 N+4 N+5 N+6
I
2
3
\
X,l
Figure 7.
2
3
2 +3
I
4
N
N+7 2N+3
4
N
N+2 N+3 N+4 N+5 N+6 N+7 N+8
i
N+I
4
N+I
2N+4
N+I
3N 3N+I 3N+2
N
N-1
Systematic element generation
Adv. Eng. Software, 198 7, Vol. 9, No. 3
145
routine works to generate a mesh pattern of a triangular or quadrilateral subarea which satisfies the condition that two edges have the same number of segments for a triangular case and a pair of opposite edges have the same number of segments for a quadrilateral case. Its procedures are summarised as following. (See Figs. 6 and 7) STEP 1 : Judgement of the configuration of the subarea. STEP 2: Determination of the number of segments on upper and lower edges of narrow quadrilateral subareas. STEP 3: Go to *PLUSA. STEP 4: Calculation o f nodes on surrounding edges.
+
+
Figure 8. ordering
Mesh generation by different node generation
The subroutine named *PLUSA is the mesh generator for a narrow quadrilateral subarea which is shown in Fig. 7. The program listing of *MESH including *PLUSA written in BASIC is given in Appendix. The Appendix also includes the explanation of input and output data of the method. It must be noted that this subroutine is valid for the mesh generation of the subarea which satisfies the condition of L < N*2 M, where L and N are the maximum and minimum number of segments of opposite edges of a quadrilateral area, and M is the number of narrow areas in Fig. 6. The results of the application of this method are shown in Fig. 10. Figure 10a shows a triangular area whose three edges are divided into four, five, and six segments, and Fig. 10b is the result of a quadrilateral area whose edges are divided into six, five, four, and seven segments.
~ m
(a)
Figure 9.
Mesh generation for different mesh pattern
STEP 12: Compare two edges adjacent to the edge with M segments, and divide the longer edge so that the area is divided into triangular and quadrilateral subareas as shown in Fig. 5. STEP 13: Apply *MESH for these subareas. STEP 14: Joint them so that they form the original one. STEP 15: Mesh generation is finished.
......... ' --__. '--.
i "--.
....~-~ . . . . . . . . . ~,. ,-" ._ ~ , £ S . ~.~-~-_ -3~_.- --- . ..... ,
"\
/
"
-7
. .-,
(b)
As obvious from the above procedure, the main part of the procedure is the subroutine called *MESH. This sub-
146 Adv. Eng. Software, 1987, Vol. 9, No. 3
Figure 10.
Application of flexible mesh generation
.--"
(2) Improvement o f node co-ordinates The mesh generator mentioned above aims mainly at the generation of the topological properties and it includes the node placement method which can set a series of nodes with the same distance between them as shown in Fig. 7. Therefore, the node co-ordinates generated by the method generally can not satisfy its user except in the cases as presented in Figs. 2 and 3. They must be modified. We can i'md a number of methods for the settlement of nodes in the region and in this paper we apply Laplacian mesh generation} In this method, the interior nodes are positioned so that the position vector, P(i), of node 'i' satisfies the equation 1
n
e(i) =-~n ff.=l Jr03 + e(k)J 'n' is the number of elements adjacent to/-node, and P(/) and P(k) are the position of adjacent nodes in the neighbouring element. The results of the application of this Laplacian mesh generation to the meshes in Fig. 10 are given in Fig. 11, and they show that the positions of the nodes are sufficiently improved by the method.
/\/\ i1
\ /;ix
\
~,
CONCLUDING REMARKS In this paper the author proposed a new flexible mesh generator which is an effective tool for generating the transition mesh pattern of a triangular and quadrilateral area. The proposed method requires less restriction compared to the conventional mesh generators for the transition region, and its introduction to the Blocking method can remove the difficulty of the setting of subareas and it can decrease the number of subareas in the area. Finally it leads to the decreasing of the amount of input data for the mesh generation. A more important conclusion obtained from this study is that the proposed method itself can be a new blocking method. Since the proposed method can generate a mesh pattern by giving an arbitrary number of segments for edges, most of the systematic mesh patterns can be prepared only by the proposed method. Thus, the Blocking Method may have only one subroutine for generating a mesh pattern. The proposed method aims to generate all of the topological properties for Finite Element Analysis, and the position of nodes should be modified especially for nodes in the transition region. In this present paper, the Laplacian method was introduced for this purpose and it was shown that the appropriate positioning of the nodes can be easily obtained by the method. The strategy presented in this paper suggests that the direct extension of the proposed method can make it possible to treat N-polygonal region as one subarea. The application of the proposed method after the subdivision of N-polygonal area into the appropriate triangular and/or quadrilateral subregions may complete the mesh generation for the area, and for this useful mesh generator an automatic subdivision method of N-polygon into subareas must be predpared.
REFERENCES
I.--'/.11
.I.--- ,, _.--%---____ "-....\ \..-\, .._-------... (a)
j
\
i- \ .,_1. ._.~
\
\, _ . . i
~-''~
_-
--'-----_.._ /
-
!1-'\ ~
" ,S.-
-"~A_---
/
, /
i ~--.. / "~_
'
-"---- -%.. /
/
I I
- " - - - . 2i - , . .
APPENDIX 1. PROGRAM LISTING FOR GENERATING ELEMENT-NODE RELATION
- _.._. ---..
---.,,./
Input-Data
- - - ....
, /
~
--~'~-,, \
"-,.......,-
(b)
Figure 11.
--~-.--_
------- ,,
! ~_
!-~, ~---"1-1
1 Yerry, M. A. and Shephard, M. S. A modified quadtree approach to finite element generation, 1EEE Computer Graphics and Applications 1983, 3 (1), 39 2 Haber, R. B., Shephard, M. S., Abel, J. F., Gallagher, R. H. and Greenberg, D. P. A general two dimensional finite element pre-processor utilizing discrete transfinite mapping, Int. J. Numer. Meth. in Eng. 1981, 17 (7), 1015 3 Okuda, O. Automatic triangular mesh generation of arbitrary planar domains utilizing the node-pattern for finite element analysis, J. Japan Society of Precision Engineering 1985, 51 (4), 802
Modification of node co-ordinates
--,--~_ --S.~
/-- .--"
N,t
Numbers o f
LL
(For Triangular Case NN~0)
Elements
on Three Edges
Output-Da~a NELM
Number of Elements
NODE
Number of Nodes
MTJ
Element-Node Relation
IEDBL
Nodes on Surrounding Edges
EXAMPLES (I) Trlangular C a s e NNffiO
MM- 3
LL=4
Adv. Eng. Software, 1987, VoL 9, No. 3
147
Result NN=O HH-3 LL=4 ICELH-IO
7
NODE-I ELEMENTI NODE INCIDENCE I 3 2 2 3 5 2 5 4 3 6 5 4 5 8 5 6 10 4 8 ? 5 9 8 5 10 9 6 11 10 IEDBL 1 3 6 11 0 11 I0 9 8 7 7 4 2 1 0 O 0 0 0 0
~
q
4
9 8
1
6
3
II
APPENDIX 2. QUADRILATERAL AREA ~N=3
MM=~
LL=5
Resuic NN=3 MM=4 LL=5 NELM=28 NODE-2 3 ELL'~ENT NODE INCIDENCE 1 2 5 2 3 6 3 4 7 2 6 5 3 7 6 8 7 5 6 9 6 7 I0 7 8 II 6 10 9 7 11 l0 8 12 11 9 i0 16 I0 11 15 9 14 13 I0 15 14 ii 16 15 12 17 16 13 14 19 14 15 20 15 16 21 16 17 22 13 19 18 14 20 19 15 21 20 16 22 21 17 23 22 IEDBL 1 2 3 4 0 0 4 8 12 17 23 0 23 22 21 20 19 18 18 13 9 5 1 0
18
19
20
21
22
23
"-.,... 2
12
3
4
APPENDIX 3. 1000 1010 1020 1030 1040 1050 1060 1070 1080 1090 1100 1110 1120 1130 1140 1150 1160 1170 1180 1190
148
REM FILE NAME=FLUSA KEN MESH PLUS REM DIM MTO(lOO.3).MTJBL(100.3).JMAP(100).IEDBL(4.100} INPUT NN INPUT MM INPUT LL LPRINT "NN="=NN: LPRINT "MM=":MM: LPRINT "LL="=LL I F NN=O THEN *SET1 ELSE *SET2 *SETI IXXX=I NN=I MM=MM-1 SOTO *STAT *SET2 IXXX=O *STAT IF LL<=NN*2"MM THEN 1200 LPRINT "INADEQUATE DISSECTION IS REQUIRED." END
Adv. Eng. Software, 1987, Vol. 9, No. 3
1200 FOR I I = l TO MM 1210 IF 2^II>=LL/NN THEN 1230 1220 NEXT I I 1230 I P = I I 1240 IF 2^IP=LL/NN THEN KK=IP ELSE KK=IP-I 1250 FOR I I = l TO KK+I 1260 JMAP(II)=NN*2^(II-1) 1270 NEXT I I 1280 I F KK=MN THEN 1380 1290 IA=INT((LL-JMAP(KK÷I))/(MM-KK)} 1300 IB=LL-(MM-KK)*IA-JMAP(KK÷I) 1310 NELM=O 1320 JNO=O 1330 JMAP(O)=O 1340 FOR II=KK÷2 TO MM+I 1350 IF II>;~M+I-IB THEN IP=l ELSE IP=O 1360 JMAP(II)=JMAP(II-1)+IA+IP 1370 NEXT I I 1380 1390 FOR I I = 1 TO MM 1400 N=JMAP(II) 1410 L=JMAP(II÷I) 1420 GOSUB *PLUSA 1430 I F I I > 1 THEN 1480 1440 FOR dd=l TO NELMI 1450 FOR KK=I TO 3 1460 MTd(JJ.KK}=MTJBL(JJ.KK) 1470 NEXT KK:NEXT dJ :88TO 1530 1480 FOR JJ=NELM+I TO NELM+NELM1 1490 FOR KK=I TO 1500 MTJ(JJ.KK)-MTJBL(JJ-NELM.KK)+JNO-JMAP(II)-I 1510 NEX'T KK 1520 NEXT dd 1530 NELH=NELM+NELMI 1540 JNO=MTJ(NELM.2) 1550 NEXT I I 1560 1570 FOR I f = I TO NN+I 1580 IEDBL(1.II)=II 1590 NEXT I I 160(} 1&I0 IEDBL(2.1)=NN+I 1620 FOR I f = 2 TO MM+I 1630 IEDBL{2.1I)'IEDBL(2.II-I)+JMAP(II)+I 1640 NEXT I I 1650 1660 IEDBL(3.1)=IEDBL(2.MM+I) 1670 FOR I I = 2 TO LL+I 1680 IEDBL(3.II)=IEDBL(3.II-I)-I 1690 NEXT I I 1700 1710 IEDBL(4,1)=IEDBL(3.LL+I) 1720 FOR I f = 2 TO MM÷I 1730 IEDBL(4.II)=IEDBL(4.II-I)-JMAP(NIq+2-II)-I 1740 NEXT I I 1750 NODE=O 1760 FOR I I = l TO MM÷I 1770 NODE=NODE+JMAP(II)+I 178,) NEXT I I 1790 IF IXXX=O THEN *OUTT 1800 FOR I I = l TO NELM 1810 FOR JJ"1 TO 3 1820 NTJ(NELM÷2-1I.JJ)=MTJ(NELM÷I-II.JJ)÷I 1830 P~XT Jd 1840 NEXT I I 1850 M T J ( 1 . I ) = I 1860 MTJ(1.2)=3 1870 MTJ(1.3)=2 1880 FOR JJ=1 TO MM+I 1890 IEDBL(2.NM+3-JJ)=IEDBL(2.MN÷2-JJ)+i 1900 NEXT dJ 1910 IEDBLt2.1)=I 1920 FOR JJ=l TO LL+I 1930 IEDBL(3.JJ)=IEDSL(3.JJ)+I 1940 NEXT JJ 1950 FOR dd=l TO MM÷I 1960 IEDBL(4.JJ)=IEDBL(4.dJ)÷I 1970 NEXT JJ 1980 IEDBL(4.MN+2)=I 1990 FOR I I = l TO 3 2000 FOR JJ=1 TO LL÷I 2010 IEDBL(II.JJ)=IEDBL(II+I.JJ) 2020 NEXT Jd 2030 NEXT I I 2040 FOR JJ=l TO MM÷2 2050 IEDBL(4.JJ)=O 2060 NEXT JJ 2070 NELM=NELM+I 2080 NODE=NODE÷t 2090 *OUTT 2100 LPRINT "NELM=":NELM 2110 LFRINT "NODE=":NGDE 2120 LPRINT "ELEMENT NODE INCIDENCE" 2130 FOR I = I TO NELM 2140 FOR J=1 TO 3 2150 LFRINT USING " # # # # " : M T J ( I . J ) : 2160 NEXT J:LPRINT 2170 NEXT I 2180 2190 LPRINT "IEDBL" 2200 FOR I I = 1 TO 4 2210 FOR JJ=l TO LL÷I 2220 LPRINT USING "~@@#"~IEDBLtII,JJ)= 2230 NEXT JJILPRINT 2240 NEXT I I
2250 2260 2270 2280 2290 2300 2310 2320 2330 2340 2350
" END *PLUER REM P L U S A REM REM DIM MTJ(IO0,3) " "INPUT N "INPUT L IC=O IF L>2*N THEN *COMM
2360 2370 FOR I=l TO N 2380 HTJBL(I.1)=I 2390 MTJBL(I.2)=I÷I 2400 MTdBL(I.3)-N~I÷I 2410 2420 MTJBL(N+I.1)=I+I 2430 MTJBL(N+I,2)=N+I+2 2440 MTdBL(N+I.3)=N+I+i 2450 NEXT I 2460 • 2470 I F L-N=O THEN *FOL 2480 *INC 2490 FOR I = I ÷ I C TO N 2500 MTdBL(I.3)=MTJBL(I.3)+I 2510 NEXT I 2520 "
2530 2540 2550 2560 2570 2580 2590 2600 26t0 2620 2630 2640 2650 2660 2670 2680 2690 2700 2710 2720 2730 2740 2750 2760 2770 2780 2790 2800
FOR I=1 TO N-IC I F N-IC=O THEN *COMM FOR J=l TO 3 MTJBL(2*N+2+IC-I.J)=MTdBL(2*N+I+IC-I.3) NEXT d NEXT I " FOR I=N÷2÷2*IC TO 2*N+I+IC FOR d = 2 TO 3 MTdBL(I.J)=MTJBL(I°J)+I NEXT d NEXT I ' MT3BL(N+2*IC+I,I)=MTJBL(N+2*IC+I.I)-I IC=IC÷I IF IC-L-N THEN *FOL ELSE *INC *FOL NELMI=2*N÷IC 'LPRINT "MTJBL":II "FOR I=1 TO 2*N+IC " FOR 0 = 1 TO 3 " LPRINT USINS " # # # # " I M T J B L ( I , d ) : ' NEXT 3 : L P R I N T "NEXT I 80TO 2800 *COMM PRINT °'INADEQUATE DISSECTI0N IS REQUIRED. '° RETURN
Adv. Eng. Software, 1987, Vol. 9, No. 3
149
Department of Ovil Engineering, School of Engineering, Okayama University, Okayama, Japan
The object of this paper is to propose a new method to generate the topological properties among the input data neceesary for the user of the Finite Element Method. The method can generate the element-node relation for trlang-l-r and quadrilateral areas which are required to subdivide into different numbers of segments on their edges. This flexibility can largely remove the complexity of the mesh generation of two-dimensional area. The algorithm of the proposed method is so simple as to be used on a microcomputer. The main part of program listing written in BASIC is attached to this paper.
INTRODUCTION The term 'Finite Element Method' generally indicates the analysis, but its program consists of the following three parts: Pre-Processor, Analysis, and Post-Processor. Preprocessor works to decrease the tedious works of the FEM user at the preparation of all the input data of Analysis part. Analysis part works in accordance with the input data which are generated at Pre-processor, and the results are governed by the data generated at Pre-processor. This fact indicates that a better pre-processor can not only save the cost for the preparation of the input data but also give better results and improve the reliability of the results. Input data required at Analysis are classified into three kinds; the topological, geometrical, and physical properties. For example, the element-node relation, the co-ordinates of nodes, and Young's modulus of elements belong to the first, the second and the third property, respectively. A slight mistake of the input of the topological property may change the structure being treated as a whole, and, on the other hand, the influence of the mistake of the latter two properties may be restricted only in a part of the structure. This fact clarifies the reason of the requirement of the automatic mesh generator. We consider the mesh generation of a two-dimensional area with complicated boundary configuration which is consisted of several subareas with different physical properties. The first problem is how we can recognise the geometrical configurations of subareas, and the second problem is how they are subdivided into finite elements. Accepted Match 1987. Discussion closes September 1987.
142
Adv. Eng. Software, 1987, Vol. 9, No. 3
The Quad-Tree Method can resolve these problems at the same time, but it requires a number of operations to fit the boundary configurations.* The same difficulty is found at the application of the Mapping Function Method. 2 The Node Pattern Method has no difficulty on this aspect, but it requires long execution-time for meshing of each subarea.a The easiest way is the Blocking Method which subdivides the area into triangular and/or quadrilateral subareas and gives a mesh pattern for each subarea. The Blocking Method is easy to make a program of and to use, but it may become stiff and can not satisfy its user unless various types of mesh patterns are not prepared. The subdivision of three edges of a triangular area into the same number of segments and to subdivide two edges of a quadrilateral area into different numbers are systematic, and their programmings are easily obtained. The combination of these systematic mesh patterns may often give troubles to its user, because the user may fail into the case where all edges of a triangular or quadrailateral subarea must be subdivided into a different number of segments.
IINPUT DATAI number
of subareas
configuration mesh pattern
of subareas of subareas
number
of edges
number
of nodes
subarea-edge number
of segments
relation on each edge
x and y coordinates
[SUBDIVISION
OF AREA INTO SIMPLE
[GENERATION
[CONNECTION
Figure 1.
of nodes
SUBAREAS]
OF MESH FOR EACH SUBAREA[
OF SUBAREAS
INTO ORIGINAL
CONFIGURATION]
Flow-chart of the Blocking Method
0141-1195/871030142-08 $2.00 © 1987 Computational Mechanics Publications
N
J
M N
N
pattern for each of A(i), and the final stage joints them together to form the original configuration. The flow-chart of the Blocking Method is illustrated in Fig. 1. Input data necessary for the Blocking Method are (1) number of subareas, (2) number of nodes and edges, (3) block-node or block-edge relation, (4) mesh pattern of each subarea, (5) number of segments on each edge, (6) x and y co-ordinates of nodes, and some other information. Nodes and edges must be given so that the boundary configuration of all subareas can be sufficiently drawn by using them. Since the input-data of (3), (4), and (5) is concerned with the number of subareas, the amount of the input data is largely governed by the number of subareas. It suggests that the decreasing of the number of subareas is important for the effective use of the Blocking Method. Since each subareas is transformed into a meshed system by the introduction of mesh patterns prepared in the method, the flexibility of the Blocking Method is wholly governed by the kinds of mesh patterns prepared there. Conventional systematic meshing methods shown in Fig. 2 can be easily prepared, but the method with only this kind of mesh generator becomes very stiff. Its flexibility can be improved by preparing the mesh generators for transition region. At the same time the ones for the transition region can decrease the number of subareas required for appropriate mesh generation of the total area. We can easily find several systematic methods for this purpose, and some of them are shown in Fig. 3. Though they are useful for the increase of the flexibility, they still require some restriction of their use. The most preferable and desirable mesh genera-
L (=N+M)
N Figure 2.
Regular mesh patterns
The Blocking Method may become flexible and easy to use when this kind of mesh generator is prepared. The purpose of this study is to propose a method which can divide all edges of a triangular and quadrilateral area into different number of segments. Since the method is derived as an extension of conventional systematic mesh generation methods, it can sufficiently work on a microcomputer. Therefore, the method presented in this paper is developed on a Japanese microcomputer, PC-9801 of the Nippon Electric Company, and the program is written in BASIC.
M
N L
(=N.2M)
AkA/t/ /tA/t M
BLOCKING METHOD Let A be a two-dimensional area with complex boundary configuration. The Blocking Method consists of three stages: at the first stage, divide A into a gathering of simple triangular and/or quadrilateral subareas by connecting nodes on the boundary and also additional nodes placed in A. We note these subareas by A(/), where i indicates the ith subarea. The second stage generates an appropriate mesh
Figure 3.
N Mesh patterns for transition region
Adv. Eng. Software, 1987, Vcl. 9, No. 3
143
tor is the one which can subdivide all edges of triangular and quadrilateral areas into different number o f segments. The Blocking Method becomes really flexible by the introduction o f this kind of meshing method.
L
K
M
FLEXIBLE MESH GENERATOR In this section we propose a new flexible mesh generator applied to a triangular or quadrilateral area which is required to subdivide all edges into different number of segments. The most preferable method is a systematic one as shown in Figs. 2 and 3, because this type of meshing method requires relatively short execution-time and their programmings have simple structure. By this reason we aim to propose a new one as the extension of them. The method consists of two steps, i.e. the generation of topological properties and the modification of the geometrical properties.
(a)
N
L
M
M
(1) Generation of topological properties Firstly we consider the mesh generation for a triangular area. Let A be a triangular area whose three edges are to be subdivided into N, M, and L segments as shown in Fig. 4a, where we assume N > M > L . By the addition of an intersection line the triangle is divided into two triangles whose two edges have the same number of segments (see Fig. 4b).
K-M b)
N
L
M M
K-M N
(c) Figure 5.
M (b)
Further subdivision of these two triangles as shown in Fig. 4c clarifies that the flexible mesh generation method for a triangular area requires a mesh generator for a quadrilateral area a pair of whose opposite edges are to be divided into the same number of segments. The case of a quadrilateral area is explained in Fig. 5. We assume K > M, and N is the minimum among them. The original area is divided into a quadrilateral area, a pair of whose opposite edges have the same number o f segments, and a triangle area whose two edges have the same number of segments (see Fig. 5b). And the addition of an intersection line for the triangle clarifies that the flexible mesh generator for a quadrilateral area can be designed by using the same mesh generator as was required for the triangular
M-1
N-M
Figure 4.
I
(•) / N-M
N-I
M-I (c) Subdivision of triangular area
144 Adv. Eng. Software, 1987, Vol. 9, No. 3
Subdivision of quadrilateral area
case.
1
Now, consider the mesh generation method for a quadrilateral area a pair of whose edges are required to have the same number o f segments. The area may be considered as a gathering of M narrow subareas as shown in Fig. 6. The above consideration leads to the conclusion that if we find a method to give a mesh pattern for a narrow quadrilateral area as shown in Fig. 6, the flexible mesh generator for triangular and quadrilateral areas is completed.
L
M
M
N
11
l,
I1
e
11
I1 d
STEP 4: (We assume M = K.) Set N N = N, MM = M, and LL = L, and go to *MESH. Go to STEP 15. STEP 5: Find edges with the maximum and minimum segments. STEP 6: Divide the area into triangular subareas as shown in Fig. 4. STEP 7: Apply *MESH for these two subareas. STEP 8: Joint these subareas into the original form. Go to STEP 15. (Quadrilateral Cases) STEP 9: If a pair of opposite edges have the same number of segments, then go to STEP 10 else to STEP 11. STEP 10: (We assume N = L.) Find the minimum number among M and K. (We assume M < K . ) Set N N = M, LL = K , and M M = N . And apply *MESH. Go to STEP 15. STEP 11: Find the edge with the minimum segments. (we assume M is the minimum.)
-i-
11
d c
I1 N+2
N+3
N+4
N+5 2N+I
2N+2
-p
b
11 11
a
+ a
i1
1
2
3
4
N+5
N+~
N
11
N Figure 6
Subdivision o f quadrilateral area
N+2 N+3 N+4
I ,' A narrow quadrilateral area must be divided so that a pair of its opposite edges have a different number of segments. A procedure for this purpose is explained in Fig. 7. At first, the area is systematically subdivided and as many nodes as the number of the difference of segments on these two edges are generated on the upper edge. The generated nodes are systematically connected to nodes on the lower edge. For example, if the upper edge is required to have more two nodes than the lower edge, the third figure from the top in Fig. 7 corresponds to it. Mesh generation according to the different ordering of the additional node generation on the upper edge is explained in Fig. 8. The application of the proposed procedure to a different mesh pattern is also illustrated in Fig. 9. The repetition of this procedure to each narrow subarea in Fig. 6 finishes t h e mesh generation of a quadrilateral area, and its application to subareas in Figs. 4 and 5 completes the design of the flexible mesh generator for triangular and quadrilateral areas. The algorithm of this method is expressed as following: STEP 1 : Input anticlockwisely the number of segments of edges of the area, namely N, M, L, and K. Set N = 0 for a triangular case. STEP 2: I f N = 0 then go to STEP 3 else to STEP 9. (Triangular Cases) STEP 3: If two edges have the same number of segments, then go to STEP 4 else to STEP 5.
l" 1
', I 2
\
X,l
I 3
N+2 N+3 N+4 N+5 N+6
I
2
3
\
X,l
Figure 7.
2
3
2 +3
I
4
N
N+7 2N+3
4
N
N+2 N+3 N+4 N+5 N+6 N+7 N+8
i
N+I
4
N+I
2N+4
N+I
3N 3N+I 3N+2
N
N-1
Systematic element generation
Adv. Eng. Software, 198 7, Vol. 9, No. 3
145
routine works to generate a mesh pattern of a triangular or quadrilateral subarea which satisfies the condition that two edges have the same number of segments for a triangular case and a pair of opposite edges have the same number of segments for a quadrilateral case. Its procedures are summarised as following. (See Figs. 6 and 7) STEP 1 : Judgement of the configuration of the subarea. STEP 2: Determination of the number of segments on upper and lower edges of narrow quadrilateral subareas. STEP 3: Go to *PLUSA. STEP 4: Calculation o f nodes on surrounding edges.
+
+
Figure 8. ordering
Mesh generation by different node generation
The subroutine named *PLUSA is the mesh generator for a narrow quadrilateral subarea which is shown in Fig. 7. The program listing of *MESH including *PLUSA written in BASIC is given in Appendix. The Appendix also includes the explanation of input and output data of the method. It must be noted that this subroutine is valid for the mesh generation of the subarea which satisfies the condition of L < N*2 M, where L and N are the maximum and minimum number of segments of opposite edges of a quadrilateral area, and M is the number of narrow areas in Fig. 6. The results of the application of this method are shown in Fig. 10. Figure 10a shows a triangular area whose three edges are divided into four, five, and six segments, and Fig. 10b is the result of a quadrilateral area whose edges are divided into six, five, four, and seven segments.
~ m
(a)
Figure 9.
Mesh generation for different mesh pattern
STEP 12: Compare two edges adjacent to the edge with M segments, and divide the longer edge so that the area is divided into triangular and quadrilateral subareas as shown in Fig. 5. STEP 13: Apply *MESH for these subareas. STEP 14: Joint them so that they form the original one. STEP 15: Mesh generation is finished.
......... ' --__. '--.
i "--.
....~-~ . . . . . . . . . ~,. ,-" ._ ~ , £ S . ~.~-~-_ -3~_.- --- . ..... ,
"\
/
"
-7
. .-,
(b)
As obvious from the above procedure, the main part of the procedure is the subroutine called *MESH. This sub-
146 Adv. Eng. Software, 1987, Vol. 9, No. 3
Figure 10.
Application of flexible mesh generation
.--"
(2) Improvement o f node co-ordinates The mesh generator mentioned above aims mainly at the generation of the topological properties and it includes the node placement method which can set a series of nodes with the same distance between them as shown in Fig. 7. Therefore, the node co-ordinates generated by the method generally can not satisfy its user except in the cases as presented in Figs. 2 and 3. They must be modified. We can i'md a number of methods for the settlement of nodes in the region and in this paper we apply Laplacian mesh generation} In this method, the interior nodes are positioned so that the position vector, P(i), of node 'i' satisfies the equation 1
n
e(i) =-~n ff.=l Jr03 + e(k)J 'n' is the number of elements adjacent to/-node, and P(/) and P(k) are the position of adjacent nodes in the neighbouring element. The results of the application of this Laplacian mesh generation to the meshes in Fig. 10 are given in Fig. 11, and they show that the positions of the nodes are sufficiently improved by the method.
/\/\ i1
\ /;ix
\
~,
CONCLUDING REMARKS In this paper the author proposed a new flexible mesh generator which is an effective tool for generating the transition mesh pattern of a triangular and quadrilateral area. The proposed method requires less restriction compared to the conventional mesh generators for the transition region, and its introduction to the Blocking method can remove the difficulty of the setting of subareas and it can decrease the number of subareas in the area. Finally it leads to the decreasing of the amount of input data for the mesh generation. A more important conclusion obtained from this study is that the proposed method itself can be a new blocking method. Since the proposed method can generate a mesh pattern by giving an arbitrary number of segments for edges, most of the systematic mesh patterns can be prepared only by the proposed method. Thus, the Blocking Method may have only one subroutine for generating a mesh pattern. The proposed method aims to generate all of the topological properties for Finite Element Analysis, and the position of nodes should be modified especially for nodes in the transition region. In this present paper, the Laplacian method was introduced for this purpose and it was shown that the appropriate positioning of the nodes can be easily obtained by the method. The strategy presented in this paper suggests that the direct extension of the proposed method can make it possible to treat N-polygonal region as one subarea. The application of the proposed method after the subdivision of N-polygonal area into the appropriate triangular and/or quadrilateral subregions may complete the mesh generation for the area, and for this useful mesh generator an automatic subdivision method of N-polygon into subareas must be predpared.
REFERENCES
I.--'/.11
.I.--- ,, _.--%---____ "-....\ \..-\, .._-------... (a)
j
\
i- \ .,_1. ._.~
\
\, _ . . i
~-''~
_-
--'-----_.._ /
-
!1-'\ ~
" ,S.-
-"~A_---
/
, /
i ~--.. / "~_
'
-"---- -%.. /
/
I I
- " - - - . 2i - , . .
APPENDIX 1. PROGRAM LISTING FOR GENERATING ELEMENT-NODE RELATION
- _.._. ---..
---.,,./
Input-Data
- - - ....
, /
~
--~'~-,, \
"-,.......,-
(b)
Figure 11.
--~-.--_
------- ,,
! ~_
!-~, ~---"1-1
1 Yerry, M. A. and Shephard, M. S. A modified quadtree approach to finite element generation, 1EEE Computer Graphics and Applications 1983, 3 (1), 39 2 Haber, R. B., Shephard, M. S., Abel, J. F., Gallagher, R. H. and Greenberg, D. P. A general two dimensional finite element pre-processor utilizing discrete transfinite mapping, Int. J. Numer. Meth. in Eng. 1981, 17 (7), 1015 3 Okuda, O. Automatic triangular mesh generation of arbitrary planar domains utilizing the node-pattern for finite element analysis, J. Japan Society of Precision Engineering 1985, 51 (4), 802
Modification of node co-ordinates
--,--~_ --S.~
/-- .--"
N,t
Numbers o f
LL
(For Triangular Case NN~0)
Elements
on Three Edges
Output-Da~a NELM
Number of Elements
NODE
Number of Nodes
MTJ
Element-Node Relation
IEDBL
Nodes on Surrounding Edges
EXAMPLES (I) Trlangular C a s e NNffiO
MM- 3
LL=4
Adv. Eng. Software, 1987, VoL 9, No. 3
147
Result NN=O HH-3 LL=4 ICELH-IO
7
NODE-I ELEMENTI NODE INCIDENCE I 3 2 2 3 5 2 5 4 3 6 5 4 5 8 5 6 10 4 8 ? 5 9 8 5 10 9 6 11 10 IEDBL 1 3 6 11 0 11 I0 9 8 7 7 4 2 1 0 O 0 0 0 0
~
q
4
9 8
1
6
3
II
APPENDIX 2. QUADRILATERAL AREA ~N=3
MM=~
LL=5
Resuic NN=3 MM=4 LL=5 NELM=28 NODE-2 3 ELL'~ENT NODE INCIDENCE 1 2 5 2 3 6 3 4 7 2 6 5 3 7 6 8 7 5 6 9 6 7 I0 7 8 II 6 10 9 7 11 l0 8 12 11 9 i0 16 I0 11 15 9 14 13 I0 15 14 ii 16 15 12 17 16 13 14 19 14 15 20 15 16 21 16 17 22 13 19 18 14 20 19 15 21 20 16 22 21 17 23 22 IEDBL 1 2 3 4 0 0 4 8 12 17 23 0 23 22 21 20 19 18 18 13 9 5 1 0
18
19
20
21
22
23
"-.,... 2
12
3
4
APPENDIX 3. 1000 1010 1020 1030 1040 1050 1060 1070 1080 1090 1100 1110 1120 1130 1140 1150 1160 1170 1180 1190
148
REM FILE NAME=FLUSA KEN MESH PLUS REM DIM MTO(lOO.3).MTJBL(100.3).JMAP(100).IEDBL(4.100} INPUT NN INPUT MM INPUT LL LPRINT "NN="=NN: LPRINT "MM=":MM: LPRINT "LL="=LL I F NN=O THEN *SET1 ELSE *SET2 *SETI IXXX=I NN=I MM=MM-1 SOTO *STAT *SET2 IXXX=O *STAT IF LL<=NN*2"MM THEN 1200 LPRINT "INADEQUATE DISSECTION IS REQUIRED." END
Adv. Eng. Software, 1987, Vol. 9, No. 3
1200 FOR I I = l TO MM 1210 IF 2^II>=LL/NN THEN 1230 1220 NEXT I I 1230 I P = I I 1240 IF 2^IP=LL/NN THEN KK=IP ELSE KK=IP-I 1250 FOR I I = l TO KK+I 1260 JMAP(II)=NN*2^(II-1) 1270 NEXT I I 1280 I F KK=MN THEN 1380 1290 IA=INT((LL-JMAP(KK÷I))/(MM-KK)} 1300 IB=LL-(MM-KK)*IA-JMAP(KK÷I) 1310 NELM=O 1320 JNO=O 1330 JMAP(O)=O 1340 FOR II=KK÷2 TO MM+I 1350 IF II>;~M+I-IB THEN IP=l ELSE IP=O 1360 JMAP(II)=JMAP(II-1)+IA+IP 1370 NEXT I I 1380 1390 FOR I I = 1 TO MM 1400 N=JMAP(II) 1410 L=JMAP(II÷I) 1420 GOSUB *PLUSA 1430 I F I I > 1 THEN 1480 1440 FOR dd=l TO NELMI 1450 FOR KK=I TO 3 1460 MTd(JJ.KK}=MTJBL(JJ.KK) 1470 NEXT KK:NEXT dJ :88TO 1530 1480 FOR JJ=NELM+I TO NELM+NELM1 1490 FOR KK=I TO 1500 MTJ(JJ.KK)-MTJBL(JJ-NELM.KK)+JNO-JMAP(II)-I 1510 NEX'T KK 1520 NEXT dd 1530 NELH=NELM+NELMI 1540 JNO=MTJ(NELM.2) 1550 NEXT I I 1560 1570 FOR I f = I TO NN+I 1580 IEDBL(1.II)=II 1590 NEXT I I 160(} 1&I0 IEDBL(2.1)=NN+I 1620 FOR I f = 2 TO MM+I 1630 IEDBL{2.1I)'IEDBL(2.II-I)+JMAP(II)+I 1640 NEXT I I 1650 1660 IEDBL(3.1)=IEDBL(2.MM+I) 1670 FOR I I = 2 TO LL+I 1680 IEDBL(3.II)=IEDBL(3.II-I)-I 1690 NEXT I I 1700 1710 IEDBL(4,1)=IEDBL(3.LL+I) 1720 FOR I f = 2 TO MM÷I 1730 IEDBL(4.II)=IEDBL(4.II-I)-JMAP(NIq+2-II)-I 1740 NEXT I I 1750 NODE=O 1760 FOR I I = l TO MM÷I 1770 NODE=NODE+JMAP(II)+I 178,) NEXT I I 1790 IF IXXX=O THEN *OUTT 1800 FOR I I = l TO NELM 1810 FOR JJ"1 TO 3 1820 NTJ(NELM÷2-1I.JJ)=MTJ(NELM÷I-II.JJ)÷I 1830 P~XT Jd 1840 NEXT I I 1850 M T J ( 1 . I ) = I 1860 MTJ(1.2)=3 1870 MTJ(1.3)=2 1880 FOR JJ=1 TO MM+I 1890 IEDBL(2.NM+3-JJ)=IEDBL(2.MN÷2-JJ)+i 1900 NEXT dJ 1910 IEDBLt2.1)=I 1920 FOR JJ=l TO LL+I 1930 IEDBL(3.JJ)=IEDSL(3.JJ)+I 1940 NEXT JJ 1950 FOR dd=l TO MM÷I 1960 IEDBL(4.JJ)=IEDBL(4.dJ)÷I 1970 NEXT JJ 1980 IEDBL(4.MN+2)=I 1990 FOR I I = l TO 3 2000 FOR JJ=1 TO LL÷I 2010 IEDBL(II.JJ)=IEDBL(II+I.JJ) 2020 NEXT Jd 2030 NEXT I I 2040 FOR JJ=l TO MM÷2 2050 IEDBL(4.JJ)=O 2060 NEXT JJ 2070 NELM=NELM+I 2080 NODE=NODE÷t 2090 *OUTT 2100 LPRINT "NELM=":NELM 2110 LFRINT "NODE=":NGDE 2120 LPRINT "ELEMENT NODE INCIDENCE" 2130 FOR I = I TO NELM 2140 FOR J=1 TO 3 2150 LFRINT USING " # # # # " : M T J ( I . J ) : 2160 NEXT J:LPRINT 2170 NEXT I 2180 2190 LPRINT "IEDBL" 2200 FOR I I = 1 TO 4 2210 FOR JJ=l TO LL÷I 2220 LPRINT USING "~@@#"~IEDBLtII,JJ)= 2230 NEXT JJILPRINT 2240 NEXT I I
2250 2260 2270 2280 2290 2300 2310 2320 2330 2340 2350
" END *PLUER REM P L U S A REM REM DIM MTJ(IO0,3) " "INPUT N "INPUT L IC=O IF L>2*N THEN *COMM
2360 2370 FOR I=l TO N 2380 HTJBL(I.1)=I 2390 MTJBL(I.2)=I÷I 2400 MTdBL(I.3)-N~I÷I 2410 2420 MTJBL(N+I.1)=I+I 2430 MTJBL(N+I,2)=N+I+2 2440 MTdBL(N+I.3)=N+I+i 2450 NEXT I 2460 • 2470 I F L-N=O THEN *FOL 2480 *INC 2490 FOR I = I ÷ I C TO N 2500 MTdBL(I.3)=MTJBL(I.3)+I 2510 NEXT I 2520 "
2530 2540 2550 2560 2570 2580 2590 2600 26t0 2620 2630 2640 2650 2660 2670 2680 2690 2700 2710 2720 2730 2740 2750 2760 2770 2780 2790 2800
FOR I=1 TO N-IC I F N-IC=O THEN *COMM FOR J=l TO 3 MTJBL(2*N+2+IC-I.J)=MTdBL(2*N+I+IC-I.3) NEXT d NEXT I " FOR I=N÷2÷2*IC TO 2*N+I+IC FOR d = 2 TO 3 MTdBL(I.J)=MTJBL(I°J)+I NEXT d NEXT I ' MT3BL(N+2*IC+I,I)=MTJBL(N+2*IC+I.I)-I IC=IC÷I IF IC-L-N THEN *FOL ELSE *INC *FOL NELMI=2*N÷IC 'LPRINT "MTJBL":II "FOR I=1 TO 2*N+IC " FOR 0 = 1 TO 3 " LPRINT USINS " # # # # " I M T J B L ( I , d ) : ' NEXT 3 : L P R I N T "NEXT I 80TO 2800 *COMM PRINT °'INADEQUATE DISSECTI0N IS REQUIRED. '° RETURN
Adv. Eng. Software, 1987, Vol. 9, No. 3
149