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

Go to the source code of this file.

Functions/Subroutines

subroutine cyclicsymmetrymodels (inpc, textpart, set, istartset, iendset, ialset, nset, tieset, tietol, co, nk, ipompc, nodempc, coefmpc, nmpc, nmpc_, ikmpc, ilmpc, mpcfree, rcs, zcs, ics, nr, nz, rcs0, zcs0, ncs_, cs, labmpc, istep, istat, n, iline, ipol, inl, ipoinp, inp, ntie, mcs, lprev, ithermal, rcscg, rcs0cg, zcscg, zcs0cg, nrcg, nzcg, jcs, kontri, straight, ne, ipkon, kon, lakon, lcs, ifacetet, inodface, ipoinpc, maxsectors, trab, ntrans, ntrans_, jobnamec, vold, cfd, mi, iaxial)
 

Function/Subroutine Documentation

◆ cyclicsymmetrymodels()

subroutine cyclicsymmetrymodels ( character*1, dimension(*)  inpc,
character*132, dimension(16)  textpart,
character*81, dimension(*)  set,
integer, dimension(*)  istartset,
integer, dimension(*)  iendset,
integer, dimension(*)  ialset,
integer  nset,
character*81, dimension(3,*)  tieset,
real*8, dimension(3,*)  tietol,
real*8, dimension(3,*)  co,
integer  nk,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer  nmpc_,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
integer  mpcfree,
real*8, dimension(*)  rcs,
real*8, dimension(*)  zcs,
integer, dimension(*)  ics,
integer, dimension(*)  nr,
integer, dimension(*)  nz,
real*8, dimension(*)  rcs0,
real*8, dimension(*)  zcs0,
integer  ncs_,
real*8, dimension(17,*)  cs,
character*20, dimension(*)  labmpc,
integer  istep,
integer  istat,
integer  n,
integer  iline,
integer  ipol,
integer  inl,
integer, dimension(2,*)  ipoinp,
integer, dimension(3,*)  inp,
integer  ntie,
integer  mcs,
integer  lprev,
integer  ithermal,
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, dimension(*)  ifacetet,
integer, dimension(*)  inodface,
integer, dimension(0:*)  ipoinpc,
integer  maxsectors,
real*8, dimension(7,*)  trab,
integer  ntrans,
integer  ntrans_,
character*132, dimension(*)  jobnamec,
real*8, dimension(0:mi(2),*)  vold,
integer  cfd,
integer, dimension(*)  mi,
integer  iaxial 
)
28 !
29 ! reading the input deck: *CYCLIC SYMMETRY MODEL
30 !
31 ! several cyclic symmetry parts can be defined for one and the
32 ! same model; for each part there must be a *CYCLIC SYMMETRY MODEL
33 ! card
34 !
35 ! cs(1,mcs): # segments in 360 degrees
36 ! cs(2,mcs): minimum node diameter
37 ! cs(3,mcs): maximum node diameter
38 ! cs(4,mcs): # nodes on the independent side (for fluids: upper bound)
39 ! cs(5,mcs): # sectors to be plotted
40 ! cs(6..8,mcs): first point on cyclic symmetry axis
41 ! cs(9..11,mcs): second point on cylic symmetry axis; turning
42 ! the slave surface clockwise about the cyclic symmetry axis
43 ! while looking from the first point to the second point one
44 ! arrives at the master surface without leaving the body
45 ! cs(12,mcs): -1 (denotes a cylindrical coordinate system)
46 ! cs(13,mcs): number of the element set
47 ! cs(14,mcs): sum of previous independent nodes
48 ! cs(15,mcs): cos(angle); angle = 2*pi/cs(1,mcs)
49 ! cs(16,mcs): sin(angle)
50 ! cs(17,mcs): number of tie constraint
51 !
52 ! notice that in this routine ics, zcs and rcs start for 1 for
53 ! each *cyclic symmetry model card (look at the pointer in
54 ! the cyclicsymmetrymodels call in calinput.f)
55 !
56  implicit none
57 !
58  logical triangulation,calcangle,nodesonaxis,check,exist
59 !
60  character*1 inpc(*),depkind,indepkind
61  character*8 lakon(*)
62  character*20 labmpc(*)
63  character*80 tie
64  character*81 set(*),depset,indepset,tieset(3,*),elset
65  character*132 textpart(16),jobnamec(*)
66 !
67  integer istartset(*),iendset(*),ialset(*),ipompc(*),ifaces,
68  & nodempc(3,*),itiecyc,ntiecyc,iaxial,nopes,nelems,m,indexe,
69  & nset,istep,istat,n,key,i,j,k,nk,nmpc,nmpc_,mpcfree,ics(*),
70  & nr(*),nz(*),jdep,jindep,l,noded,ikmpc(*),ilmpc(*),lcs(*),
71  & kflag,node,ncsnodes,ncs_,iline,ipol,inl,ipoinp(2,*),nneigh,
72  & inp(3,*),itie,iset,ipos,mcs,lprev,ntie,ithermal,ncounter,
73  & nrcg(*),nzcg(*),jcs(*),kontri(3,*),ne,ipkon(*),kon(*),nodei,
74  & ifacetet(*),inodface(*),ipoinpc(0:*),maxsectors,id,jfaces,
75  & noden(2),ntrans,ntrans_,cfd,mi(*),ifaceq(8,6),ifacet(6,4),
76  & ifacew1(4,5),ifacew2(8,5),idof
77 !
78  real*8 tolloc,co(3,*),coefmpc(*),rcs(*),zcs(*),rcs0(*),zcs0(*),
79  & csab(7),xn,yn,zn,dd,xap,yap,zap,tietol(3,*),cs(17,*),xsectors,
80  & gsectors,x3,y3,z3,phi,rcscg(*),rcs0cg(*),zcscg(*),zcs0cg(*),
81  & straight(9,*),x1,y1,z1,x2,y2,z2,zp,rp,dist,trab(7,*),
82  & vold(0:mi(2),*),calculated_angle,user_angle
83 !
84 ! nodes per face for hex elements
85 !
86  data ifaceq /4,3,2,1,11,10,9,12,
87  & 5,6,7,8,13,14,15,16,
88  & 1,2,6,5,9,18,13,17,
89  & 2,3,7,6,10,19,14,18,
90  & 3,4,8,7,11,20,15,19,
91  & 4,1,5,8,12,17,16,20/
92 !
93 ! nodes per face for tet elements
94 !
95  data ifacet /1,3,2,7,6,5,
96  & 1,2,4,5,9,8,
97  & 2,3,4,6,10,9,
98  & 1,4,3,8,10,7/
99 !
100 ! nodes per face for linear wedge elements
101 !
102  data ifacew1 /1,3,2,0,
103  & 4,5,6,0,
104  & 1,2,5,4,
105  & 2,3,6,5,
106  & 3,1,4,6/
107 !
108 ! nodes per face for quadratic wedge elements
109 !
110  data ifacew2 /1,3,2,9,8,7,0,0,
111  & 4,5,6,10,11,12,0,0,
112  & 1,2,5,4,7,14,10,13,
113  & 2,3,6,5,8,15,11,14,
114  & 3,1,4,6,9,13,12,15/
115 !
116  if(istep.gt.0) then
117  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
118  write(*,*) ' *CYCLIC SYMMETRY MODEL should'
119  write(*,*) ' be placed before all step definitions'
120  call exit(201)
121  endif
122 !
123  check=.true.
124  gsectors=1
125  elset='
126  & '
127  tie='
128  & '
129  do i=2,n
130  if(textpart(i)(1:2).eq.'N=') then
131  read(textpart(i)(3:22),'(f20.0)',iostat=istat) xsectors
132  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
133  &"*CYCLIC SYMMETRY MODEL%")
134  elseif(textpart(i)(1:8).eq.'CHECK=NO') then
135  check=.false.
136  elseif(textpart(i)(1:7).eq.'NGRAPH=') then
137  read(textpart(i)(8:27),'(f20.0)',iostat=istat) gsectors
138  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
139  &"*CYCLIC SYMMETRY MODEL%")
140  elseif(textpart(i)(1:4).eq.'TIE=') then
141  read(textpart(i)(5:84),'(a80)',iostat=istat) tie
142  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
143  &"*CYCLIC SYMMETRY MODEL%")
144  elseif(textpart(i)(1:6).eq.'ELSET=') then
145  read(textpart(i)(7:86),'(a80)',iostat=istat) elset
146  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
147  &"*CYCLIC SYMMETRY MODEL%")
148  elset(81:81)=' '
149  ipos=index(elset,' ')
150  elset(ipos:ipos)='E'
151  else
152  write(*,*)
153  & '*WARNING reading *CYCLIC SYMMETRY MODEL:'
154  write(*,*) ' parameter not recognized:'
155  write(*,*) ' ',
156  & textpart(i)(1:index(textpart(i),' ')-1)
157  call inputwarning(inpc,ipoinpc,iline,
158  &"*CYCLIC SYMMETRY MODEL%")
159  endif
160  enddo
161 !
162  if(xsectors.le.0) then
163  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
164  write(*,*) ' the required parameter N'
165  write(*,*) ' is lacking on the *CYCLIC SYMMETRY MODEL'
166  write(*,*) ' keyword card or has a value <=0'
167  call exit(201)
168  endif
169  if(gsectors.lt.1) then
170  write(*,*) '*WARNING reading *CYCLIC SYMMETRY MODEL:'
171  write(*,*) ' cannot plot less than'
172  write(*,*) ' one sector: one sector will be plotted'
173  gsectors=1
174  endif
175  if(gsectors.gt.xsectors) then
176  write(*,*) '*WARNING reading *CYCLIC SYMMETRY MODEL:'
177  write(*,*) ' cannot plot more than'
178  write(*,*) ' ',xsectors,'sectors;',
179  & xsectors,' sectors will'
180  write(*,*) ' be plotted'
181  gsectors=xsectors
182  endif
183 !
184  maxsectors=max(maxsectors,int(xsectors+0.5d0))
185 !
186  mcs=mcs+1
187  cs(2,mcs)=-0.5
188  cs(3,mcs)=-0.5
189  cs(14,mcs)=lprev+0.5
190 !
191 ! determining the tie constraint
192 !
193  itie=0
194  ntiecyc=0
195  do i=1,ntie
196  if((tieset(1,i)(1:80).eq.tie).and.
197  & (tieset(1,i)(81:81).eq.'P')) then
198  itie=i
199 !
200 ! fluid periodicity (translational)
201 !
202  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
203  & ipoinp,inp,ipoinpc)
204 !
205 ! translation vector from the slave surface to the
206 ! master surface
207 !
208  do j=1,3
209  read(textpart(j)(1:20),'(f20.0)',iostat=istat)
210  & cs(5+j,mcs)
211  enddo
212  cs(17,mcs)=itie+0.5
213 !
214 ! counting the number of faces on the master side (should be the
215 ! same as on the slave side)
216 !
217  indepset=tieset(3,itie)
218  do j=1,nset
219  if(set(j).eq.indepset) exit
220  enddo
221 !
222 ! max 8 nodes per face
223 !
224  cs(4,mcs)=(iendset(j)-istartset(j)+1)*8+0.5d0
225 !
226  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
227  & ipoinp,inp,ipoinpc)
228 !
229  return
230  elseif((tieset(1,i)(1:80).eq.tie).and.
231  & (tieset(1,i)(81:81).eq.'Z')) then
232  itie=i
233 !
234 ! fluid periodicity (rotational)
235 !
236  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
237  & ipoinp,inp,ipoinpc)
238 !
239 ! translation vector from the slave surface to the
240 ! master surface
241 !
242  do j=1,6
243  read(textpart(j)(1:20),'(f20.0)',iostat=istat)
244  & cs(5+j,mcs)
245  enddo
246  cs(1,mcs)=xsectors
247  cs(17,mcs)=itie+0.5
248 !
249 ! counting the number of faces on the master side (should be the
250 ! same as on the slave side)
251 !
252  indepset=tieset(3,itie)
253  do j=1,nset
254  if(set(j).eq.indepset) exit
255  enddo
256 !
257 ! max 8 nodes per face
258 !
259  cs(4,mcs)=(iendset(j)-istartset(j)+1)*8+0.5d0
260 !
261  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
262  & ipoinp,inp,ipoinpc)
263 !
264  return
265  elseif((tieset(1,i)(1:80).eq.tie).and.
266  & (tieset(1,i)(81:81).ne.'C').and.
267  & (tieset(1,i)(81:81).ne.'T')) then
268  itie=i
269  exit
270  elseif((tieset(1,i)(81:81).ne.'C').and.
271  & (tieset(1,i)(81:81).ne.'T')) then
272  ntiecyc=ntiecyc+1
273  itiecyc=i
274  endif
275  enddo
276  if(itie.eq.0) then
277  if(ntiecyc.eq.1) then
278  itie=itiecyc
279  else
280  write(*,*)
281  & '*ERROR reading *CYCLIC SYMMETRY MODEL:'
282  write(*,*) ' tie constraint is nonexistent'
283  call inputerror(inpc,ipoinpc,iline,
284  &"*CYCLIC SYMMETRY MODEL%")
285  endif
286  endif
287 !
288  cs(1,mcs)=xsectors
289  cs(5,mcs)=gsectors+0.5
290  cs(17,mcs)=itie+0.5
291  depset=tieset(2,itie)
292  indepset=tieset(3,itie)
293  tolloc=tietol(1,itie)
294 !
295 ! determining the element set
296 !
297  iset=0
298  if(elset.eq.' ') then
299  write(*,*) '*INFO reading *CYCLIC SYMMETRY MODEL:'
300  write(*,*) ' no element set given'
301  call inputinfo(inpc,ipoinpc,iline,
302  &"*CYCLIC SYMMETRY MODEL%")
303  else
304  do i=1,nset
305  if(set(i).eq.elset) then
306  iset=i
307  exit
308  endif
309  enddo
310  if(iset.eq.0) then
311  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
312  write(*,*) ' element set does not'
313  write(*,*) ' exist; '
314  call inputerror(inpc,ipoinpc,iline,
315  &"*CYCLIC SYMMETRY MODEL%")
316  endif
317  endif
318  cs(13,mcs)=iset+0.5
319 !
320  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
321  & ipoinp,inp,ipoinpc)
322 !
323  if((istat.lt.0).or.(key.eq.1)) then
324  write(*,*)'*ERROR reading *CYCLIC SYMMETRY MODEL:'
325  write(*,*) ' definition of the cyclic'
326  write(*,*) ' symmetry model is not complete'
327  call exit(201)
328  endif
329 !
330  ntrans=ntrans+1
331  if(ntrans.gt.ntrans_) then
332  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
333  write(*,*) ' increase ntrans_'
334  call exit(201)
335  endif
336 !
337  do i=1,6
338  read(textpart(i)(1:20),'(f20.0)',iostat=istat) csab(i)
339  trab(i,ntrans)=csab(i)
340  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
341  &"*CYCLIC SYMMETRY MODEL%")
342  enddo
343 !
344 ! cyclic coordinate system
345 !
346  csab(7)=-1.d0
347 !
348 ! marker for cyclic symmetry axis
349 !
350  trab(7,ntrans)=2
351 !
352 ! check whether depset and indepset exist
353 ! determine the kind of set (nodal or facial)
354 !
355  ipos=index(depset,' ')
356  depkind='S'
357  depset(ipos:ipos)=depkind
358  do i=1,nset
359  if(set(i).eq.depset) exit
360  enddo
361  if(i.gt.nset) then
362  depkind='T'
363  depset(ipos:ipos)=depkind
364  do i=1,nset
365  if(set(i).eq.depset) exit
366  enddo
367  if(i.gt.nset) then
368  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
369  write(*,*) ' surface ',depset
370  write(*,*) ' has not yet been defined.'
371  call exit(201)
372  endif
373  endif
374  jdep=i
375 !
376  ipos=index(indepset,' ')
377  indepkind='S'
378  indepset(ipos:ipos)=indepkind
379  do i=1,nset
380  if(set(i).eq.indepset) exit
381  enddo
382  if(i.gt.nset) then
383  indepkind='T'
384  indepset(ipos:ipos)=indepkind
385  do i=1,nset
386  if(set(i).eq.indepset) exit
387  enddo
388  if(i.gt.nset) then
389  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
390  write(*,*) ' surface ',indepset
391  write(*,*) ' has not yet been defined.'
392  call exit(201)
393  endif
394  endif
395  jindep=i
396 !
397 ! unit vector along the rotation axis (xn,yn,zn)
398 !
399  xn=csab(4)-csab(1)
400  yn=csab(5)-csab(2)
401  zn=csab(6)-csab(3)
402  dd=dsqrt(xn*xn+yn*yn+zn*zn)
403  xn=xn/dd
404  yn=yn/dd
405  zn=zn/dd
406 !
407 ! defining the indepset as a 2-D data field (axes: r=radial
408 ! coordinate, z=axial coordinate): needed to allocate a node
409 ! of the depset to a node of the indepset for the cyclic
410 ! symmetry equations
411 !
412  l=0
413  do j=istartset(jindep),iendset(jindep)
414  if(ialset(j).gt.0) then
415  if(indepkind.eq.'T') then
416 !
417 ! facial independent surface
418 !
419  ifaces=ialset(j)
420  nelems=int(ifaces/10)
421  jfaces=ifaces - nelems*10
422  indexe=ipkon(nelems)
423 !
424  if(lakon(nelems)(4:5).eq.'20') then
425  nopes=8
426  elseif(lakon(nelems)(4:4).eq.'2') then
427  nopes=9
428  elseif(lakon(nelems)(4:4).eq.'8') then
429  nopes=4
430  elseif(lakon(nelems)(4:5).eq.'10') then
431  nopes=6
432  elseif(lakon(nelems)(4:4).eq.'4') then
433  nopes=3
434  endif
435 !
436  if(lakon(nelems)(4:4).eq.'6') then
437  if(jfaces.le.2) then
438  nopes=3
439  else
440  nopes=4
441  endif
442  endif
443  if(lakon(nelems)(4:5).eq.'15') then
444  if(jfaces.le.2) then
445  nopes=6
446  else
447  nopes=8
448  endif
449  endif
450  else
451  nopes=1
452  endif
453 !
454  do m=1,nopes
455  if(indepkind.eq.'T') then
456  if((lakon(nelems)(4:4).eq.'2').or.
457  & (lakon(nelems)(4:4).eq.'8')) then
458  node=kon(indexe+ifaceq(m,jfaces))
459  elseif((lakon(nelems)(4:4).eq.'4').or.
460  & (lakon(nelems)(4:5).eq.'10')) then
461  node=kon(indexe+ifacet(m,jfaces))
462  elseif(lakon(nelems)(4:4).eq.'6') then
463  node=kon(indexe+ifacew1(m,jfaces))
464  elseif(lakon(nelems)(4:5).eq.'15') then
465  node=kon(indexe+ifacew2(m,jfaces))
466  endif
467  call nident(ics,node,l,id)
468  exist=.false.
469  if(id.gt.0) then
470  if(ics(id).eq.node) then
471  exist=.true.
472  endif
473  endif
474  if(exist) cycle
475  l=l+1
476  if(lprev+l.gt.ncs_) then
477  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
478  write(*,*) ' increase ncs_'
479  call exit(201)
480  endif
481  do k=l,id+2,-1
482  ics(k)=ics(k-1)
483  zcs(k)=zcs(k-1)
484  rcs(k)=rcs(k-1)
485  enddo
486 !
487  xap=co(1,node)-csab(1)
488  yap=co(2,node)-csab(2)
489  zap=co(3,node)-csab(3)
490 !
491  ics(id+1)=node
492  zcs(id+1)=xap*xn+yap*yn+zap*zn
493  rcs(id+1)=dsqrt((xap-zcs(id+1)*xn)**2+
494  & (yap-zcs(id+1)*yn)**2+
495  & (zap-zcs(id+1)*zn)**2)
496  else
497  l=l+1
498  if(lprev+l.gt.ncs_) then
499  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
500  write(*,*) ' increase ncs_'
501  call exit(201)
502  endif
503  node =ialset(j)
504 !
505  xap=co(1,node)-csab(1)
506  yap=co(2,node)-csab(2)
507  zap=co(3,node)-csab(3)
508 !
509  ics(l)=node
510  zcs(l)=xap*xn+yap*yn+zap*zn
511  rcs(l)=dsqrt((xap-zcs(l)*xn)**2+
512  & (yap-zcs(l)*yn)**2+
513  & (zap-zcs(l)*zn)**2)
514  endif
515  enddo
516  else
517  k=ialset(j-2)
518  do
519  k=k-ialset(j)
520  if(k.ge.ialset(j-1)) exit
521  l=l+1
522  if(l.gt.ncs_) then
523  write(*,*) '*ERROR reading *CYCLIC SYMMETRY MODEL:'
524  write(*,*) ' increase ncs_'
525  call exit(201)
526  endif
527  node=k
528 !
529  xap=co(1,node)-csab(1)
530  yap=co(2,node)-csab(2)
531  zap=co(3,node)-csab(3)
532 !
533  ics(l)=node
534  zcs(l)=xap*xn+yap*yn+zap*zn
535  rcs(l)=dsqrt((xap-zcs(l)*xn)**2+
536  & (yap-zcs(l)*yn)**2+
537  & (zap-zcs(l)*zn)**2)
538  enddo
539  endif
540  enddo
541 !
542  ncsnodes=l
543 !
544 ! initialization of near2d
545 !
546  do i=1,ncsnodes
547  nr(i)=i
548  nz(i)=i
549  rcs0(i)=rcs(i)
550  zcs0(i)=zcs(i)
551  enddo
552  kflag=2
553  call dsort(rcs,nr,ncsnodes,kflag)
554  call dsort(zcs,nz,ncsnodes,kflag)
555 !
556 ! check whether a tolerance was defined. If not, a tolerance
557 ! is calculated as 0.5 % of the mean of the distance of every
558 ! independent node to its nearest neighbour
559 !
560  if(tolloc.lt.0.d0) then
561  nneigh=2
562  dist=0.d0
563  do i=1,ncsnodes
564  nodei=ics(i)
565 !
566  xap=co(1,nodei)-csab(1)
567  yap=co(2,nodei)-csab(2)
568  zap=co(3,nodei)-csab(3)
569 !
570  zp=xap*xn+yap*yn+zap*zn
571  rp=dsqrt((xap-zp*xn)**2+(yap-zp*yn)**2+(zap-zp*zn)**2)
572 !
573  call near2d(rcs0,zcs0,rcs,zcs,nr,nz,rp,zp,ncsnodes,noden,
574  & nneigh)
575 !
576  dist=dist+dsqrt((co(1,nodei)-co(1,noden(2)))**2+
577  & (co(2,nodei)-co(2,noden(2)))**2+
578  & (co(3,nodei)-co(3,noden(2)))**2)
579  enddo
580  tolloc=1.d-10*dist/ncsnodes
581  write(*,*) '*INFO reading *CYCLIC SYMMETRY MODEL:'
582  write(*,*) ' no tolerance was defined'
583  write(*,*) ' in the *TIE option; a tolerance of ',
584  & tolloc
585  write(*,*) ' will be used'
586  write(*,*)
587  endif
588 !
589 ! calculating the angle between dependent and independent
590 ! side and check for nodes on the axis
591 !
592 ! this angle may be different from 2*pi/xsectors: in that way
593 ! the user can simulate fractional nodal diameters
594 !
595 ! (x2,y2,z2): unit vector on the dependent side and orthogonal
596 ! to the rotation axis
597 ! (x3,y3,z3): unit vector on the independent side and orthogonal
598 ! to the rotation axis
599 ! (x1,y1,z1)=(x2,y2,z2)x(x3,y3,z3)
600 ! points in the same direction of xn if the independent
601 ! side is on the clockwise side of the dependent side if
602 ! looking in the direction of xn
603 !
604  calcangle=.false.
605  nodesonaxis=.false.
606 !
607  nneigh=1
608  loop1: do i=istartset(jdep),iendset(jdep)
609  if(ialset(i).gt.0) then
610 !
611 ! check whether dependent side is node based or
612 ! face based
613 !
614  if(depkind.eq.'T') then
615  ifaces=ialset(i)
616  nelems=int(ifaces/10)
617  jfaces=ifaces - nelems*10
618  indexe=ipkon(nelems)
619 !
620  if(lakon(nelems)(4:5).eq.'20') then
621  nopes=8
622  elseif(lakon(nelems)(4:4).eq.'2') then
623  nopes=9
624  elseif(lakon(nelems)(4:4).eq.'8') then
625  nopes=4
626  elseif(lakon(nelems)(4:5).eq.'10') then
627  nopes=6
628  elseif(lakon(nelems)(4:4).eq.'4') then
629  nopes=3
630  endif
631 !
632  if(lakon(nelems)(4:4).eq.'6') then
633  if(jfaces.le.2) then
634  nopes=3
635  else
636  nopes=4
637  endif
638  endif
639  if(lakon(nelems)(4:5).eq.'15') then
640  if(jfaces.le.2) then
641  nopes=6
642  else
643  nopes=8
644  endif
645  endif
646  else
647  nopes=1
648  endif
649 !
650  do m=1,nopes
651  if(depkind.eq.'T') then
652  if((lakon(nelems)(4:4).eq.'2').or.
653  & (lakon(nelems)(4:4).eq.'8')) then
654  noded=kon(indexe+ifaceq(m,jfaces))
655  elseif((lakon(nelems)(4:4).eq.'4').or.
656  & (lakon(nelems)(4:5).eq.'10')) then
657  noded=kon(indexe+ifacet(m,jfaces))
658  elseif(lakon(nelems)(4:4).eq.'6') then
659  noded=kon(indexe+ifacew1(m,jfaces))
660  elseif(lakon(nelems)(4:5).eq.'15') then
661  noded=kon(indexe+ifacew2(m,jfaces))
662  endif
663  else
664  if(i.gt.istartset(jdep)) then
665  if(ialset(i).eq.ialset(i-1)) cycle loop1
666  endif
667  noded=ialset(i)
668  endif
669 c enddo
670 !
671  xap=co(1,noded)-csab(1)
672  yap=co(2,noded)-csab(2)
673  zap=co(3,noded)-csab(3)
674 !
675  zp=xap*xn+yap*yn+zap*zn
676  rp=dsqrt((xap-zp*xn)**2+(yap-zp*yn)**2+(zap-zp*zn)**2)
677 !
678  if((.not.calcangle).and.(rp.gt.1.d-10)) then
679  x2=(xap-zp*xn)/rp
680  y2=(yap-zp*yn)/rp
681  z2=(zap-zp*zn)/rp
682  endif
683 !
684  call near2d(rcs0,zcs0,rcs,zcs,nr,nz,rp,zp,ncsnodes,node,
685  & nneigh)
686 !
687  nodei=ics(node)
688  if(nodei.lt.0) cycle
689  if(nodei.eq.noded) then
690  ics(node)=-nodei
691  nodesonaxis=.true.
692  cycle
693  endif
694 !
695  xap=co(1,nodei)-csab(1)
696  yap=co(2,nodei)-csab(2)
697  zap=co(3,nodei)-csab(3)
698 !
699  zp=xap*xn+yap*yn+zap*zn
700  rp=dsqrt((xap-zp*xn)**2+(yap-zp*yn)**2+(zap-zp*zn)**2)
701 !
702  if((.not.calcangle).and.(rp.gt.1.d-10)) then
703  x3=(xap-zp*xn)/rp
704  y3=(yap-zp*yn)/rp
705  z3=(zap-zp*zn)/rp
706 !
707  x1=y2*z3-y3*z2
708  y1=x3*z2-x2*z3
709  z1=x2*y3-x3*y2
710 !
711  phi=(x1*xn+y1*yn+z1*zn)/dabs(x1*xn+y1*yn+z1*zn)*
712  & dacos(x2*x3+y2*y3+z2*z3)
713  if(check) then
714  calculated_angle=dacos(x2*x3+y2*y3+z2*z3)
715  user_angle=6.28318531d0/cs(1,mcs)
716  if(dabs(calculated_angle-user_angle)/
717  & calculated_angle.gt.0.01d0) then
718  write(*,*)
719  & '*ERROR reading *CYCLIC SYMMETRY MODEL'
720  write(*,*) ' number of segments does not'
721  write(*,*) ' agree with the geometry'
722  write(*,*) ' angle based on N:',
723  & user_angle*57.29577951d0
724  write(*,*)' angle based on the geometry:',
725  & calculated_angle*57.29577951d0
726  call exit(201)
727  endif
728  else
729  write(*,*) '*INFO in cyclicsymmetrymodels: angle'
730  write(*,*)' check is deactivated by the user;'
731  write(*,*) ' the real geometry is used for'
732  write(*,*) ' the calculation of the segment'
733  write(*,*) ' angle'
734  write(*,*)
735  endif
736  calcangle=.true.
737  endif
738  enddo
739 !
740  else
741  k=ialset(i-2)
742  do
743  k=k-ialset(i)
744  if(k.ge.ialset(i-1)) exit
745  noded=k
746 !
747  xap=co(1,noded)-csab(1)
748  yap=co(2,noded)-csab(2)
749  zap=co(3,noded)-csab(3)
750 !
751  zp=xap*xn+yap*yn+zap*zn
752  rp=dsqrt((xap-zp*xn)**2+(yap-zp*yn)**2+(zap-zp*zn)**2)
753 !
754  if((.not.calcangle).and.(rp.gt.1.d-10)) then
755  x2=(xap-zp*xn)/rp
756  y2=(yap-zp*yn)/rp
757  z2=(zap-zp*zn)/rp
758  endif
759 !
760  call near2d(rcs0,zcs0,rcs,zcs,nr,nz,rp,zp,ncsnodes,node,
761  & nneigh)
762 !
763  nodei=ics(node)
764  if(nodei.lt.0) cycle
765  if(nodei.eq.noded) then
766  ics(node)=-nodei
767  nodesonaxis=.true.
768  cycle
769  endif
770 !
771  xap=co(1,nodei)-csab(1)
772  yap=co(2,nodei)-csab(2)
773  zap=co(3,nodei)-csab(3)
774 !
775  zp=xap*xn+yap*yn+zap*zn
776  rp=dsqrt((xap-zp*xn)**2+(yap-zp*yn)**2+(zap-zp*zn)**2)
777 !
778  if((.not.calcangle).and.(rp.gt.1.d-10)) then
779  x3=(xap-zp*xn)/rp
780  y3=(yap-zp*yn)/rp
781  z3=(zap-zp*zn)/rp
782 !
783  x1=y2*z3-y3*z2
784  y1=x3*z2-x2*z3
785  z1=x2*y3-x3*y2
786 !
787  phi=(x1*xn+y1*yn+z1*zn)/dabs(x1*xn+y1*yn+z1*zn)*
788  & dacos(x2*x3+y2*y3+z2*z3)
789  if(check) then
790  calculated_angle=dacos(x2*x3+y2*y3+z2*z3)
791  user_angle=6.28318531d0/cs(1,mcs)
792  if(dabs(calculated_angle-user_angle)
793  & /calculated_angle.gt.0.01d0) then
794  write(*,*)
795  & '*ERROR reading *CYCLIC SYMMETRY MODEL'
796  write(*,*) ' number of segments does not'
797  write(*,*) ' agree with the geometry'
798  write(*,*) ' angle based on N:',
799  & user_angle*57.29577951d0
800  write(*,*) ' angle based on the geometry:'
801  & ,calculated_angle*57.29577951d0
802  call exit(201)
803  endif
804  endif
805  calcangle=.true.
806  endif
807 !
808  enddo
809  endif
810 !
811  enddo loop1
812 !
813 ! allocating a node of the depset to each node of the indepset
814 !
815  ncounter=0
816  triangulation=.false.
817 !
818 ! opening a file to store the nodes which are not connected
819 !
820  open(40,file='WarnNodeMissCyclicSymmetry.nam',status='unknown')
821  write(40,*) '*NSET,NSET=WarnNodeCyclicSymmetry'
822  write(*,*) '*INFO in cyclicsymmetrymodels:'
823  write(*,*) ' failed nodes (if any) are stored in file'
824  write(*,*) ' WarnNodeMissCyclicSymmetry.nam'
825  write(*,*) ' This file can be loaded into'
826  write(*,*) ' an active cgx-session by typing'
827  write(*,*)
828  & ' read WarnNodeMissCyclicSymmetry.nam inp'
829  write(*,*)
830 !
831  loop2: do i=istartset(jdep),iendset(jdep)
832  if(ialset(i).gt.0) then
833 !
834 ! check whether dependent side is node based or
835 ! face based
836 !
837  if(depkind.eq.'T') then
838  ifaces=ialset(i)
839  nelems=int(ifaces/10)
840  jfaces=ifaces - nelems*10
841  indexe=ipkon(nelems)
842 !
843  if(lakon(nelems)(4:5).eq.'20') then
844  nopes=8
845  elseif(lakon(nelems)(4:4).eq.'2') then
846  nopes=9
847  elseif(lakon(nelems)(4:4).eq.'8') then
848  nopes=4
849  elseif(lakon(nelems)(4:5).eq.'10') then
850  nopes=6
851  elseif(lakon(nelems)(4:4).eq.'4') then
852  nopes=3
853  endif
854 !
855  if(lakon(nelems)(4:4).eq.'6') then
856  if(jfaces.le.2) then
857  nopes=3
858  else
859  nopes=4
860  endif
861  endif
862  if(lakon(nelems)(4:5).eq.'15') then
863  if(jfaces.le.2) then
864  nopes=6
865  else
866  nopes=8
867  endif
868  endif
869  else
870  nopes=1
871  endif
872 !
873  do m=1,nopes
874  if(depkind.eq.'T') then
875  if((lakon(nelems)(4:4).eq.'2').or.
876  & (lakon(nelems)(4:4).eq.'8')) then
877  noded=kon(indexe+ifaceq(m,jfaces))
878  elseif((lakon(nelems)(4:4).eq.'4').or.
879  & (lakon(nelems)(4:5).eq.'10')) then
880  noded=kon(indexe+ifacet(m,jfaces))
881  elseif(lakon(nelems)(4:4).eq.'6') then
882  noded=kon(indexe+ifacew1(m,jfaces))
883  elseif(lakon(nelems)(4:5).eq.'15') then
884  noded=kon(indexe+ifacew2(m,jfaces))
885  endif
886  else
887  if(i.gt.istartset(jdep)) then
888  if(ialset(i).eq.ialset(i-1)) cycle loop2
889  endif
890  noded=ialset(i)
891  endif
892 c enddo
893 !
894 ! check whether cyclic MPC's have already been
895 ! generated (e.g. for nodes belonging to several
896 ! faces for face based dependent surfaces)
897 !
898  idof=8*(noded-1)+1
899  call nident(ikmpc,idof,nmpc,id)
900  if(id.gt.0) then
901  if(ikmpc(id).eq.idof) then
902  if(labmpc(ilmpc(id))(1:6).eq.'CYCLIC') cycle
903  endif
904  endif
905 !
906  call generatecycmpcs(tolloc,co,nk,ipompc,nodempc,
907  & coefmpc,nmpc,ikmpc,ilmpc,mpcfree,rcs,zcs,ics,
908  & nr,nz,rcs0,zcs0,labmpc,
909  & mcs,triangulation,csab,xn,yn,zn,phi,noded,
910  & ncsnodes,rcscg,rcs0cg,zcscg,zcs0cg,nrcg,
911  & nzcg,jcs,lcs,kontri,straight,ne,ipkon,kon,lakon,
912  & ifacetet,inodface,ncounter,jobnamec,vold,cfd,mi,
913  & indepset)
914  enddo
915 !
916  else
917  k=ialset(i-2)
918  do
919  k=k-ialset(i)
920  if(k.ge.ialset(i-1)) exit
921  noded=k
922 !
923 ! check whether cyclic MPC's have already been
924 ! generated (e.g. for nodes belonging to several
925 ! faces for face based dependent surfaces)
926 !
927  idof=8*(noded-1)+1
928  call nident(ikmpc,idof,nmpc,id)
929  if(id.gt.0) then
930  if(ikmpc(id).eq.idof) then
931  if(labmpc(ilmpc(id))(1:6).eq.'CYCLIC') cycle
932  endif
933  endif
934 !
935  call generatecycmpcs(tolloc,co,nk,ipompc,nodempc,
936  & coefmpc,nmpc,ikmpc,ilmpc,mpcfree,rcs,zcs,ics,
937  & nr,nz,rcs0,zcs0,labmpc,
938  & mcs,triangulation,csab,xn,yn,zn,phi,noded,
939  & ncsnodes,rcscg,rcs0cg,zcscg,zcs0cg,nrcg,
940  & nzcg,jcs,lcs,kontri,straight,ne,ipkon,kon,lakon,
941  & ifacetet,inodface,ncounter,jobnamec,vold,cfd,mi,
942  & indepset)
943  enddo
944  endif
945 !
946  enddo loop2
947 !
948  close(40)
949 !
950  if(ncounter.ne.0) then
951  write(*,*) '*WARNING reading *CYCLIC SYMMETRY MODEL:'
952  write(*,*) ' for at least one dependent'
953  write(*,*) ' node in a cyclic symmetry definition no '
954  write(*,*) ' independent counterpart was found.'
955  write(*,*) ' Failed nodes are stored in file '
956  write(*,*) ' WarnNodeMissCyclicSymmetry.nam'
957 c next line was commented on 19/04/2012
958 c call exit(201)
959  endif
960 !
961 ! sorting ics
962 ! ics contains the master (independent) nodes
963 !
964  kflag=1
965  call isortii(ics,nr,ncsnodes,kflag)
966  cs(4,mcs)=ncsnodes+0.5
967  lprev=lprev+ncsnodes
968 !
969 ! check orientation of (xn,yn,zn) (important for copying of base
970 ! sector in arpackcs)
971 !
972  if(phi.lt.0.d0) then
973  csab(4)=2.d0*csab(1)-csab(4)
974  csab(5)=2.d0*csab(2)-csab(5)
975  csab(6)=2.d0*csab(3)-csab(6)
976  endif
977 !
978  do i=1,7
979  cs(5+i,mcs)=csab(i)
980  enddo
981 !
982  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
983  & ipoinp,inp,ipoinpc)
984 !
985  return
generatecycmpcs
subroutine generatecycmpcs(tolloc, co, nk, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, mpcfree, rcs, zcs, ics, nr, nz, rcs0, zcs0, labmpc, mcs, triangulation, csab, xn, yn, zn, phi, noded, ncsnodes, rcscg, rcs0cg, zcscg, zcs0cg, nrcg, nzcg, jcs, lcs, kontri, straight, ne, ipkon, kon, lakon, ifacetet, inodface, ncounter, jobnamec, vold, cfd, mi, indepset)
Definition: generatecycmpcs.f:26
z1
static double * z1
Definition: filtermain.c:48
max
#define max(a, b)
Definition: cascade.c:32
inputwarning
subroutine inputwarning(inpc, ipoinpc, iline, text)
Definition: inputwarning.f:20
x1
static double * x1
Definition: filtermain.c:48
isortii
subroutine isortii(ix, iy, n, kflag)
Definition: isortii.f:6
dist
static double * dist
Definition: radflowload.c:42
inputinfo
subroutine inputinfo(inpc, ipoinpc, iline, text)
Definition: inputinfo.f:20
getnewline
subroutine getnewline(inpc, textpart, istat, n, key, iline, ipol, inl, ipoinp, inp, ipoinpc)
Definition: getnewline.f:21
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
dsort
subroutine dsort(dx, iy, n, kflag)
Definition: dsort.f:6
inputerror
subroutine inputerror(inpc, ipoinpc, iline, text)
Definition: inputerror.f:20
near2d
subroutine near2d(xo, yo, x, y, nx, ny, xp, yp, n, neighbor, k)
Definition: near2d.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)