1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 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 selectcyclicsymmetrymodess(inpc,textpart,cs,ics,
20     &  tieset,istartset,
21     &  iendset,ialset,ipompc,nodempc,coefmpc,nmpc,nmpc_,ikmpc,ilmpc,
22     &  mpcfree,mcs,set,nset,labmpc,istep,istat,n,iline,ipol,inl,
23     &  ipoinp,inp,nmethod,key,ipoinpc,ier)
24!
25!     reading the input deck: *SELECT CYCLIC SYMMETRY MODES
26!
27      implicit none
28!
29      character*1 inpc(*)
30      character*20 labmpc(*)
31      character*81 set(*),leftset,tieset(3,*)
32      character*132 textpart(16)
33!
34      integer istep,istat,n,key,i,ns(5),ics(*),istartset(*),ier,
35     &  iendset(*),ialset(*),id,ipompc(*),nodempc(3,*),nmpc,nmpc_,
36     &  ikmpc(*),ilmpc(*),mpcfree,i1(2),i2(2),i3,i4,i5,j,k,
37     &  mpcfreeold,idof,node,ileft,nset,irepeat,ipoinpc(0:*),
38     &  mpc,iline,ipol,inl,ipoinp(2,*),inp(3,*),mcs,lprev,ij,nmethod
39!
40      real*8 coefmpc(*),csab(7),x1(2),x2(2),x3,x4,x5,dd,xn,yn,zn,
41     &  cs(17,*)
42!
43!     irepeat indicates whether the step was preceded by another
44!     cyclic symmetry step (irepeat=1) or not (irepeat=0)
45!
46      data irepeat /0/
47      save irepeat
48!
49      if(istep.eq.0) then
50         write(*,*)'*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
51         write(*,*)'       *SELECT CYCLIC SYMMETRY MODES'
52         write(*,*)'       should be placed within a step definition'
53         ier=1
54         return
55      endif
56!
57!     check whether in case of cyclic symmetry the frequency procedure
58!     is chosen
59!
60      if((nmethod.ne.2).and.(nmethod.ne.13)) then
61         write(*,*) '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
62         write(*,*) '       the only valid procedures'
63         write(*,*) '       for cyclic symmetry calculations'
64         write(*,*) '       with nodal diameters are *FREQUENCY'
65         write(*,*) '       and *GREEN'
66         ier=1
67         return
68      endif
69!
70      ns(2)=0
71      ns(3)=0
72!
73      do i=2,n
74         if(textpart(i)(1:5).eq.'NMIN=') then
75            read(textpart(i)(6:15),'(i10)',iostat=istat) ns(2)
76            if(istat.gt.0) then
77               call inputerror(inpc,ipoinpc,iline,
78     &              "*SELECT CYCLIC SYMMETRY MODES%",ier)
79               return
80            endif
81         elseif(textpart(i)(1:5).eq.'NMAX=') then
82            read(textpart(i)(6:15),'(i10)',iostat=istat) ns(3)
83            if(istat.gt.0) then
84               call inputerror(inpc,ipoinpc,iline,
85     &              "*SELECT CYCLIC SYMMETRY MODES%",ier)
86               return
87            endif
88         else
89            write(*,*) '*WARNING reading *SELECT CYCLIC SYMMETRY MODES:'
90            write(*,*) '         parameter not recognized:'
91            write(*,*) '         ',
92     &                 textpart(i)(1:index(textpart(i),' ')-1)
93            call inputwarning(inpc,ipoinpc,iline,
94     &"*SELECT CYCLIC SYMMETRY MODES%")
95         endif
96      enddo
97!
98!     check the input
99!
100      if(ns(2).lt.0) then
101         ns(2)=0
102         write(*,*) '*WARNING reading *SELECT CYCLIC SYMMETRY MODES:'
103         write(*,*) '         minimum nodal'
104         write(*,*) '         diameter must be nonnegative'
105      endif
106      if(ns(3).lt.ns(2)) then
107         write(*,*) '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
108         write(*,*) '       maximum nodal'
109         write(*,*) '       diameter should not exceed minimal one'
110         ier=1
111         return
112      endif
113!
114!     loop over all cyclic symmetry parts
115!
116      do ij=1,mcs
117         ns(1)=int(cs(1,ij))
118         ns(4)=int(cs(4,ij))
119         leftset=tieset(2,int(cs(17,ij)))
120         lprev=int(cs(14,ij))
121         do i=1,7
122            csab(i)=cs(5+i,ij)
123         enddo
124!
125!     check whether cyclic symmetry axis is part of the structure
126!
127c         do i=1,nset
128c            if(set(i).eq.leftset) exit
129c         enddo
130c         ileft=i
131         call cident81(set,leftset,nset,id)
132         ileft=nset+1
133         if(id.gt.0) then
134           if(leftset.eq.set(id)) then
135             ileft=id
136           endif
137         endif
138!
139!     if this step was preceded by a cyclic symmetry step:
140!     check for MPC's for nodes on the cyclic symmetry axis
141!     and delete them
142!
143         if(irepeat.eq.1) then
144            do i=1,ns(4)
145               node=ics(lprev+i)
146               if(node.lt.0) then
147                  node=-node
148                  do k=1,3
149                     idof=8*(node-1)+k
150                     call nident(ikmpc,idof,nmpc,id)
151                     if(id.gt.0) then
152                        if(ikmpc(id).eq.idof) then
153c                           write(*,*) 'removing MPC',node,k
154                           mpc=ilmpc(id)
155                           call mpcrem(mpc,mpcfree,nodempc,nmpc,ikmpc,
156     &                           ilmpc,labmpc,coefmpc,ipompc)
157                        endif
158                     endif
159                  enddo
160               endif
161            enddo
162         endif
163!
164         do i=1,ns(4)
165            node=ics(lprev+i)
166            if(node.lt.0) then
167               node=-node
168               if(ns(2).ne.ns(3)) then
169                  if((ns(2).eq.0).or.(ns(2).eq.1)) then
170                     write(*,*) '*ERROR: axis of cyclic symmetry'
171                     write(*,*) '        is part of the structure;'
172                     write(*,*) '        nodal diameters 0, 1, and'
173                     write(*,*) '        those above must be each in'
174                     write(*,*) '        separate steps.'
175                     ier=1
176                     return
177                  endif
178               endif
179!
180!     specifying special MPC's for nodes on the axis
181!
182!     normal along the axis
183!
184               xn=csab(4)-csab(1)
185               yn=csab(5)-csab(2)
186               zn=csab(6)-csab(3)
187               dd=dsqrt(xn*xn+yn*yn+zn*zn)
188               xn=xn/dd
189               yn=yn/dd
190               zn=zn/dd
191!
192!     nodal diameter 0
193!
194               if(ns(2).eq.0) then
195                  if(dabs(xn).gt.1.d-10) then
196                     i1(1)=2
197                     i1(2)=3
198                     i2(1)=1
199                     i2(2)=1
200                     x1(1)=xn
201                     x1(2)=xn
202                     x2(1)=-yn
203                     x2(2)=-zn
204                  elseif(dabs(yn).gt.1.d-10) then
205                     i1(1)=1
206                     i1(2)=3
207                     i2(1)=2
208                     i2(2)=2
209                     x1(1)=yn
210                     x1(2)=yn
211                     x2(1)=-xn
212                     x2(2)=-zn
213                  elseif(dabs(zn).gt.1.d-10) then
214                     i1(1)=1
215                     i1(2)=2
216                     i2(1)=3
217                     i2(2)=3
218                     x1(1)=zn
219                     x1(2)=zn
220                     x2(1)=-xn
221                     x2(2)=-yn
222                  endif
223!
224!     generating two MPC's expressing that the nodes cannot
225!     move in planes perpendicular to the cyclic symmetry
226!     axis
227!
228                  do k=1,2
229                     idof=8*(node-1)+i1(k)
230                     call nident(ikmpc,idof,nmpc,id)
231                     if(id.gt.0) then
232                        if(ikmpc(id).eq.idof) then
233                           write(*,*)
234     &                 '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
235                           write(*,*) '       node',node,
236     &                          ' on cyclic symmetry'
237                           write(*,*) '       axis is used in other MPC'
238                           ier=1
239                           return
240                        endif
241                     endif
242                     nmpc=nmpc+1
243                     ipompc(nmpc)=mpcfree
244                     labmpc(nmpc)='                    '
245!
246!     updating ikmpc and ilmpc
247!
248                     do j=nmpc,id+2,-1
249                        ikmpc(j)=ikmpc(j-1)
250                        ilmpc(j)=ilmpc(j-1)
251                     enddo
252                     ikmpc(id+1)=idof
253                     ilmpc(id+1)=nmpc
254!
255                     nodempc(1,mpcfree)=node
256                     nodempc(2,mpcfree)=i1(k)
257                     coefmpc(mpcfree)=x1(k)
258                     mpcfree=nodempc(3,mpcfree)
259                     if(mpcfree.eq.0) then
260                        write(*,*)
261     &                  '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
262                        write(*,*) '       increase memmpc_'
263                        ier=1
264                        return
265                     endif
266                     nodempc(1,mpcfree)=node
267                     nodempc(2,mpcfree)=i2(k)
268                     coefmpc(mpcfree)=x2(k)
269                     mpcfreeold=mpcfree
270                     mpcfree=nodempc(3,mpcfree)
271                     if(mpcfree.eq.0) then
272                        write(*,*)
273     &                  '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
274                        write(*,*) '       increase memmpc_'
275                        ier=1
276                        return
277                     endif
278                     nodempc(3,mpcfreeold)=0
279                  enddo
280               elseif(ns(2).eq.1) then
281!
282!     nodal diameter 1
283!
284                  if(dabs(xn).gt.1.d-10) then
285                     i3=1
286                     i4=2
287                     i5=3
288                     x3=xn
289                     x4=yn
290                     x5=zn
291                  elseif(dabs(yn).gt.1.d-10) then
292                     i3=2
293                     i4=2
294                     i5=3
295                     x3=yn
296                     x4=xn
297                     x5=zn
298                  else
299                     i3=3
300                     i4=1
301                     i5=2
302                     x3=zn
303                     x4=xn
304                     x5=yn
305                  endif
306!
307!     generating one MPC expressing that the nodes should
308!     not move along the axis
309!
310                  idof=8*(node-1)+i3
311                  call nident(ikmpc,idof,nmpc,id)
312                  if(id.gt.0) then
313                     if(ikmpc(id).eq.idof) then
314                        write(*,*)
315     &                  '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
316                        write(*,*) '       node',node,
317     &                       ' on cyclic symmetry'
318                        write(*,*) '       axis is used in other MPC'
319                        ier=1
320                        return
321                     endif
322                  endif
323                  nmpc=nmpc+1
324                  ipompc(nmpc)=mpcfree
325                  labmpc(nmpc)='                    '
326!
327!     updating ikmpc and ilmpc
328!
329                  do j=nmpc,id+2,-1
330                     ikmpc(j)=ikmpc(j-1)
331                     ilmpc(j)=ilmpc(j-1)
332                  enddo
333                  ikmpc(id+1)=idof
334                  ilmpc(id+1)=nmpc
335!
336                  nodempc(1,mpcfree)=node
337                  nodempc(2,mpcfree)=i3
338                  coefmpc(mpcfree)=x3
339                  mpcfree=nodempc(3,mpcfree)
340                  if(mpcfree.eq.0) then
341                     write(*,*)
342     &               '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
343                     write(*,*) '       increase memmpc_'
344                     ier=1
345                     return
346                  endif
347                  nodempc(1,mpcfree)=node
348                  nodempc(2,mpcfree)=i4
349                  coefmpc(mpcfree)=x4
350                  mpcfree=nodempc(3,mpcfree)
351                  if(mpcfree.eq.0) then
352                     write(*,*)
353     &               '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
354                     write(*,*) '       increase memmpc_'
355                     ier=1
356                     return
357                  endif
358                  nodempc(1,mpcfree)=node
359                  nodempc(2,mpcfree)=i5
360                  coefmpc(mpcfree)=x5
361                  mpcfreeold=mpcfree
362                  mpcfree=nodempc(3,mpcfree)
363                  if(mpcfree.eq.0) then
364                     write(*,*)
365     &               '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
366                     write(*,*) '       increase memmpc_'
367                     ier=1
368                     return
369                  endif
370                  nodempc(3,mpcfreeold)=0
371               else
372                  do k=1,3
373                     idof=8*(node-1)+k
374                     call nident(ikmpc,idof,nmpc,id)
375                     if(id.gt.0) then
376                        if(ikmpc(id).eq.idof) then
377                           write(*,*)
378     &                   '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
379                           write(*,*) '       node',node,
380     &                          ' on cyclic symmetry'
381                           write(*,*) '       axis is used in other MPC'
382                           ier=1
383                           return
384                        endif
385                     endif
386                     nmpc=nmpc+1
387                     ipompc(nmpc)=mpcfree
388                     labmpc(nmpc)='                    '
389!
390!     updating ikmpc and ilmpc
391!
392                     do j=nmpc,id+2,-1
393                        ikmpc(j)=ikmpc(j-1)
394                        ilmpc(j)=ilmpc(j-1)
395                     enddo
396                     ikmpc(id+1)=idof
397                     ilmpc(id+1)=nmpc
398!
399                     nodempc(1,mpcfree)=node
400                     nodempc(2,mpcfree)=k
401                     coefmpc(mpcfree)=1.d0
402                     mpcfreeold=mpcfree
403                     mpcfree=nodempc(3,mpcfree)
404                     if(mpcfree.eq.0) then
405                        write(*,*)
406     &                  '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:'
407                        write(*,*) '       increase memmpc_'
408                        ier=1
409                        return
410                     endif
411                     nodempc(3,mpcfreeold)=0
412                  enddo
413               endif
414            endif
415         enddo
416!
417         cs(2,ij)=ns(2)+0.5d0
418         cs(3,ij)=ns(3)+0.5d0
419      enddo
420!
421      if(irepeat.eq.0) irepeat=1
422!
423      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
424     &     ipoinp,inp,ipoinpc)
425!
426c      do j=1,nmpc
427c         call writempc(ipompc,nodempc,coefmpc,labmpc,j)
428c      enddo
429!
430      return
431      end
432
433