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

Go to the source code of this file.

Functions/Subroutines

subroutine precfd (ne, ipkon, kon, lakon, ipnei, neifa, neiel, ipoface, nodface, ielfa, nflnei, nface, ifaext, nfaext, isolidsurf, nsolidsurf, set, nset, istartset, iendset, ialset, vel, vold, mi, neij, nef, nactdoh, ipkonf, lakonf, ielmatf, ielmat, ielorienf, ielorien, norien, cs, mcs, tieset, x, y, z, xo, yo, zo, nx, ny, nz, co, ifatie)
 

Function/Subroutine Documentation

◆ precfd()

subroutine precfd ( integer  ne,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
integer, dimension(*)  ipnei,
integer, dimension(*)  neifa,
integer, dimension(*)  neiel,
integer, dimension(*)  ipoface,
integer, dimension(5,*)  nodface,
integer, dimension(4,*)  ielfa,
integer  nflnei,
integer  nface,
integer, dimension(*)  ifaext,
integer  nfaext,
integer, dimension(*)  isolidsurf,
integer  nsolidsurf,
character*81, dimension(*)  set,
integer  nset,
integer, dimension(*)  istartset,
integer, dimension(*)  iendset,
integer, dimension(*)  ialset,
real*8, dimension(nef,0:7)  vel,
real*8, dimension(0:mi(2),*)  vold,
integer, dimension(*)  mi,
integer, dimension(*)  neij,
integer  nef,
integer, dimension(*)  nactdoh,
integer, dimension(*)  ipkonf,
character*8, dimension(*)  lakonf,
integer, dimension(mi(3),*)  ielmatf,
integer, dimension(mi(3),*)  ielmat,
integer, dimension(mi(3),*)  ielorienf,
integer, dimension(mi(3),*)  ielorien,
integer  norien,
real*8, dimension(17,*)  cs,
integer  mcs,
character*81, dimension(3,*)  tieset,
real*8, dimension(*)  x,
real*8, dimension(*)  y,
real*8, dimension(*)  z,
real*8, dimension(*)  xo,
real*8, dimension(*)  yo,
real*8, dimension(*)  zo,
integer, dimension(*)  nx,
integer, dimension(*)  ny,
integer, dimension(*)  nz,
real*8, dimension(3,*)  co,
integer, dimension(*)  ifatie 
)
25 !
26 ! gathering topological information for computational fluid
27 ! dynamics calculations
28 !
29  implicit none
30 !
31  character*8 lakon(*),lakonf(*)
32  character*81 set(*),noset,tieset(3,*),slavset,mastset
33 !
34  integer ne,ipkon(*),ipnei(*),ipoface(*),nodface(5,*),neifa(*),
35  & ielfa(4,*),nflnei,nface,i,j,k,index,indexe,neiel(*),ithree,
36  & nfaext,ifaext(*),isolidsurf(*),nsolidsurf,indexf,nneigh,
37  & nset,istartset(*),iendset(*),ialset(*),iaux,kflag,ifour,iel2,
38  & ifaceq(8,6),ifacet(7,4),ifacew(8,5),kon(*),nodes(4),iel1,j2,
39  & indexold,ifree,ifreenew,ifreenei,mi(*),neij(*),ifreenei2,nef,
40  & nactdoh(*),ipkonf(*),ielmatf(mi(3),*),ielmat(mi(3),*),nf(5),
41  & nope,ielorien(mi(3),*),ielorienf(mi(3),*),norien,jopposite8(6),
42  & jopposite6(5),itie,nx(*),ny(*),nz(*),noden(1),nelemm,nelems,
43  & n,mcs,l,jfacem,jfaces,islav,imast,ifaces,ifacem,ifatie(*),
44  & nodeinface,nodeoutface,nopes,jop
45 !
46  real*8 vel(nef,0:7),vold(0:mi(2),*),coords(3),cs(17,*),x(*),
47  & y(*),z(*),xo(*),yo(*),zo(*),co(3,*),a(3),b(3),xn(3),p(3),
48  & q(3),c(3,3),dot,dc,ds,dd,theta,pi
49 !
50 ! nodes belonging to the element faces
51 !
52  data ifaceq /4,3,2,1,11,10,9,12,
53  & 5,6,7,8,13,14,15,16,
54  & 1,2,6,5,9,18,13,17,
55  & 2,3,7,6,10,19,14,18,
56  & 3,4,8,7,11,20,15,19,
57  & 4,1,5,8,12,17,16,20/
58  data ifacet /1,3,2,7,6,5,11,
59  & 1,2,4,5,9,8,12,
60  & 2,3,4,6,10,9,13,
61  & 1,4,3,8,10,7,14/
62  data ifacew /1,3,2,9,8,7,0,0,
63  & 4,5,6,10,11,12,0,0,
64  & 1,2,5,4,7,14,10,13,
65  & 2,3,6,5,8,15,11,14,
66  & 4,6,3,1,12,15,9,13/
67  data jopposite6 /2,1,0,0,0/
68  data jopposite8 /2,1,5,6,3,4/
69  data nf /3,3,4,4,4/
70 !
71 ! transfer element information from ipkon, lakon, ielmat
72 ! and ielorien to ipkonf, lakonf, ielmatf and ielorienf.
73 ! The latter fields were created
74 ! for CFD applications in which the elements are numerated in
75 ! ascending order (without gaps)
76 !
77  do i=1,ne
78 c write(*,*) 'precfd ',i,nactdoh(i)
79  if(nactdoh(i).ne.0) then
80  ipkonf(nactdoh(i))=ipkon(i)
81  lakonf(nactdoh(i))=lakon(i)
82  do j=1,mi(3)
83  ielmatf(j,nactdoh(i))=ielmat(j,i)
84  if(norien.gt.0) ielorienf(j,nactdoh(i))=ielorien(j,i)
85  enddo
86  endif
87  enddo
88 !
89  kflag=1
90  ithree=3
91  ifour=4
92 !
93 ! determining the external element faces of the fluid mesh
94 ! the faces are catalogued by the three lowes nodes numbers
95 ! in ascending order. ipoface(i) points to a face for which
96 ! node i is the lowest node and nodface(1,ipoface(i)) and
97 ! nodface(2,ipoface(i)) are the next lower ones.
98 ! nodface(3,ipoface(i)) contains the element number,
99 ! nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i))
100 ! is a pointer to the next surface for which node i is the
101 ! lowest node; if there are no more such surfaces the pointer
102 ! has the value zero
103 ! An external element face is one which belongs to one element
104 ! only
105 !
106  ifree=1
107  ifreenei=0
108 !
109  do i=1,6*nef
110  nodface(5,i)=i+1
111  enddo
112  do i=1,nef
113  indexe=ipkonf(i)
114  if(lakonf(i)(4:4).eq.'8') then
115 !
116 ! hex element
117 !
118  ipnei(i)=ifreenei
119  do j=1,6
120  do k=1,4
121  nodes(k)=kon(indexe+ifaceq(k,j))
122  enddo
123  call isortii(nodes,iaux,ifour,kflag)
124  indexold=0
125  index=ipoface(nodes(1))
126  do
127 !
128 ! adding a surface which has not been
129 ! catalogued so far
130 !
131  if(index.eq.0) then
132  ifreenew=nodface(5,ifree)
133  nodface(1,ifree)=nodes(2)
134  nodface(2,ifree)=nodes(3)
135  nodface(3,ifree)=i
136  nodface(4,ifree)=j
137  nodface(5,ifree)=ipoface(nodes(1))
138  ipoface(nodes(1))=ifree
139  ifreenei=ifreenei+1
140  neifa(ifreenei)=ifree
141  ielfa(1,ifree)=i
142  ielfa(4,ifree)=j
143  ifree=ifreenew
144  exit
145  endif
146 !
147 ! removing a surface which has already
148 ! been catalogued
149 !
150  if((nodface(1,index).eq.nodes(2)).and.
151  & (nodface(2,index).eq.nodes(3))) then
152  ifreenei=ifreenei+1
153 !
154 ! completing the facial info in neifa
155 !
156  neifa(ifreenei)=index
157  ielfa(2,index)=i
158 !
159 ! the neighboring elements to the face are i and iel2
160 ! with local face number for the face of j and j2
161 !
162  iel2=ielfa(1,index)
163  j2=ielfa(4,index)
164 !
165 ! completing the neighboring info for (element i,side j)
166 !
167  neiel(ifreenei)=iel2
168  neij(ifreenei)=j2
169 !
170 ! completing the neighboring info for (element iel2,side j2)
171 !
172  ifreenei2=ipnei(iel2)+j2
173  neiel(ifreenei2)=i
174  neij(ifreenei2)=j
175  exit
176  endif
177  indexold=index
178  index=nodface(5,index)
179  enddo
180  enddo
181  else if(lakonf(i)(4:4).eq.'6') then
182 !
183 ! wedge element
184 !
185  ipnei(i)=ifreenei
186  do j=1,5
187  do k=1,nf(j)
188  nodes(k)=kon(indexe+ifacew(k,j))
189  enddo
190  call isortii(nodes,iaux,nf(j),kflag)
191  indexold=0
192  index=ipoface(nodes(1))
193  do
194 !
195 ! adding a surface which has not been
196 ! catalogued so far
197 !
198  if(index.eq.0) then
199  ifreenew=nodface(5,ifree)
200  nodface(1,ifree)=nodes(2)
201  nodface(2,ifree)=nodes(3)
202  nodface(3,ifree)=i
203  nodface(4,ifree)=j
204  nodface(5,ifree)=ipoface(nodes(1))
205  ipoface(nodes(1))=ifree
206  ifreenei=ifreenei+1
207  neifa(ifreenei)=ifree
208  ielfa(1,ifree)=i
209  ielfa(4,ifree)=j
210  ifree=ifreenew
211  exit
212  endif
213 !
214 ! removing a surface which has already
215 ! been catalogued
216 !
217  if((nodface(1,index).eq.nodes(2)).and.
218  & (nodface(2,index).eq.nodes(3))) then
219  ifreenei=ifreenei+1
220 !
221 ! completing the facial info in neifa
222 !
223  neifa(ifreenei)=index
224  ielfa(2,index)=i
225 !
226 ! the neighboring elements to the face are i and iel2
227 ! with local face number for the face of j and j2
228 !
229  iel2=ielfa(1,index)
230  j2=ielfa(4,index)
231 !
232 ! completing the neighboring info for (element i,side j)
233 !
234  neiel(ifreenei)=iel2
235  neij(ifreenei)=j2
236 !
237 ! completing the neighboring info for (element iel2,side j2)
238 !
239  ifreenei2=ipnei(iel2)+j2
240  neiel(ifreenei2)=i
241  neij(ifreenei2)=j
242  exit
243  endif
244  indexold=index
245  index=nodface(5,index)
246  enddo
247  enddo
248  else
249 !
250 ! tet element
251 !
252  ipnei(i)=ifreenei
253  do j=1,4
254  do k=1,3
255  nodes(k)=kon(indexe+ifacet(k,j))
256  enddo
257  call isortii(nodes,iaux,ithree,kflag)
258  indexold=0
259  index=ipoface(nodes(1))
260  do
261 !
262 ! adding a surface which has not been
263 ! catalogued so far
264 !
265  if(index.eq.0) then
266  ifreenew=nodface(5,ifree)
267  nodface(1,ifree)=nodes(2)
268  nodface(2,ifree)=nodes(3)
269  nodface(3,ifree)=i
270  nodface(4,ifree)=j
271  nodface(5,ifree)=ipoface(nodes(1))
272  ipoface(nodes(1))=ifree
273  ifreenei=ifreenei+1
274  neifa(ifreenei)=ifree
275  ielfa(1,ifree)=i
276  ielfa(4,ifree)=j
277  ifree=ifreenew
278  exit
279  endif
280 !
281 ! removing a surface which has already
282 ! been catalogued
283 !
284  if((nodface(1,index).eq.nodes(2)).and.
285  & (nodface(2,index).eq.nodes(3))) then
286  ifreenei=ifreenei+1
287 !
288 ! completing the facial info in neifa
289 !
290  neifa(ifreenei)=index
291  ielfa(2,index)=i
292 !
293 ! the neighboring elements to the face are i and iel2
294 ! with local face number for the face of j and j2
295 !
296  iel2=ielfa(1,index)
297  j2=ielfa(4,index)
298 !
299 ! completing the neighboring info for (element i,side j)
300 !
301  neiel(ifreenei)=iel2
302  neij(ifreenei)=j2
303 !
304 ! completing the neighboring info for (element iel2,side j2)
305 !
306  ifreenei2=ipnei(iel2)+j2
307  neiel(ifreenei2)=i
308  neij(ifreenei2)=j
309  exit
310  endif
311  indexold=index
312  index=nodface(5,index)
313  enddo
314  enddo
315  endif
316  enddo
317 !
318  nflnei=ifreenei
319  ipnei(nef+1)=ifreenei
320  nface=ifree-1
321 !
322 ! check for cyclic symmetry (translational or rotational)
323 !
324  do i=1,mcs
325  itie=int(cs(17,i))
326  if((tieset(1,itie)(81:81).ne.'P').and.
327  & (tieset(1,itie)(81:81).ne.'Z')) cycle
328 !
329  slavset=tieset(2,itie)
330  mastset=tieset(3,itie)
331 !
332  do j=1,nset
333  if(set(j).eq.slavset) exit
334  enddo
335  islav=j
336 !
337  do j=1,nset
338  if(set(j).eq.mastset) exit
339  enddo
340  imast=j
341 !
342 ! catalogueing the center of gravity of the master faces
343 !
344  n=0
345  do j=istartset(imast),iendset(imast)
346  ifacem=ialset(j)
347 !
348  nelemm=int(ifacem/10)
349  jfacem=ifacem-nelemm*10
350 !
351 ! taking the renumbering of the elements into account
352 !
353  nelemm=nactdoh(nelemm)
354 !
355  indexe=ipkonf(nelemm)
356 !
357 ! determination of the face center of the master face
358 !
359  if(lakonf(nelemm)(4:4).eq.'8') then
360  do l=1,3
361  coords(l)=0.d0
362  enddo
363  do k=1,4
364  nodes(k)=kon(indexe+ifaceq(k,jfacem))
365  do l=1,3
366  coords(l)=coords(l)+co(l,nodes(k))
367  enddo
368  enddo
369  do l=1,3
370  coords(l)=coords(l)/4.d0
371  enddo
372  elseif(lakonf(nelemm)(4:4).eq.'6') then
373  do l=1,3
374  coords(l)=0.d0
375  enddo
376  do k=1,nf(jfacem)
377  nodes(k)=kon(indexe+ifacew(k,jfacem))
378  do l=1,3
379  coords(l)=coords(l)+co(l,nodes(k))
380  enddo
381  enddo
382  do l=1,3
383  coords(l)=coords(l)/nf(jfacem)
384  enddo
385  else
386  do l=1,3
387  coords(l)=0.d0
388  enddo
389  do k=1,3
390  nodes(k)=kon(indexe+ifaceq(k,jfacem))
391  do l=1,3
392  coords(l)=coords(l)+co(l,nodes(k))
393  enddo
394  enddo
395  do l=1,3
396  coords(l)=coords(l)/3.d0
397  enddo
398  endif
399 !
400  if(j.eq.istartset(imast)) then
401  if(tieset(1,itie)(81:81).eq.'Z') then
402 !
403 ! check the correct oprientation of the axis vector
404 !
405  if(lakonf(nelemm)(4:4).eq.'8') then
406  nope=8
407  nopes=4
408  elseif(lakonf(nelemm)(4:4).eq.'6') then
409  nope=6
410  nopes=nf(jfacem)
411  else
412  nope=4
413  nopes=3
414  endif
415 !
416 ! two nodes on the axis
417 !
418  do k=1,3
419  a(k)=cs(5+k,i)
420  b(k)=cs(8+k,i)
421  enddo
422 !
423 ! looking for a master node not on the axis
424 !
425  do k=1,nopes
426  if((co(1,nodes(k))-a(1))**2+
427  & (co(2,nodes(k))-a(2))**2+
428  & (co(3,nodes(k))-a(3))**2.gt.1.d-20) exit
429  enddo
430  nodeinface=nodes(k)
431 !
432 ! looking for a node belonging to the same element
433 ! but not in the master face
434 !
435  loop:do k=1,nope
436  nodeoutface=kon(indexe+k)
437  do l=1,nopes
438  if(nodes(l).eq.nodeoutface) cycle loop
439  enddo
440  exit
441  enddo loop
442 !
443 ! normalized vector along the axis
444 !
445  dd=dsqrt((b(1)-a(1))**2+
446  & (b(2)-a(2))**2+
447  & (b(3)-a(3))**2)
448  do k=1,3
449  xn(k)=(b(k)-a(k))/dd
450  enddo
451 !
452 ! vector connecting a on the axis with
453 ! nodeinface and nodeoutface
454 !
455  do k=1,3
456  p(k)=co(k,nodeinface)-a(k)
457  q(k)=co(k,nodeoutface)-a(k)
458  enddo
459 !
460 ! n.(pxq) must be negative for the axis vector
461 ! to have the right orientation (= you turn from the
462 ! slave surface to the master surface through the
463 ! body in clockwise direction while looking in
464 ! the direction of xn)
465 !
466  dot=xn(1)*(p(2)*q(3)-q(2)*p(3))
467  & +xn(2)*(p(3)*q(1)-q(3)*p(1))
468  & +xn(3)*(p(1)*q(2)-q(1)*p(2))
469 !
470 ! revert the direction of the axis if not correct
471 !
472  if(dot.gt.0.d0) then
473  do k=1,3
474  b(k)=2.d0*a(k)-b(k)
475  cs(8+k,i)=b(k)
476  xn(k)=-xn(k)
477  enddo
478  endif
479 !
480 ! angle from the master to the slave surface
481 !
482  pi=4.d0*datan(1.d0)
483  theta=-2.d0*pi/cs(1,i)
484 !
485 ! rotation matrix rotating a vector in the master
486 ! surface into a vector in the slave surface
487 !
488  dc=dcos(theta)
489  ds=dsin(theta)
490 !
491 ! C-matrix from Guido Dhondt, The Finite Element
492 ! Method for Three-Dimensional Thermomechanical
493 ! Applications p 158
494 !
495  c(1,1)=dc+(1.d0-dc)*xn(1)*xn(1)
496  c(1,2)= (1.d0-dc)*xn(1)*xn(2)-ds*xn(3)
497  c(1,3)= (1.d0-dc)*xn(1)*xn(3)+ds*xn(2)
498  c(2,1)= (1.d0-dc)*xn(2)*xn(1)+ds*xn(3)
499  c(2,2)=dc+(1.d0-dc)*xn(2)*xn(2)
500  c(2,3)= (1.d0-dc)*xn(2)*xn(3)-ds*xn(1)
501  c(3,1)= (1.d0-dc)*xn(3)*xn(1)-ds*xn(2)
502  c(3,2)= (1.d0-dc)*xn(3)*xn(2)+ds*xn(1)
503  c(3,3)=dc+(1.d0-dc)*xn(3)*xn(3)
504  endif
505  endif
506 !
507 ! storing the face centers (translated from the master
508 ! face into the slave face) in the vectors x,y and z
509 !
510  n=n+1
511  if(tieset(1,itie)(81:81).eq.'P') then
512  x(n)=coords(1)-cs(6,i)
513  y(n)=coords(2)-cs(7,i)
514  z(n)=coords(3)-cs(8,i)
515  else
516 !
517 ! vector orthogonal to axis (in the master surface)
518 !
519  do k=1,3
520  p(k)=coords(k)-a(k)
521  enddo
522  dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3)
523  do k=1,3
524  p(k)=p(k)-dd*xn(k)
525  enddo
526 !
527 ! vector rotated into the slave surface
528 !
529  do k=1,3
530  q(k)=c(k,1)*p(1)+c(k,2)*p(2)+c(k,3)*p(3)
531  enddo
532 !
533 ! rotated node in the slave surface
534 !
535  x(n)=coords(1)+q(1)-p(1)
536  y(n)=coords(2)+q(2)-p(2)
537  z(n)=coords(3)+q(3)-p(3)
538  endif
539  enddo
540 !
541 ! preparing the fields for near3d
542 !
543  do j=1,n
544  nx(j)=j
545  ny(j)=j
546  nz(j)=j
547  xo(j)=x(j)
548  yo(j)=y(j)
549  zo(j)=z(j)
550  enddo
551  kflag=2
552  call dsort(x,nx,n,kflag)
553  call dsort(y,ny,n,kflag)
554  call dsort(z,nz,n,kflag)
555 !
556 ! loop over all slave faces
557 !
558  do j=istartset(islav),iendset(islav)
559  ifaces=ialset(j)
560 !
561  nelems=int(ifaces/10)
562  jfaces=ifaces-nelems*10
563 !
564 ! taking the renumbering of the elements into account
565 !
566  nelems=nactdoh(nelems)
567 !
568  indexe=ipkonf(nelems)
569 !
570 ! determination of the face center of the slave face
571 !
572  if(lakonf(nelems)(4:4).eq.'8') then
573  do l=1,3
574  coords(l)=0.d0
575  enddo
576  do k=1,4
577  nodes(k)=kon(indexe+ifaceq(k,jfaces))
578  do l=1,3
579  coords(l)=coords(l)+co(l,nodes(k))
580  enddo
581  enddo
582  do l=1,3
583  coords(l)=coords(l)/4.d0
584  enddo
585  elseif(lakonf(nelems)(4:4).eq.'6') then
586  do l=1,3
587  coords(l)=0.d0
588  enddo
589  do k=1,nf(jfaces)
590  nodes(k)=kon(indexe+ifacew(k,jfaces))
591  do l=1,3
592  coords(l)=coords(l)+co(l,nodes(k))
593  enddo
594  enddo
595  do l=1,3
596  coords(l)=coords(l)/nf(jfaces)
597  enddo
598  else
599  do l=1,3
600  coords(l)=0.d0
601  enddo
602  do k=1,3
603  nodes(k)=kon(indexe+ifaceq(k,jfaces))
604  do l=1,3
605  coords(l)=coords(l)+co(l,nodes(k))
606  enddo
607  enddo
608  do l=1,3
609  coords(l)=coords(l)/3.d0
610  enddo
611  endif
612 !
613  nneigh=1
614  call near3d(xo,yo,zo,x,y,z,nx,ny,nz,coords(1),
615  & coords(2),coords(3),n,noden,nneigh)
616 !
617  ifacem=ialset(istartset(imast)+noden(1)-1)
618 !
619  nelemm=int(ifacem/10)
620  jfacem=ifacem-nelemm*10
621 !
622 ! taking the renumbering of the elements into account
623 !
624  nelemm=nactdoh(nelemm)
625 !
626  ielfa(2,neifa(ipnei(nelems)+jfaces))=nelemm
627  ielfa(2,neifa(ipnei(nelemm)+jfacem))=nelems
628 !
629  neiel(ipnei(nelems)+jfaces)=nelemm
630  neiel(ipnei(nelemm)+jfacem)=nelems
631 !
632  neij(ipnei(nelems)+jfaces)=jfacem
633  neij(ipnei(nelemm)+jfacem)=jfaces
634 !
635 ! field indicating whether a face (slave: +, master: -)
636 ! is part of a cyclic symmetry condition i
637 !
638  ifatie(neifa(ipnei(nelems)+jfaces))=i
639  ifatie(neifa(ipnei(nelemm)+jfacem))=-i
640  enddo
641  enddo
642 !
643 ! catalogueing external faces
644 !
645  nfaext=0
646 !
647  do i=1,nface
648  if(ielfa(2,i).ne.0) cycle
649  nfaext=nfaext+1
650  ifaext(nfaext)=i
651  iel1=ielfa(1,i)
652  indexf=ipnei(iel1)
653  if(lakonf(iel1)(4:4).eq.'8') then
654 !
655 ! hexes
656 !
657  j=ielfa(4,i)
658  j=jopposite8(j)
659  elseif(lakonf(iel1)(4:4).eq.'6') then
660 !
661 ! wedges
662 !
663  j=ielfa(4,i)
664 c j=jopposite6(j)
665 c if(j.eq.0) cycle
666  jop=jopposite6(j)
667  if(jop.eq.0) then
668  jop=j+1
669  if(jop.gt.5) jop=3
670  if(neiel(indexf+jop).eq.0) then
671  jop=j-1
672  if(jop.lt.3) jop=5
673  if(neiel(indexf+jop).eq.0) cycle
674  endif
675  j=jop
676  endif
677  else
678  cycle
679  endif
680 c indexf=ipnei(iel1)
681  ielfa(3,i)=neiel(indexf+j)
682  enddo
683 !
684 c do i=1,nef
685 c write(*,*)'precfd neiel',i,ipnei(i),(neiel(ipnei(i)+j),j=1,6)
686 c write(*,*)'precfd neij',i,ipnei(i),(neij(ipnei(i)+j),j=1,6)
687 c write(*,*)'precfd neifa',i,ipnei(i),(neifa(ipnei(i)+j),j=1,6)
688 c enddo
689 c do k=1,nface
690 c write(*,*) 'precfd ielfa',k,(ielfa(j,k),j=1,4),ifatie(k)
691 c enddo
692 c do k=1,nfaext
693 c write(*,*) 'precfd ifaext',k,ifaext(k)
694 c enddo
695 c write(*,*)
696 !
697 ! faces belonging to solid surfaces
698 !
699  nsolidsurf=0
700  noset(1:13)='SOLIDSURFACET'
701  do i=1,nset
702  if(set(i)(1:13).eq.noset(1:13)) exit
703  enddo
704  if(i.gt.nset) then
705  write(*,*) '*WARNING in precfd: facial surface SOLID SURFACE '
706  write(*,*) ' has not been defined.'
707  write(*,*)
708  else
709 c nsolidsurf=0
710  do j=istartset(i),iendset(i)
711  nsolidsurf=nsolidsurf+1
712  isolidsurf(nsolidsurf)=ialset(j)
713  enddo
714  call isortii(isolidsurf,iaux,nsolidsurf,kflag)
715  endif
716 c do i=1,nsolidsurf
717 c write(*,*) 'precfd solidsurf ',i,isolidsurf(i)
718 c enddo
719 !
720 ! initial conditions: element values
721 ! vel was initialized at allocation
722 !
723  do i=1,nef
724  indexe=ipkonf(i)
725  if(lakonf(i)(4:4).eq.'8') then
726  nope=8
727  elseif(lakonf(i)(4:4).eq.'6') then
728  nope=6
729  else
730  nope=4
731  endif
732  do j=0,4
733  do k=1,nope
734  vel(i,j)=vel(i,j)+vold(j,kon(indexe+k))
735  enddo
736  vel(i,j)=vel(i,j)/nope
737  enddo
738  enddo
739 c do i=1,nef
740 c write(*,*) 'precfd vel ',i,(vel(j,i),j=1,3)
741 c enddo
742 !
743  return
near3d
subroutine near3d(xo, yo, zo, x, y, z, nx, ny, nz, xp, yp, zp, n, neighbor, k)
Definition: near3d.f:20
isortii
subroutine isortii(ix, iy, n, kflag)
Definition: isortii.f:6
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
dsort
subroutine dsort(dx, iy, n, kflag)
Definition: dsort.f:6
Hosted by OpenAircraft.com, (Michigan UAV, LLC)