1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19! S.W. Sloan, Adv.Eng.Software,1987,9(1),34-55. 20! Permission for use with the GPL license granted by Prof. Scott 21! Sloan on 17. Nov. 2013 22! 23 subroutine delaun(numpts,n,x,y,list,stack,v,e,numtri) 24! 25 implicit none 26! 27 integer v(3,*),n,i,t,list(*),numtri,p,e(3,*),maxstk,topstk, 28 & v1,v2,v3,l,r,pop,a,b,c,erl,era,erb,edg,triloc,numpts, 29 & tstrt,tstop,stack(*) 30! 31 real*8 x(*),y(*),xp,yp,c00000,c00100 32! 33 logical swap 34! 35 parameter(c00000=0.d0, 36 & c00100=100000.d0) 37! 38 v1=numpts+1 39 v2=numpts+2 40 v3=numpts+3 41 v(1,1)=v1 42 v(2,1)=v2 43 v(3,1)=v3 44 e(1,1)=0 45 e(2,1)=0 46 e(3,1)=0 47! 48 x(v1)=-c00100 49 x(v2)= c00100 50 x(v3)= c00000 51 y(v1)=-c00100 52 y(v2)=-c00100 53 y(v3)= c00100 54! 55 numtri=1 56 topstk=0 57 maxstk=numpts 58 do 100 i=1,n 59 p=list(i) 60 xp=x(p) 61 yp=y(p) 62 t=triloc(xp,yp,x,y,v,e,numtri) 63 a=e(1,t) 64 b=e(2,t) 65 c=e(3,t) 66 v1=v(1,t) 67 v2=v(2,t) 68 v3=v(3,t) 69 v(1,t)=p 70 v(2,t)=v1 71 v(3,t)=v2 72 e(1,t)=numtri+2 73 e(2,t)=a 74 e(3,t)=numtri+1 75! 76 numtri=numtri+1 77 v(1,numtri)=p 78 v(2,numtri)=v2 79 v(3,numtri)=v3 80 e(1,numtri)=t 81 e(2,numtri)=b 82 e(3,numtri)=numtri+1 83 numtri=numtri+1 84 v(1,numtri)=p 85 v(2,numtri)=v3 86 v(3,numtri)=v1 87 e(1,numtri)=numtri-1 88 e(2,numtri)=c 89 e(3,numtri)=t 90! 91 if(a.ne.0) then 92 call push(t,maxstk,topstk,stack) 93 end if 94 if(b.ne.0) then 95 e(edg(b,t,e),b)=numtri-1 96 call push(numtri-1,maxstk,topstk,stack) 97 end if 98 if(c.ne.0) then 99 e(edg(c,t,e),c)=numtri 100 call push(numtri,maxstk,topstk,stack) 101 end if 102! 103 50 if(topstk.gt.0) then 104 l=pop(topstk,stack) 105 r=e(2,l) 106! 107 erl=edg(r,l,e) 108 era=mod(erl,3)+1 109 erb=mod(era,3)+1 110 v1=v(erl,r) 111 v2=v(era,r) 112 v3=v(erb,r) 113 if(swap(x(v1),y(v1),x(v2),y(v2),x(v3),y(v3),xp,yp)) then 114 a=e(era,r) 115 b=e(erb,r) 116 c=e(3,l) 117 v(3,l)=v3 118 e(2,l)=a 119 e(3,l)=r 120 v(1,r)=p 121 v(2,r)=v3 122 v(3,r)=v1 123 e(1,r)=l 124 e(2,r)=b 125 e(3,r)=c 126 if(a.ne.0) then 127 e(edg(a,r,e),a)=l 128 call push(l,maxstk,topstk,stack) 129 end if 130 if(b.ne.0) then 131 call push(r,maxstk,topstk,stack) 132 end if 133 if(c.ne.0) then 134 e(edg(c,l,e),c)=r 135 end if 136 end if 137 goto 50 138 end if 139 100 continue 140 if(numtri.ne.2*n+1) then 141 write(6,'("o***error in subroutine delaun***")') 142 write(6,'(" ***incorrect number of triangls formed***")') 143 call exit(201) 144 end if 145 do 120 t=1,numtri 146 if((v(1,t).gt.numpts).or. 147 & (v(2,t).gt.numpts).or. 148 & (v(3,t).gt.numpts))then 149 do 110 i=1,3 150 a=e(i,t) 151 if(a.ne.0) then 152 e(edg(a,t,e),a)=0 153 end if 154 110 continue 155 goto 125 156 end if 157 120 continue 158 125 tstrt=t+1 159 tstop=numtri 160 numtri=t-1 161 do 200 t=tstrt,tstop 162 if((v(1,t).gt.numpts).or. 163 & (v(2,t).gt.numpts).or. 164 & (v(3,t).gt.numpts))then 165 do 130 i=1,3 166 a=e(i,t) 167 if(a.ne.0) then 168 e(edg(a,t,e),a)=0 169 end if 170 130 continue 171 else 172 numtri=numtri+1 173 do 140 i=1,3 174 a=e(i,t) 175 e(i,numtri)=a 176 v(i,numtri)=v(i,t) 177 if(a.ne.0) then 178 e(edg(a,t,e),a)=numtri 179 end if 180 140 continue 181 endif 182 200 continue 183 end 184 185