A triangulation algorithm for arbitrary planar domains

A triangulation algorithm for arbitrary planar domains J. M . Nelson Shock Physics Group, The Boeing Co., Seattle, Washington, U.S.A. (Received 1 August 1977)

For finite element problems which involve irregular regions it is useful to have an automatic scheme capable of forming consistent, high-quality triangulations economically and efficiently. Existing computer algorithms tend either to be lacking in flexibility or inordinately time consuming. This algorithm was developed for the FEHPOL system and has significant advantages over previous methods. A flow chart and computer listing are included. Introduction Although several algorithms have been available for the triangulation of given regions, notably the modified Suhara-Fukudu algorithm and the method of isoparametric coordinates, the necessity to represent the typically, irregularly shaped regions associated with coastal and estuarial problems and the desire for fast interactive response requires a more flexible and less time-consuming algorithm. The method of isoparametric coordinates, while quite effective for regularly shaped regions, is insufficiently flexible for arbitrary domains. The modified Suhara-Fukudu algorithm, while (in the present author's opinion) the best of those in common usage, is unnecessarily complicated and time consuming. In addition, the triangulizations formed by it are not always the most desirable. The algorithm requires a large amount of checking to see if a potential element overlaps a previously generated element or boundary side. By choosing a better suitability criterion, the overlap check may be completely avoided, and the boundary check limited to the most suitable of the potential elements. Because the time required for triangulation is largely due to these consistency checks, elimination of them would be a substantial economy. In choosing a suitability criterion for triangular elements one may consider the following problem: given any two of several points as a baseline, and a half-space on which to seek one of the other points to form a triangle, what is the locus of points of equal suitability such that the same connections result when approached from any other baseline? For instance, consider the points A, B, C and D in Figure 1. If A B is the baseline and D is more suitable than C, we want the-point A to be more suitable than C when the baseline is BD. If this is the case, we can avoid checking to see if the most suitable triangle with

0307-904X 78 0203-0151 $02.00

(" 1978IPC BusinessPress.

baseline BD crosses element ABD. Since, without this criterion, every element must be checked for overlap, it is a considerable time saver. The locus of points equally suitable to form a triangular element with AB, must be symmetrical about the perpendicular bisector of AB since the specification of A and B is orderless. Similarly, the locus of points of equal suitability for BD is symmetric about its own perpendicular bisector. If the locus of points equally suitable for BD is different from those equally suitable for AB then it is possible to pick a point between the two loci such that one of the two baselines will prefer it to A or D. Therefore, they are the same locus. Being the same locus and symmetrical about two arbitrary lines, the locus is a circle. Therefore, the most suitable triangulation is Considered to be the one which, when circumscribed by a circle, the centre of the circle is farthest behind the baseline. That is to say, the signed distance of the centre of the circle from the baseline is a minimum (see Figure 2),

A

,/

\ \

\

\

\

,/

J J ,/

/

/

\ \

J

/

---

/

\

,/

\

/

/

J

\ .'e

D

Figure 1

Appl. Math. Modelling, 1978, Vol 2, September 151

Triangulation algorithm for arbitrary planar domains: J. M. Nelson //----

~

/Less ~

C

! \

\

/1t

\

~

M

o

s

suitable

\

t

baseline is a boundary segment or already occurs in two existing elements, unstack a saved baseline. When the stack is empty the routine is finished. There should be 2(c - 2 + m) - t elements in the final mesh, where c is the number of vertices (corner nodes), m is the number of boundaries (1 for simply-connected regions) and t is the total number of boundary nodes.

suitable

Mostsuitable triangulation

Figure 2

If any n nodes are found to lie on the circle defining the most suitable triangulation, they are all triangulized in fight-hand order since any triangulation combinations would be equally satisfactory. The order of triangulation is as follows. First, starting with any boundary segment, find the most suitable node(s) which can be connected according to the criterion without the resulting element crossing any region boundary. Of the two new baselines created by this connexion(s), stack one baseline and let the other be the new current baseline. Whenever the current

Figure 3 Interactiveresponsetime lessthan 0.1 sec on

CDC 6600

1t5"1

/\

sL~

\

I01 • I'~

/ \/

-O57~

02S

Figure 4

152

Interactive responsetime, 13 sec on Texas9806 minicomputer

Appl. Math. Modelling, 1978, Vol 2, September

Triangulation algorithm for arbitrary planar domains: J. M. Nelson Besides being more efficient than the modified Suhara-Fukudu scheme, the algorithm produces high quality triangulations. There is, of course, no agreedupon criterion for element quality over an entire mesh, however it is clear that the method used does avoid long, thin elements whenever possible. It also avoids the characteristic problem of being forced into a poor connexion by previously generated elements since the triangulations are the same independent of the order of procedure. The routine performing of these calculations has been operational for nearly three years and has proved its effectiveness and versatility in several applications. A flow chart and computer listing are included in the next pages together with some examples of triangulized regions (Figures 3-7).

References 1

Adey, IL et al, Proc. Int. Conf. AppL Numer. Modelling, Univ. of Southampton, July 1977 Brebbia, C. A. and J. J. Connor, "Finite Element Techniques for Fluid Flow', 2nd Edition, Butterworths, London, 1977 Cavendish, J. C. Int. J. Numer. Meth. Eng, 1974, 8, 679 Nelson, J. M. PhD Thesis, Univ. Southampton, 1976. Zienkiewicz, O. C. and Phillips, D. V. hzt. J. Numer. Meth. Eng. 1971, 3, 519

2 3 4 5

Figure5 Mesh generation requiring 12 CPU sec on CDC 6600

SUBROUTINE C C C C C C C C C C C C C C C C C C C C C C

FPCNTR(XI,YI,~2,Y2,X3,Y3,XC,YC,RSQ)

SEGMENT NAME ................................ FPCNTR SEGMENT TYPE ............................ SUBROUTINE SEGMENT LENGTH IN C A R D S . . . . . . . . . . . . . . . . . . . . . . . . 22 DATE OF LAST ALTERATION ................... 20/07/76 NUMBER OF LOCAL INTEGERS ...................... 0 NUMBER OF LOCAL REALS ......................... 8 FEHPOL

ROUTINES

FORTRAN

LIBRARY

-

OTHER

LIBRARY

CALLED ROUTINES

ROUTINES

CALLING

ROUTINES

COMMONS

REFERENCED

REMARKS: THIS ROUTINE IS (X2,Y2),(X3,Y3).

o

o

o

o

o

o

o

o

o

e

CALLED

CALLED

o

o

s

o

o

o

o

o

o

a

o

o

o

o

...............

.................

............................ ............................

CALLED WITH IT RETURNS

NONE NONE

NONE FPINTP FPMESH NONE

THREE POINTS (XI,X1) (XC,YC):THE CENTRE OF

Appl. Math. Modelling, 1978, Vol 2, September

153

Triangu~tiona~or~hm ~ra~itra~p~nardomam~J.M. Nelson C C C C

THE CIRCLE DEFINED BY T H E T H R E E P O I N T S , AND RSQ= THE RADIUS SQUARED. IF T H E Y A R E C O L I N E A R , RSQ=A LARGE NUMBER.

10

20

X2SQ:X2*X2 ~2SQ:Y2*Y2 A:X2-XI A:A+A B:Y2-YI B:B+B C=Y2SQ-YI*YI+X2SQ-XI*XI D=X3-X2 D=D+D E=Y3-Y2 E=E+E F:Y3*Y3-Y2SQ+X3*X3-X2SQ X2SQ=D*B-E*A IF(X2SQ)IO,20,10 XC=(F*B-C*E)/X2SQ YC=(D*C-F*A)/X2SQ RSQ=(XC-X1)**2+(YC-YI)**2 RETURN RSQ:I0000000000000, RETURN

END SUBROUTINE NL,NB,NP) C C C C C C C C C C C C C C C C C C C C C C C C C C

FPMESH(S1,S2,S3,AA,CY,NS,NR,X,Y,M1,M2,NC,

NAME ................................ FPMESH SEGMENT SEGMENT LENGTH IN C A R D S . . . . . . . . . . . . . . . . . . . . . . . . . 126 D A T E OF L A S T A L T E R A T I O N . . . . . . . . . . . . . . . . . . . 21/07/76 NUMBER OF L O C A L I N T E G E R S . . . . . . . . . . . . . . . . . . . . . . 18 N U M B E R OF L O C A L R E A L S . . . . . . . . . . . . . . . . . . . . . . . . . . 13 FEHPOL

FORTRAN OTHER

ROUTINES

CALLED ......................

LIBRARY LIBRARY

ROUTINES ROUTINES

CALLING

ROUTINES

COMMONS

REFERENCED

FPSSOL FPC~TR

CALLED ...............

SQRT IABS NONE

CALLED .................

............................

FPCNCT

............................

NONE

REMARKS: THIS ROUTINE PERFORMS THE CONNECTION OF A F I N I T E ELEMENT MESH GIVEN THE LOCATION OF A L L T H E N O D E S AND IDENTIFICATION OF T H E B O U N D A R Y NODES. ARGUMENTS: SI,S2,S3=WORKAREA

FOR

STACKS,

154 Appl. Math. Modelling, 1978, Vol 2, September

LENGTH

AT

LEAST

MI

"o

E E

0 o.

B" .ca

.m o

O~

o"

m

.

Ill

t~ 0

.

~z: t-4 ~ I ~:11z~, c., H .,,0 L-z'J I

*

*

II

*.-~ *

~

~

*

* *

*

* *

~ H

*

*

*

*

I-~ ~-] ~Z~ 0

I~

~:1

*

*

*

*

~

* *

.

*

*

L"~ Jl

*

~J

H *

*

~

*

," ~

* *

I * *

I

~ * *

*

* *

I ~ O

~ .

*

~ ' ~

Z O

* ~ H ~ * * *

.

~*

*

*

*

*

0 H ~ 0 0 H

*

O

*

~

*

*

*

*

~

~

* ~

~ I

II

*

~ O ~ @ ~ ~ ~ ~ O ~ ~ * ~ 0 ~ , ~ ~ ~ H H ~ ~

~

~

* * * * * * *

* *

~ ~ I

~

* * *

~

*

*

* *

V *

*

*

* *

*

* "o O 'x~ * O ~ H

~

*

D-J t 0

*

*

*

*

*

* . . *

*

I

~0(n I ~ 0 0 LTJ ~'rj C: L--' =Z: Z H

*

IZ~

* . . *

*

*I-4 * *

0

H

,-~ L-'

O r~ C')

,...]

O ~J

,.-]

