← Back to branch summary

~gabriel1984sibiu/calculix/ccx

« back to all changes in this revision

Viewing changes to CalculiX/ccx_2.11/src/cyclicsymmetrymodels.f

  • Committer: Grevutiu Gabriel
  • Date: 2016-12-30 12:06:41 UTC
  • Revision ID: gabriel1984sibiu@gmail.com-20161230120641-kzmhfy8mn00w3mhg
New upstream version

Show diffs side-by-side

added added

removed removed

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