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

Go to the source code of this file.

Functions/Subroutines

subroutine initialcfd (nef, ipkonf, konf, lakonf, co, coel, cofa, nface, ielfa, area, ipnei, neiel, xxn, xxi, xle, xlen, xlet, xrlfa, cosa, volume, neifa, xxj, cosb, dmin, ifatie, cs, tieset, icyclic, c, neij, physcon, isolidsurf, nsolidsurf, dy, xxni, xxnj, xxicn, nflnei, iturbulent, rf)
 

Function/Subroutine Documentation

◆ initialcfd()

subroutine initialcfd ( integer  nef,
integer, dimension(*)  ipkonf,
integer, dimension(*)  konf,
character*8, dimension(*)  lakonf,
real*8, dimension(3,*)  co,
real*8, dimension(3,*)  coel,
real*8, dimension(3,*)  cofa,
integer  nface,
integer, dimension(4,*)  ielfa,
real*8, dimension(*)  area,
integer, dimension(*)  ipnei,
integer, dimension(*)  neiel,
real*8, dimension(3,*)  xxn,
real*8, dimension(3,*)  xxi,
real*8, dimension(*)  xle,
real*8, dimension(*)  xlen,
real*8, dimension(*)  xlet,
real*8, dimension(3,*)  xrlfa,
real*8, dimension(*)  cosa,
real*8, dimension(*)  volume,
integer, dimension(*)  neifa,
real*8, dimension(3,*)  xxj,
real*8, dimension(*)  cosb,
real*8  dmin,
integer, dimension(*)  ifatie,
real*8, dimension(17,*)  cs,
character*81, dimension(3,*)  tieset,
integer  icyclic,
real*8, dimension(3,3)  c,
integer, dimension(*)  neij,
real*8, dimension(*)  physcon,
integer, dimension(*)  isolidsurf,
integer  nsolidsurf,
real*8, dimension(*)  dy,
real*8, dimension(3,*)  xxni,
real*8, dimension(3,*)  xxnj,
real*8, dimension(3,*)  xxicn,
integer  nflnei,
integer  iturbulent,
real*8, dimension(3,*)  rf 
)
24 !
25 ! calculating geometric variables of the cells and their faces
26 !
27  implicit none
28 !
29  character*8 lakonf(*)
30  character*81 tieset(3,*)
31 !
32  integer nef,ipkonf(*),konf(*),nface,ielfa(4,*),ipnei(*),neiel(*),
33  & ifaceq(8,6),i,j,k,indexe,kflag,index1,index2,j1,j2,nope,
34  & nodes(4),iel1,iel2,iel3,iface,indexf,neifa(*),nf(5),ifacet(7,4),
35  & ifacew(8,5),numfaces,ied4(2,6),ied6(2,9),ied8(2,12),ifatie(*),
36  & ics,itie,neighface,ifirst_occurrence,icyclic,neij(*),jface,
37  & isolidsurf(*),nsolidsurf,iopp8(4,6),iopp6(3,5),iopp4(1,4),
38  & node,nelem,nflnei,iturbulent
39 !
40  real*8 co(3,*),coel(3,*),cofa(3,*),area(*),xxn(3,*),xxi(3,*),
41  & xle(*),xlen(*),xlet(*),xrlfa(3,*),cosa(*),xsj2(3),xi,et,
42  & shp2(7,4),xs2(3,7),xl2(3,8),xl13,volume(*),dxsj2,xl(3,8),
43  & xxj(3,*),cosb(*),dmin,cs(17,*),xn(3),theta,pi,dc,ds,dd,
44  & c(3,3),diff(3),p(3),q(3),a(3),physcon(*),dy(*),xs(3,3),
45  & aa,bb,cc,dist,xxni(3,*),xxnj(3,*),xxicn(3,*),rf(3,*),x13(3)
46 !
47 ! nodes belonging to the cell faces
48 !
49  data ifaceq /4,3,2,1,11,10,9,12,
50  & 5,6,7,8,13,14,15,16,
51  & 1,2,6,5,9,18,13,17,
52  & 2,3,7,6,10,19,14,18,
53  & 3,4,8,7,11,20,15,19,
54  & 4,1,5,8,12,17,16,20/
55  data ifacet /1,3,2,7,6,5,11,
56  & 1,2,4,5,9,8,12,
57  & 2,3,4,6,10,9,13,
58  & 1,4,3,8,10,7,14/
59  data ifacew /1,3,2,9,8,7,0,0,
60  & 4,5,6,10,11,12,0,0,
61  & 1,2,5,4,7,14,10,13,
62  & 2,3,6,5,8,15,11,14,
63  & 4,6,3,1,12,15,9,13/
64  data nf /3,3,4,4,4/
65  data ied4 /1,2,2,3,3,1,1,4,2,4,3,4/
66  data ied6 /1,2,2,3,3,1,4,5,5,6,6,4,1,4,2,5,3,6/
67  data ied8 /1,2,2,3,3,4,4,1,5,6,6,7,7,8,8,1,1,5,2,6,3,7,4,8/
68  data iopp4 /4,3,1,2/
69  data iopp6 /4,5,6,1,3,2,3,6,0,1,4,0,2,5,0/
70  data iopp8 /5,6,7,8,4,3,2,1,3,4,8,7,4,1,5,8,1,2,6,5,2,3,7,6/
71 !
72  ifirst_occurrence=1
73  icyclic=0
74 !
75 ! coordinates of the center of the cells
76 !
77  do i=1,nef
78  if(ipkonf(i).lt.0) cycle
79  if(lakonf(i)(1:1).ne.'F') cycle
80  indexe=ipkonf(i)
81  if(lakonf(i)(4:4).eq.'8') then
82  nope=8
83  else if(lakonf(i)(4:4).eq.'6') then
84  nope=6
85  else
86  nope=4
87  endif
88  do j=1,3
89  do k=1,nope
90  coel(j,i)=coel(j,i)+co(j,konf(indexe+k))
91  enddo
92  coel(j,i)=coel(j,i)/nope
93  enddo
94  enddo
95 !
96  kflag=2
97 !
98 ! loop over all faces
99 !
100  do i=1,nface
101 !
102 ! check for cyclic symmetry
103 !
104  if(ifatie(i).ne.0) then
105  ics=abs(ifatie(i))
106  itie=int(cs(17,ics))
107  if(tieset(1,itie)(81:81).eq.'P') then
108  if(ifirst_occurrence.eq.1) then
109  do k=1,3
110  diff(k)=-cs(5+k,ics)
111  enddo
112  ifirst_occurrence=0
113  endif
114  elseif(tieset(1,itie)(81:81).eq.'Z') then
115  if(ifirst_occurrence.eq.1) then
116  icyclic=1
117  pi=4.d0*datan(1.d0)
118 !
119 ! normal along the cyclic symmetry axis such that
120 ! the slave surface is rotated clockwise through the
121 ! body into the master surface while looking in
122 ! the direction of xn
123 !
124  do k=1,3
125  a(k)=cs(5+k,ics)
126  xn(k)=cs(8+k,ics)-a(k)
127  enddo
128  dd=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3))
129  do k=1,3
130  xn(k)=xn(k)/dd
131  enddo
132 !
133 ! angle from the master to the slave surface
134 !
135  theta=-2.d0*pi/cs(1,ics)
136 !
137 ! rotation matrix rotating a vector in the master
138 ! surface into a vector in the slave surface
139 !
140  dc=dcos(theta)
141  ds=dsin(theta)
142 !
143 ! C-matrix from Guido Dhondt, The Finite Element
144 ! Method for Three-Dimensional Thermomechanical
145 ! Applications p 158
146 !
147  c(1,1)=dc+(1.d0-dc)*xn(1)*xn(1)
148  c(1,2)= (1.d0-dc)*xn(1)*xn(2)-ds*xn(3)
149  c(1,3)= (1.d0-dc)*xn(1)*xn(3)+ds*xn(2)
150  c(2,1)= (1.d0-dc)*xn(2)*xn(1)+ds*xn(3)
151  c(2,2)=dc+(1.d0-dc)*xn(2)*xn(2)
152  c(2,3)= (1.d0-dc)*xn(2)*xn(3)-ds*xn(1)
153  c(3,1)= (1.d0-dc)*xn(3)*xn(1)-ds*xn(2)
154  c(3,2)= (1.d0-dc)*xn(3)*xn(2)+ds*xn(1)
155  c(3,3)=dc+(1.d0-dc)*xn(3)*xn(3)
156  ifirst_occurrence=0
157  endif
158  else
159  write(*,*) '*ERROR in initialcfd'
160  write(*,*) ' kind of cyclic symmetry'
161  write(*,*) ' not known'
162  stop
163  endif
164  endif
165 !
166  iel1=ielfa(1,i)
167  indexe=ipkonf(iel1)
168  j1=ielfa(4,i)
169  if(lakonf(iel1)(4:4).eq.'8') then
170 !
171 ! hexahedral element
172 !
173 ! coordinates of the face centers
174 !
175  do j=1,4
176  nodes(j)=konf(indexe+ifaceq(j,j1))
177  do k=1,3
178  xl2(k,j)=co(k,nodes(j))
179  cofa(k,i)=cofa(k,i)+xl2(k,j)
180  enddo
181  enddo
182  do k=1,3
183  cofa(k,i)=cofa(k,i)/4.d0
184  enddo
185 !
186  xi=0.d0
187  et=0.d0
188  call shape4q(xi,et,xl2,xsj2,xs2,shp2,kflag)
189 !
190 ! area of the face
191 !
192  dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
193  & xsj2(3)*xsj2(3))
194  area(i)=4.d0*dxsj2
195 !
196  neighface=6
197 !
198  else if(lakonf(iel1)(4:4).eq.'6') then
199 !
200 ! wedge element
201 !
202 ! coordinates of the face centers
203 !
204  do j=1,nf(j1)
205  nodes(j)=konf(indexe+ifacew(j,j1))
206  do k=1,3
207  xl2(k,j)=co(k,nodes(j))
208  cofa(k,i)=cofa(k,i)+xl2(k,j)
209  enddo
210  enddo
211  do k=1,3
212  cofa(k,i)=cofa(k,i)/nf(j1)
213  enddo
214 !
215  xi=0.d0
216  et=0.d0
217  if(nf(j1).eq.3) then
218  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,kflag)
219  else
220  call shape4q(xi,et,xl2,xsj2,xs2,shp2,kflag)
221  endif
222 !
223 ! area of the face
224 !
225  dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
226  & xsj2(3)*xsj2(3))
227  if(nf(j1).eq.3) then
228  area(i)=dxsj2/2.d0
229  else
230  area(i)=4.d0*dxsj2
231  endif
232 !
233  neighface=5
234 !
235  else
236 !
237 ! tetrahedral element
238 !
239 ! coordinates of the face centers
240 !
241  do j=1,3
242  nodes(j)=konf(indexe+ifacet(j,j1))
243  do k=1,3
244  xl2(k,j)=co(k,nodes(j))
245  cofa(k,i)=cofa(k,i)+xl2(k,j)
246  enddo
247  enddo
248  do k=1,3
249  cofa(k,i)=cofa(k,i)/3.d0
250  enddo
251 !
252  xi=0.d0
253  et=0.d0
254  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,kflag)
255 !
256 ! area of the face
257 !
258  dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
259  & xsj2(3)*xsj2(3))
260  area(i)=dxsj2/2.d0
261 !
262  neighface=4
263 !
264  endif
265 !
266 ! normal and xi-vector on face viewed from cell 1
267 !
268  index1=ipnei(iel1)+j1
269  do k=1,3
270  xxn(k,index1)=xsj2(k)/dxsj2
271  xxi(k,index1)=cofa(k,i)-coel(k,iel1)
272  rf(k,i)=xxi(k,index1)
273  enddo
274 !
275 ! distance from face center to the center of cell 1
276 !
277  xle(index1)=dsqrt(xxi(1,index1)**2+xxi(2,index1)**2+
278  & xxi(3,index1)**2)
279  do k=1,3
280  xxi(k,index1)=xxi(k,index1)/xle(index1)
281  enddo
282 !
283 ! angle between the normal and the xi-vector
284 !
285  cosa(index1)=xxn(1,index1)*xxi(1,index1)+
286  & xxn(2,index1)*xxi(2,index1)+
287  & xxn(3,index1)*xxi(3,index1)
288 !
289  iel2=ielfa(2,i)
290 !
291 ! check whether there is an adjacent cell
292 !
293  if(iel2.ne.0) then
294  index2=ipnei(iel2)+neij(index1)
295 !
296 ! normal and xi-vector on face viewed from cell 2
297 !
298  if(ifatie(i).eq.0) then
299 !
300 ! genuine neighbor
301 !
302  do k=1,3
303  xxi(k,index2)=cofa(k,i)-coel(k,iel2)
304  xxn(k,index2)=-xxn(k,index1)
305  enddo
306 !
307  xle(index2)=dsqrt(xxi(1,index2)**2+xxi(2,index2)**2+
308  & xxi(3,index2)**2)
309  do k=1,3
310  xxi(k,index2)=xxi(k,index2)/xle(index2)
311  enddo
312 !
313 ! angle between the normal and the xi-vector: xxn.xxi
314 !
315  cosa(index2)=xxn(1,index2)*xxi(1,index2)+
316  & xxn(2,index2)*xxi(2,index2)+
317  & xxn(3,index2)*xxi(3,index2)
318 !
319 ! distance from the face center to the center of the
320 ! adjacent cell
321 !
322  xlen(index1)=xle(index2)
323  xlen(index2)=xle(index1)
324 !
325  do k=1,3
326  xxj(k,index2)=coel(k,iel1)-coel(k,iel2)
327  enddo
328 !
329 ! distance between the cell center and the center of the
330 ! adjacent cell
331 !
332  xlet(index1)=dsqrt(xxj(1,index2)**2+xxj(2,index2)**2
333  & +xxj(3,index2)**2)
334  xlet(index2)=xlet(index1)
335 !
336 ! xxj is the unit vector connecting neighboring cell centers
337 !
338  do k=1,3
339  xxj(k,index2)=xxj(k,index2)/xlet(index2)
340  xxj(k,index1)=-xxj(k,index2)
341  enddo
342 !
343 ! xxn.xxj
344 !
345  cosb(index2)=xxn(1,index1)*xxj(1,index1)+
346  & xxn(2,index1)*xxj(2,index1)+
347  & xxn(3,index1)*xxj(3,index1)
348  cosb(index1)=cosb(index2)
349 !
350 c xrlfa(1,i)=xle(index2)/(xle(index1)+xle(index2))
351 c xrlfa(2,i)=xle(index1)/(xle(index1)+xle(index2))
352  else
353 !
354 ! cyclic symmetry face: some quantities are
355 ! calculated on the cyclic symmetric face
356 !
357  xlen(index2)=xle(index1)
358 !
359 ! rotational cyclic symmetry
360 !
361 ! vector from the axis to the center of the
362 ! cyclic symmetric cell and orthogonal to the
363 ! axis
364 !
365  if(tieset(1,itie)(81:81).eq.'Z') then
366  do k=1,3
367  p(k)=coel(k,iel2)-a(k)
368  enddo
369  dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3)
370  do k=1,3
371  p(k)=p(k)-dd*xn(k)
372  enddo
373 !
374  if(ifatie(i).gt.0) then
375 !
376 ! vector rotated in the direction of the slave surface
377 ! (iel2 is adjacent to the master surface)
378 !
379  do k=1,3
380  q(k)=c(k,1)*p(1)+c(k,2)*p(2)+c(k,3)*p(3)
381  enddo
382  else
383 !
384 ! vector rotated in the direction of the master surface
385 ! (iel2 is adjacent to the slave surface)
386 !
387  do k=1,3
388  q(k)=c(1,k)*p(1)+c(2,k)*p(2)+c(3,k)*p(3)
389  enddo
390  endif
391 !
392 ! vector connecting the center of the cyclic
393 ! symmetry cell with the center of its rotated
394 ! ghost cell
395 !
396  do k=1,3
397  diff(k)=q(k)-p(k)
398  enddo
399  endif
400 !
401  do k=1,3
402  xxj(k,index1)=coel(k,iel2)-coel(k,iel1)
403  & +diff(k)
404  enddo
405 !
406 ! distance between the cell center and the center of the
407 ! adjacent cell
408 !
409  xlet(index1)=dsqrt(xxj(1,index1)**2+xxj(2,index1)**2
410  & +xxj(3,index1)**2)
411 !
412 ! xxj is the unit vector connecting neighboring cell centers
413 !
414  do k=1,3
415  xxj(k,index1)=xxj(k,index1)/xlet(index1)
416  enddo
417 !
418 ! xxn.xxj
419 !
420  cosb(index1)=xxn(1,index1)*xxj(1,index1)+
421  & xxn(2,index1)*xxj(2,index1)+
422  & xxn(3,index1)*xxj(3,index1)
423  endif
424 !
425 ! calculating rf = the shortest vector between the
426 ! center of the face and the line connecting
427 ! the centers of the adjacent elements. The
428 ! corresponding point on this line is called p
429 ! calculating xrlfa = the ratio of the segments created by
430 ! the point p to the total distance between
431 ! the centers of the adjacent elements of
432 ! the face
433 !
434  dd=rf(1,i)*xxj(1,index1)+
435  & rf(2,i)*xxj(2,index1)+
436  & rf(3,i)*xxj(3,index1)
437 !
438  do k=1,3
439  rf(k,i)=rf(k,i)-dd*xxj(k,index1)
440  enddo
441 !
442  xrlfa(2,i)=dd/xlet(index1)
443  xrlfa(1,i)=1.d0-xrlfa(2,i)
444  else
445 !
446 ! xxi and xxj coincide
447 !
448  do k=1,3
449  xxj(k,index1)=xxi(k,index1)
450  enddo
451  cosb(index1)=cosa(index1)
452 !
453 ! external face: determining the cell next to the
454 ! adjacent cell
455 !
456  iel3=ielfa(3,i)
457  if(iel3.eq.0) cycle
458 c xl13=dsqrt((coel(1,iel1)-coel(1,iel3))**2+
459 c & (coel(2,iel1)-coel(2,iel3))**2+
460 c & (coel(3,iel1)-coel(3,iel3))**2)
461 c xrlfa(1,i)=(xl13+xle(index1))/xl13
462 c xrlfa(3,i)=1.d0-xrlfa(1,i)
463 !
464 ! unit vector pointing from the center in element iel3
465 ! to the center in element iel1
466 !
467  do k=1,3
468  x13(k)=coel(k,iel1)-coel(k,iel3)
469  enddo
470 !
471  xl13=dsqrt(x13(1)*x13(1)+x13(2)*x13(2)+x13(3)*x13(3))
472 !
473  do k=1,3
474  x13(k)=x13(k)/xl13
475  enddo
476 !
477  dd=rf(1,i)*x13(1)+rf(2,i)*x13(2)+rf(3,i)*x13(3)
478 !
479  do k=1,3
480  rf(k,i)=rf(k,i)-dd*x13(k)
481  enddo
482 !
483  xrlfa(3,i)=-dd/xl13
484  xrlfa(1,i)=1.d0-xrlfa(3,i)
485 !
486  endif
487  enddo
488 !
489 ! for cyclic symmetric faces xrlfa has not been filled yet
490 !
491 c if(ifirst_occurrence.eq.0) then
492 c do i=1,nface
493 c if(ifatie(i).ne.0) then
494 c index1=ipnei(ielfa(1,i))+ielfa(4,i)
495 c xrlfa(1,i)=xlen(index1)/(xle(index1)+xlen(index1))
496 c xrlfa(2,i)=1.d0-xrlfa(1,i)
497 c endif
498 c enddo
499 c endif
500 !
501 ! calculation of the volume of the elements
502 !
503  do i=1,nef
504  if(ipkonf(i).lt.0) cycle
505  if(lakonf(i)(1:1).ne.'F') cycle
506  indexf=ipnei(i)
507  volume(i)=0.d0
508  do j=1,ipnei(i+1)-ipnei(i)
509  iface=neifa(indexf+j)
510  volume(i)=volume(i)+
511  & area(iface)*cofa(1,iface)*xxn(1,indexf+j)
512  enddo
513  enddo
514 !
515 ! calculation of the minimum length within the cells
516 !
517  dmin=1.d30
518  do i=1,nef
519  indexe=ipkonf(i)
520  read(lakonf(i)(4:4),'(i1)') nope
521  do j=1,nope
522  do k=1,3
523  xl(k,j)=co(k,konf(indexe+j))
524  enddo
525  enddo
526  if(nope.eq.4) then
527  do j=1,6
528  dmin=min(dmin,(xl(1,ied4(1,j))-xl(1,ied4(2,j)))**2+
529  & (xl(2,ied4(1,j))-xl(2,ied4(2,j)))**2+
530  & (xl(3,ied4(1,j))-xl(3,ied4(2,j)))**2)
531  enddo
532  elseif(nope.eq.6) then
533  do j=1,9
534  dmin=min(dmin,(xl(1,ied6(1,j))-xl(1,ied6(2,j)))**2+
535  & (xl(2,ied6(1,j))-xl(2,ied6(2,j)))**2+
536  & (xl(3,ied6(1,j))-xl(3,ied6(2,j)))**2)
537  enddo
538  else
539  do j=1,12
540  dmin=min(dmin,(xl(1,ied8(1,j))-xl(1,ied8(2,j)))**2+
541  & (xl(2,ied8(1,j))-xl(2,ied8(2,j)))**2+
542  & (xl(3,ied8(1,j))-xl(3,ied8(2,j)))**2)
543  enddo
544  endif
545  enddo
546  dmin=dsqrt(dmin)
547 !
548 ! calculate the distance to the nearest node for solid surface
549 ! faces
550 !
551  if(iturbulent.gt.0) then
552 !
553  if(dabs(physcon(5)).le.0.d0) then
554  write(*,*) '*ERROR in initialcfd: velocity at infinity'
555  write(*,*) ' is nonpositive;'
556  write(*,*) ' wrong *VALUES AT INFINITY'
557  call exit(201)
558  endif
559 !
560  if(dabs(physcon(7)).le.0.d0) then
561  write(*,*) '*ERROR in initialcfd: density at infinity'
562  write(*,*) ' is nonpositive;'
563  write(*,*) ' wrong *VALUES AT INFINITY'
564  call exit(201)
565  endif
566 !
567  if(dabs(physcon(8)).le.0.d0) then
568  write(*,*) '*ERROR in initialcfd: length of the '
569  write(*,*) ' computational domain is nonpositive;'
570  write(*,*) ' wrong *VALUES AT INFINITY'
571  call exit(201)
572  endif
573 !
574  do i=1,nsolidsurf
575  iface=isolidsurf(i)
576  nelem=int(iface/10)
577  indexe=ipkonf(nelem)
578  jface=iface-nelem*10
579 !
580 ! xl contains the coordinates of the nodes belonging
581 ! to the face
582 !
583  if(lakonf(nelem)(4:4).eq.'8') then
584  do j=1,4
585  node=konf(indexe+ifaceq(j,jface))
586  do k=1,3
587  xl(k,j)=co(k,node)
588  enddo
589  enddo
590  elseif(lakonf(nelem)(4:4).eq.'6') then
591  if(jface.le.2) then
592  do j=1,3
593  node=konf(indexe+ifacew(j,jface))
594  do k=1,3
595  xl(k,j)=co(k,node)
596  enddo
597  enddo
598  else
599  do j=1,4
600  node=konf(indexe+ifacew(j,jface))
601  do k=1,3
602  xl(k,j)=co(k,node)
603  enddo
604  enddo
605  endif
606  else
607  do j=1,3
608  node=konf(indexe+ifacet(j,jface))
609  do k=1,3
610  xl(k,j)=co(k,node)
611  enddo
612  enddo
613  endif
614 !
615 ! determine the plane through the face (exact for
616 ! 3-node face, approximate for a 4-node face)
617 !
618  if((lakonf(nelem)(4:4).eq.'8').or.
619  & ((lakonf(nelem)(4:4).eq.'6').and.(jface.gt.2))) then
620 !
621 ! computation of the local derivative of the global coordinates
622 ! (xs)
623 !
624  do j=1,3
625  xs(j,1)=-xl(j,1)+xl(j,2)+xl(j,3)-xl(j,4)
626  xs(j,2)=-xl(j,1)-xl(j,2)+xl(j,3)+xl(j,4)
627  enddo
628 !
629 ! computation of the jacobian vector for xi,et=0
630 !
631  aa=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2)
632  bb=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1)
633  cc=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2)
634  dd=dsqrt(aa*aa+bb*bb+cc*cc)
635  aa=aa/dd
636  bb=bb/dd
637  cc=cc/dd
638  dd=-(aa*(xl(1,1)+xl(1,2)+xl(1,3)+xl(1,4))
639  & +bb*(xl(2,1)+xl(2,2)+xl(2,3)+xl(2,4))
640  & +cc*(xl(3,1)+xl(3,2)+xl(3,3)+xl(3,4)))/4.d0
641  else
642 !
643 ! computation of the local derivative of the global coordinates
644 ! (xs)
645 !
646  do j=1,3
647  xs(j,1)=-xl(j,1)+xl(j,2)
648  xs(j,2)=-xl(j,1)+xl(j,3)
649  enddo
650 !
651 ! computation of the jacobian vector (unique for triangle)
652 !
653  aa=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2)
654  bb=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1)
655  cc=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2)
656  dd=dsqrt(aa*aa+bb*bb+cc*cc)
657  aa=aa/dd
658  bb=bb/dd
659  cc=cc/dd
660  dd=-(aa*xl(1,1)+bb*xl(2,1)+cc*xl(3,1))
661  endif
662 !
663 ! determine the shortest distance within the element
664 ! from the solid surface face
665 !
666  dist=1.d30
667  if(lakonf(nelem)(4:4).eq.'8') then
668  do j=1,4
669  node=konf(indexe+iopp8(j,jface))
670  dist=min(dist,-(aa*co(1,node)+bb*co(2,node)
671  & +cc*co(3,node)+dd))
672  enddo
673  elseif(lakonf(nelem)(4:4).eq.'6') then
674  if(jface.le.2) then
675  do j=1,3
676  node=konf(indexe+iopp6(j,jface))
677  dist=min(dist,-(aa*co(1,node)+bb*co(2,node)
678  & +cc*co(3,node)+dd))
679  enddo
680  else
681  do j=1,2
682  node=konf(indexe+iopp6(j,jface))
683  dist=min(dist,-(aa*co(1,node)+bb*co(2,node)
684  & +cc*co(3,node)+dd))
685  enddo
686  endif
687  else
688  node=konf(indexe+iopp8(1,jface))
689  dist=min(dist,-(aa*co(1,node)+bb*co(2,node)
690  & +cc*co(3,node)+dd))
691  endif
692 !
693 ! 60.d0/(0.075*delta(y)**2)
694 !
695  dy(i)=800.d0/(dist*dist)
696  enddo
697  endif
698 !
699  do i=1,nflnei
700  do k=1,3
701  xxni(k,i)=xxn(k,i)-xxi(k,i)
702  xxnj(k,i)=xxn(k,i)-xxj(k,i)
703  xxicn(k,i)=xxi(k,i)-cosa(i)*xxn(k,i)
704  enddo
705  enddo
706 c write(*,*) 'initialcfd neifa,neiel'
707 c do i=1,6*nef
708 c write(*,*) (i-1)/6+1,i-6*((i-1)/6),neifa(i),neiel(i)
709 c enddo
710 c write(*,*) 'initialcfd xle,xlen,xlet'
711 c do i=1,6*nef
712 c write(*,*) (i-1)/6+1,i-6*((i-1)/6),xle(i),xlen(i),xlet(i)
713 c enddo
714 c write(*,*) 'initialcfd xxn'
715 c do i=1,6*nef
716 c write(*,*) (i-1)/6+1,i-6*((i-1)/6),(xxn(j,i),j=1,3)
717 c enddo
718 c write(*,*) 'initialcfd xxi'
719 c do i=1,6*nef
720 c write(*,*) (i-1)/6+1,i-6*((i-1)/6),(xxi(j,i),j=1,3)
721 c enddo
722 c write(*,*) 'initialcfd xxj'
723 c do i=1,6*nef
724 c write(*,*) (i-1)/6+1,i-6*((i-1)/6),(xxj(j,i),j=1,3)
725 c enddo
726 c write(*,*) 'initialcfd cosa,cosb'
727 c do i=1,6*nef
728 c write(*,*) (i-1)/6+1,i-6*((i-1)/6),cosa(i),cosb(i)
729 c enddo
730 !
731  return
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
min
#define min(a, b)
Definition: cascade.c:31
stop
subroutine stop()
Definition: stop.f:20
dist
static double * dist
Definition: radflowload.c:42
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
nodes
subroutine nodes(inpc, textpart, co, nk, nk_, set, istartset, iendset, ialset, nset, nset_, nalset, nalset_, istep, istat, n, iline, ipol, inl, ipoinp, inp, ipoinpc)
Definition: nodes.f:22
Hosted by OpenAircraft.com, (Michigan UAV, LLC)