CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
triangulate.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine triangulate (ics, rcs0, zcs0, ncsnodes, rcscg, rcs0cg, zcscg, zcs0cg, nrcg, nzcg, jcs, kontri, straight, ne, ipkon, kon, lakon, lcs, netri, ifacetet, inodface)
 

Function/Subroutine Documentation

◆ triangulate()

subroutine triangulate ( integer, dimension(*)  ics,
real*8, dimension(*)  rcs0,
real*8, dimension(*)  zcs0,
integer  ncsnodes,
real*8, dimension(*)  rcscg,
real*8, dimension(*)  rcs0cg,
real*8, dimension(*)  zcscg,
real*8, dimension(*)  zcs0cg,
integer, dimension(*)  nrcg,
integer, dimension(*)  nzcg,
integer, dimension(*)  jcs,
integer, dimension(3,*)  kontri,
real*8, dimension(9,*)  straight,
integer  ne,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
integer, dimension(*)  lcs,
integer  netri,
integer, dimension(*)  ifacetet,
integer, dimension(*)  inodface 
)
22 !
23 ! generate a triangulation of the independent side (= right side)
24 !
25 ! the element faces of the independent side are identified and
26 ! triangulated. The nodes belonging to the faces are stored in
27 ! field inodface, face after face. For a triangle i the value
28 ! ifacetet(i) points to the last node in field inodface of the
29 ! face the triangle belongs to.
30 !
31  implicit none
32 !
33  character*8 lakon(*)
34 !
35  integer jcs(*),l,j,ics(*),nodef(8),ifacetet(*),
36  & nrcg(*),node,ncsnodes,id,ifaceq(8,6),ifacet(6,4),
37  & ifacew1(4,5),iface(8,6),nodelem(20),nnodelem,nzcg(*),
38  & itrifac3(3,1),itrifac4(3,2),itrifac6(3,4),itrifac8(3,6),
39  & itrifac(3,6),ifacew2(8,5),lcs(*),inodface(*),nnodface,
40  & k,kflag,i,ne,ipkon(*),kon(*),indexe,nope,nface,nodface,jface,
41  & netri,ntrifac,kontri(3,*),m
42 !
43  real*8 straight(9,*),zcscg(*),rcscg(*),zcs0cg(*),
44  & rcs0cg(*),cgl(2),col(2,3),rcs0(*),zcs0(*)
45 !
46 ! nodes per face for hex elements
47 !
48  data ifaceq /4,3,2,1,11,10,9,12,
49  & 5,6,7,8,13,14,15,16,
50  & 1,2,6,5,9,18,13,17,
51  & 2,3,7,6,10,19,14,18,
52  & 3,4,8,7,11,20,15,19,
53  & 4,1,5,8,12,17,16,20/
54 !
55 ! nodes per face for tet elements
56 !
57  data ifacet /1,3,2,7,6,5,
58  & 1,2,4,5,9,8,
59  & 2,3,4,6,10,9,
60  & 1,4,3,8,10,7/
61 !
62 ! nodes per face for linear wedge elements
63 !
64  data ifacew1 /1,3,2,0,
65  & 4,5,6,0,
66  & 1,2,5,4,
67  & 2,3,6,5,
68  & 3,1,4,6/
69 !
70 ! nodes per face for quadratic wedge elements
71 !
72  data ifacew2 /1,3,2,9,8,7,0,0,
73  & 4,5,6,10,11,12,0,0,
74  & 1,2,5,4,7,14,10,13,
75  & 2,3,6,5,8,15,11,14,
76  & 3,1,4,6,9,13,12,15/
77 !
78 ! triangulation for three-node face
79 !
80  data itrifac3 /1,2,3/
81 !
82 ! triangulation for four-node face
83 !
84  data itrifac4 /1,2,4,2,3,4/
85 !
86 ! triangulation for six-node face
87 !
88  data itrifac6 /1,4,6,4,2,5,6,5,3,4,5,6/
89 !
90 ! triangulation for eight-node face
91 !
92  data itrifac8 /1,5,8,5,2,6,7,6,3,8,7,4,8,5,7,5,6,7/
93 !
94 ! pointer into field inodface
95 !
96  nnodface=0
97 !
98 ! sort the nodes on the independent side
99 !
100  do j=1,ncsnodes
101  jcs(j)=abs(ics(j))
102  lcs(j)=j
103  enddo
104 !
105  kflag=2
106  call isortii(jcs,lcs,ncsnodes,kflag)
107 !
108  netri=0
109 !
110 ! check the elements adjacent to the independent nodes
111 !
112  do i=1,ne
113  indexe=ipkon(i)
114  if(lakon(i)(4:4).eq.'2') then
115  nope=20
116  nface=6
117  nodface=8
118  elseif(lakon(i)(4:4).eq.'8') then
119  nope=8
120  nface=6
121  nodface=4
122  elseif(lakon(i)(4:5).eq.'10') then
123  nope=10
124  nface=4
125  nodface=6
126  elseif(lakon(i)(4:4).eq.'4') then
127  nope=4
128  nface=4
129  nodface=3
130  elseif(lakon(i)(4:5).eq.'15') then
131  nope=15
132  nface=5
133  nodface=8
134  elseif(lakon(i)(4:4).eq.'6') then
135  nope=6
136  nface=5
137  nodface=4
138  else
139  cycle
140  endif
141 !
142 ! check which nodes of the element belong to the independent set
143 !
144  nnodelem=0
145  do j=1,nope
146  nodelem(j)=0
147  node=kon(indexe+j)
148  call nident(jcs,node,ncsnodes,id)
149  if(id.le.0) cycle
150  if(jcs(id).ne.node) cycle
151  nodelem(j)=node
152  nnodelem=nnodelem+1
153  enddo
154  if(nnodelem.eq.0) cycle
155 !
156  if(nface.eq.4) then
157  do j=1,nface
158  do k=1,nodface
159  iface(k,j)=ifacet(k,j)
160  enddo
161  enddo
162  elseif(nface.eq.5) then
163  if(nope.eq.6) then
164  do j=1,nface
165  do k=1,nodface
166  iface(k,j)=ifacew1(k,j)
167  enddo
168  enddo
169  elseif(nope.eq.15) then
170  do j=1,nface
171  do k=1,nodface
172  iface(k,j)=ifacew2(k,j)
173  enddo
174  enddo
175  endif
176  elseif(nface.eq.6) then
177  do j=1,nface
178  do k=1,nodface
179  iface(k,j)=ifaceq(k,j)
180  enddo
181  enddo
182  endif
183 !
184 ! check which face of the element belongs to the independent side
185 !
186  jface=0
187  loop: do m=1,nface
188  nnodelem=nodface
189 !
190 ! several faces of one and the same element may belong
191 ! to the master surface
192 !
193  do k=1,nodface
194  if(iface(k,m).eq.0) then
195  nnodelem=k-1
196  exit
197  endif
198  if(nodelem(iface(k,m)).eq.0) cycle loop
199  enddo
200  jface=m
201 c exit
202 c enddo loop
203 c if(jface.eq.0) cycle
204 !
205 ! store the node numbers in a local face field
206 !
207  do k=1,nnodelem
208  nodef(k)=nodelem(iface(k,jface))
209  inodface(nnodface+k)=nodef(k)
210  enddo
211  nnodface=nnodface+nnodelem
212 !
213 ! number of triangles
214 !
215  if(nnodelem.eq.3) then
216  ntrifac=1
217  do j=1,ntrifac
218  do k=1,3
219  itrifac(k,j)=itrifac3(k,j)
220  enddo
221  enddo
222  elseif(nnodelem.eq.4) then
223  ntrifac=2
224  do j=1,ntrifac
225  do k=1,3
226  itrifac(k,j)=itrifac4(k,j)
227  enddo
228  enddo
229  elseif(nnodelem.eq.6) then
230  ntrifac=4
231  do j=1,ntrifac
232  do k=1,3
233  itrifac(k,j)=itrifac6(k,j)
234  enddo
235  enddo
236  elseif(nnodelem.eq.8) then
237  ntrifac=6
238  do j=1,ntrifac
239  do k=1,3
240  itrifac(k,j)=itrifac8(k,j)
241  enddo
242  enddo
243  endif
244 !
245  do j=1,ntrifac
246 !
247 ! new triangle
248 !
249  netri=netri+1
250  do l=1,2
251  cgl(l)=0.d0
252  enddo
253  do k=1,3
254  node=nodef(itrifac(k,j))
255  kontri(k,netri)=node
256  call nident(jcs,node,ncsnodes,id)
257  col(1,k)=rcs0(lcs(id))
258  col(2,k)=zcs0(lcs(id))
259  do l=1,2
260  cgl(l)=cgl(l)+col(l,k)
261  enddo
262  enddo
263 !
264 ! center of gravity of the triangle
265 !
266  rcscg(netri)=cgl(1)/3.d0
267  zcscg(netri)=cgl(2)/3.d0
268 !
269 ! determining the equations of the straight lines bordering
270 ! the triangle
271 !
272  call straighteq2d(col,straight(1,netri))
273 !
274  ifacetet(netri)=nnodface
275 !
276  enddo
277  enddo loop
278  enddo
279 !
280  if(netri.eq.0) then
281  write(*,*) '*ERROR in triangulate: no faces found on the'
282  write(*,*) ' independent side'
283  call exit(201)
284  endif
285 !
286 ! initialization of near2d
287 !
288  do i=1,netri
289  nrcg(i)=i
290  nzcg(i)=i
291  rcs0cg(i)=rcscg(i)
292  zcs0cg(i)=zcscg(i)
293  enddo
294 !
295  kflag=2
296  call dsort(rcscg,nrcg,netri,kflag)
297  call dsort(zcscg,nzcg,netri,kflag)
298 !
299  return
straighteq2d
subroutine straighteq2d(col, straight)
Definition: straighteq2d.f:20
isortii
subroutine isortii(ix, iy, n, kflag)
Definition: isortii.f:6
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
dsort
subroutine dsort(dx, iy, n, kflag)
Definition: dsort.f:6
Hosted by OpenAircraft.com, (Michigan UAV, LLC)