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 modaldynamics(inpc,textpart,nmethod,tinc,tper,iexpl,
20     &  istep,istat,n,iline,ipol,inl,ipoinp,inp,iperturb,isolver,
21     &  cs,mcs,ipoinpc,idrct,ctrl,tmin,tmax,nforc,nload,nbody,iprestr,
22     &  t0,t1,ithermal,nk,vold,veold,xmodal,set,nset,mi,cyclicsymmetry,
23     &  ier)
24!
25!     reading the input deck: *MODAL DYNAMIC
26!
27      implicit none
28!
29      logical steadystate,nodalset
30!
31      character*1 inpc(*)
32      character*20 solver
33      character*81 set(*),noset
34      character*132 textpart(16)
35!
36      integer nmethod,istep,istat,n,key,iexpl,iline,ipol,inl,
37     &  ipoinp(2,*),inp(3,*),iperturb(*),isolver,i,mcs,ipoinpc(0:*),
38     &  idrct,nforc,nload,nbody,iprestr,ithermal(*),j,nk,ipos,nset,
39     &  mi(*),cyclicsymmetry,ier
40!
41      real*8 tinc,tper,cs(17,*),ctrl(*),tmin,tmax,t0(*),t1(*),
42     &  vold(0:mi(2),*),veold(0:mi(2),*),xmodal(*)
43!
44      iexpl=0
45c      iperturb(1)=0
46      iperturb(2)=0
47      idrct=1
48      tmin=0.d0
49      tmax=0.d0
50      steadystate=.false.
51      if((mcs.ne.0).and.(cs(2,1).ge.0.d0)) then
52         cyclicsymmetry=1
53      endif
54      nodalset=.false.
55!
56      if(istep.lt.1) then
57         write(*,*)
58     &      '*ERROR reading *MODAL DYNAMIC: *MODAL DYNAMIC can only'
59         write(*,*) '  be used within a STEP'
60         ier=1
61         return
62      endif
63!
64!     default solver
65!
66      solver='                    '
67      if(isolver.eq.0) then
68         solver(1:7)='SPOOLES'
69      elseif(isolver.eq.2) then
70         solver(1:16)='ITERATIVESCALING'
71      elseif(isolver.eq.3) then
72         solver(1:17)='ITERATIVECHOLESKY'
73      elseif(isolver.eq.4) then
74         solver(1:3)='SGI'
75      elseif(isolver.eq.5) then
76         solver(1:5)='TAUCS'
77      elseif(isolver.eq.7) then
78         solver(1:7)='PARDISO'
79      elseif(isolver.eq.8) then
80         solver(1:6)='PASTIX'
81      endif
82!
83      do i=2,n
84         if(textpart(i)(1:7).eq.'SOLVER=') then
85            read(textpart(i)(8:27),'(a20)') solver
86         elseif(textpart(i)(1:9).eq.'DIRECT=NO') then
87            idrct=0
88         elseif(textpart(i)(1:7).eq.'DELTMX=') then
89            read(textpart(i)(8:27),'(f20.0)',iostat=istat) ctrl(27)
90         elseif(textpart(i)(1:11).eq.'STEADYSTATE') then
91            steadystate=.true.
92         else
93            write(*,*)
94     &      '*WARNING reading *MODAL DYNAMIC: parameter not recognized:'
95            write(*,*) '         ',
96     &                 textpart(i)(1:index(textpart(i),' ')-1)
97            call inputwarning(inpc,ipoinpc,iline,
98     &"*MODAL DYNAMIC%")
99         endif
100      enddo
101!
102      if(solver(1:7).eq.'SPOOLES') then
103         isolver=0
104      elseif(solver(1:16).eq.'ITERATIVESCALING') then
105         write(*,*)
106     &      '*WARNING reading *MODAL DYNAMIC: the iterative scaling'
107         write(*,*) '         procedure is not available for modal'
108         write(*,*) '         dynamic calculations; the default solver'
109         write(*,*) '         is used'
110      elseif(solver(1:17).eq.'ITERATIVECHOLESKY') then
111         write(*,*)
112     &       '*WARNING reading *MODAL DYNAMIC: the iterative scaling'
113         write(*,*) '         procedure is not available for modal'
114         write(*,*) '         dynamic calculations; the default solver'
115         write(*,*) '         is used'
116      elseif(solver(1:3).eq.'SGI') then
117         isolver=4
118      elseif(solver(1:5).eq.'TAUCS') then
119         isolver=5
120      elseif(solver(1:7).eq.'PARDISO') then
121         isolver=7
122      elseif(solver(1:6).eq.'PARDISO') then
123         isolver=8
124      else
125         write(*,*) '*WARNING reading *MODAL DYNAMIC: unknown solver;'
126         write(*,*) '         the default solver is used'
127      endif
128!
129c      if((isolver.eq.2).or.(isolver.eq.3)) then
130c        write(*,*) '*ERROR reading *MODAL DYNAMIC: the default solver ',
131c     & solver
132c         write(*,*) '       cannot be used for modal dynamic'
133c         write(*,*) '       calculations '
134c         ier=1
135c         return
136c      endif
137!
138      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
139     &     ipoinp,inp,ipoinpc)
140      if((istat.lt.0).or.(key.eq.1)) then
141         write(*,*)
142     &       '*ERROR reading *MODAL DYNAMIC: definition not complete'
143         write(*,*) '       '
144         call inputerror(inpc,ipoinpc,iline,
145     &        "*MODAL DYNAMIC%",ier)
146         return
147      endif
148      read(textpart(1)(1:20),'(f20.0)',iostat=istat)tinc
149      if(istat.gt.0) then
150         call inputerror(inpc,ipoinpc,iline,
151     &        "*MODAL DYNAMIC%",ier)
152         return
153      endif
154      read(textpart(2)(1:20),'(f20.0)',iostat=istat)tper
155      if(istat.gt.0) then
156         call inputerror(inpc,ipoinpc,iline,
157     &        "*MODAL DYNAMIC%",ier)
158         return
159      endif
160      read(textpart(3)(1:20),'(f20.0)',iostat=istat) tmin
161      if(istat.gt.0) then
162         call inputerror(inpc,ipoinpc,iline,
163     &        "*MODAL DYNAMIC%",ier)
164         return
165      endif
166      read(textpart(4)(1:20),'(f20.0)',iostat=istat) tmax
167      if(istat.gt.0) then
168         call inputerror(inpc,ipoinpc,iline,
169     &        "*MODAL DYNAMIC%",ier)
170         return
171      endif
172!
173      if(steadystate) then
174!
175!        modal dynamics calculation till steady state
176!
177         if(tper.le.0.d0) then
178            write(*,*) '*ERROR reading *MODAL DYNAMIC: relative error'
179            write(*,*) '       is nonpositive'
180            ier=1
181            return
182         endif
183         tper=-tper
184         if(tinc.le.0.d0) then
185            write(*,*)
186     &         '*ERROR reading *MODAL DYNAMIC: initial increment'
187            write(*,*) '       size is nonpositive'
188            ier=1
189            return
190         endif
191         if(tmin.lt.0.d0) then
192            tmin=1.d-10
193         endif
194         if(tmax.lt.1.d-10) then
195            tmax=1.d+30
196         endif
197      else
198!
199!        transient modal dynamics calculation
200!
201         if(tper.lt.0.d0) then
202            write(*,*)
203     &        '*ERROR reading *MODAL DYNAMIC: step size is negative'
204            ier=1
205            return
206         elseif(tper.le.0.d0) then
207            tper=1.d0
208         endif
209         if(tinc.lt.0.d0) then
210            write(*,*)
211     &            '*ERROR reading *MODAL DYNAMIC: initial increment size
212     &is negative'
213            ier=1
214            return
215         elseif(tinc.le.0.d0) then
216            tinc=tper
217         endif
218         if(tinc.gt.tper) then
219            write(*,*)
220     &            '*ERROR reading *MODAL DYNAMIC: initial increment size
221     &exceeds step size'
222            ier=1
223            return
224         endif
225!
226         if(idrct.ne.1) then
227            if(tmin.lt.1.d-10*tper) then
228               tmin=min(tinc,1.d-10*tper)
229            endif
230            if(tmax.lt.1.d-10) then
231               tmax=1.d+30
232            endif
233         endif
234      endif
235!
236!     removing the present loading
237!
238c      nforc=0
239c      nload=0
240c      nbody=0
241c      iprestr=0
242c      if((ithermal.eq.1).or.(ithermal.eq.3)) then
243c         do j=1,nk
244c            t1(j)=t0(j)
245c         enddo
246c      endif
247!
248!     resetting fields vold and veold after a frequency or
249!     buckling step
250!
251      if((nmethod.eq.2).or.(nmethod.eq.3)) then
252         do i=1,nk
253            do j=1,3
254               vold(j,i)=0.d0
255               veold(j,i)=0.d0
256            enddo
257         enddo
258      endif
259!
260      nmethod=4
261!
262!     correction for cyclic symmetric structures:
263!     if the present step was not preceded by a frequency step
264!     no nodal diameter has been selected. To make sure that
265!     mastructcs is called instead of mastruct a fictitious
266!     minimum nodal diameter is stored
267!
268      if((cyclicsymmetry.eq.1).and.(mcs.ne.0).and.(cs(2,1).lt.0.d0))
269     &       cs(2,1)=0.d0
270!
271      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
272     &     ipoinp,inp,ipoinpc)
273!
274      return
275      end
276