I-4

L-'

DO

O? :Z

~D

,-] I.-4 O

~

~

~

~

~

~

~

~

~

O

0

~

,

~

~

~

~

~

~

0

~

0

~

~

~

~

~ ~

II II ~ ~

~ ~

~ ~

~

~

0

0

~

0

Z

0 ~ 0 ~

0 0 ~ ~ .

0 ~

~

0

~

~

0 ~

~

~0

0

~

~

~

0

H - O ~ ~

~

~

0

~

~

II II

~

~

~',

0"

(b "0

o_

<

.co

-.I

m. 5"

o

E

E

2

"o

>

o'1

O

ilZ~'II

~

o

O

O

~

H

II ~

H

II O

~

II H

H

I

n H

H

"

•

I 0

0

O

0

H

H

0

O

HH

n H

H

W:

We r.~

lie Wc

IIc w:

I~lk}

l(cO IIC m

Wc

II:

Wc

Wc~

11¢ ~

m

~

IkO

I~E~

IN

I~

IK

~tD

0 1~ 0 wc 0 0

We ~0

i~

I~1-~

o

W~ O

Wc

~ II We

IW

I~1-4

Wc L'~ I~

O 0

II "11

~0 I~

"

~

H

m

~

II O

H

O 0

0

~

~

~

I1~

0

O 0

0

~

O

H

~ H

H

O O O 0 0 ~ ~ 0 0 0 0

.

.

I

n ~

H

0 (D m wc ~N 0 W~ W~ W, {n ~W

O 0

~ 0

~

.

~ ~ ~

~ • ~

•

.

~

.

~

II H

H

•

H

H

,

.

H v v V -

II ~

~

~

H

~ n ~ H

H

O

O

0 0 0 0

0

0

II

~

~

~

II

II

"

~"

"

We

..,.&

V

r~O

I'0

tn

%.n

UJ

C3

L~

--~

~0 +

I

f~ '.-" Ca

~ - - ~ "--" ~-w ~-," II C ~ " Ca ~--" Cn

-

he C~ Ca

~,~ ~

o

5'

3

DJ

W

E

"11

w

)I¢

0

~C~O

~

~ - ~

X¢

C~

i

8

L"z~ m

I0

L"~

O H

WC

H

S

R

__i

o

O~

o

<

00

o

E

E

"0

~+

II

~

~

H

~ ~

0

O

H

~

~

~ 0 ~

~ ~0

H ~ O

H

O ~

O ~

~

~

~

•

•

•

~

0

0

H

~

~

~

~ o

0

0

~

+

II

•

0

0

I

H

0

~ H O ~ ~

•

O

O

~

0

~

. . . . . .

~

H H H H H ~ ~ ~ ~ ~ 0 -

~

o ~

-

-

O

0

H

~ ' O ~

~

~

H

O

~

•

~

O ~

~

II

H 0

~

~+

~

H

o

o

~

.~

~

H

,

H ~

0 ~

~ ~

~ II

0 ,

-

~ ~

0

0

H

H

-

O ~

~ 0

~ 0

o ~

-

0

0

0

0

~

-

0 ~ 0

II

~

O ~ I I ~

~

~

H

,

~ ~

H

~ 0 0 0 0 ~

~

~

H

~ - 0 0 0 0 ~

0

~ ~

P ~

,

~

.

H

~

~

X ~ v v v v H

~

~ 0

~ ~

0

O

~ • ~

~ ~H

I

~

~

II

H II

II

~0

o

0

0

0

c') cJ c )

o~ o

0

0

v

m

I

L,UL~

L~

" t~

--, .I=

~00 0o.

"

"

~-r'

"

~

"

"~H

Z

0

H

I~

~

~

[./~ C }

I

~

L'-'O

0

~

~ 0

~

0

Z

Z

L-'

I-t

L-~

0

~

~

~

.

II

II

•

•

~

II

0

¢-i

~

~

0

O0

C) 0

°

H

I~

~'~ ~.4

~-~ I.-.I -..~ I

•

m

II

~

~

II

|1

II

II

q

o-

8"

"o

t.0

o

<

t..o ,,,4

m

o D. e_

E

E

"'o

f,yl

e

o

6

~

@

s

o

•

o

e

o

•

4

o

@

~

*

o

*

•

e

@

o

o

e

•

o

o

e

e

e

.

~

~

~

o

o

o

o

o

e

e

o

e

o

@

@

e

. •

*

s

~

o

$

e

.

o

e

o

o

o

@

~

~ 0 0 ~ ~

e

o

o

,

e

o

e

o

o

o

~

o

,

*

•

~ , ~

~

~ . .

@

6

@

~

.

.

~

~

o

*

*

e

o

o

e

~ ~ .

e

e

J

@

@

o

$

•

e

•

•

. •

~

~

~

~

O

~

O

~

O0

~

• .

0

O

~

•

Z

H

0

O

~

~

~

II

O

O

--~

O

~

II

~

~

~

~

m

m

•

~

8

8

•

1~ O

~

m

,...-,

~ O

•

~C*

~

$

•

~

N

I--4

~

II

~

II

II

O

O

I--4

I

I

m I~

m

m

I~

~

C~ C~ C~

II

I-4 I--4

[~1

II

II

.. I"4

II

--'~ ,.~ )-4

II

,'~

,~'~ I I

I ~ --~ ~tl ~-'~ O ,~

~

I'-I I-4

0

I-I .O ÷

II

i--4 O

t-I I-4

O~ H

t-t ll~, 0

II

I

v

v

H O

1"4 O

A

X

I-I

t~

v

v

I-I O

I-4 O

A

~

I.~ ,

kn H

v

II

0

0

O

•

~

~

~u I--4 ~u

l-4 I-4 ~

t.~ O

~

o~

1-40

0,

"5

.~

o'1

o"

3

tO

O0

to

co

e.

o o..

0

e-

0

t,O

o

tO

r"

t:::

~

~

~

:>< :><: 0 L"~ ~ - , , , • • --,1

0000 O 0

O 0

~-~ ~-~ O 0

O 0 0o0o

~

•

O 0

•

~

~-] "

~

•

:~,

~

t

~

~..~ C.~

r.~

~

O

I-'t 0

tt

C.~'~

I

O

tt

,.~-..,

I

~

It

:>, :~

I

tt

•

I

It

0

L~J~

• ~.

0

0

•

°

~

0

t

•

--~

~

I! 0

~

0

~ O ~ O .

~

o

It

C~

CJ

¢-3 C~

O0 ~

~:~

~,.

•

•

,

•

"~"J ~-8

•

•

c3 D-J ~ •

~

L'~

~

C~

•

•

.

•

•

•

•

O0 • • •

I-'4

0

t-'

C'~

•

•

•

•

•

t:~

O c= i--3 I-4

~t~

tzJ

0

C'~ C~ ¢."3 C-3 ~."3 C )

•

C~

O ~

I--4

I-4

L-~..j L,-rj ~

•

~:~

i - ] t.-3 ::~ c=~ O D3 ::s:/ t-3 ~

O

0

I-8 ~:x:~ I-.-.I

L'rJ,.

~'J , " ~ ~-3 ,v,~

C3

o z

•

•

t"

L'T'J

0

I-4 t~ m:;I

0

~rj

C'3 C~

.<..

"0

8.