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 steadystatedynamicss(inpc,textpart,nmethod,
20     &  iexpl,istep,istat,n,iline,ipol,inl,ipoinp,inp,iperturb,isolver,
21     &  xmodal,cs,mcs,ipoinpc,nforc,nload,nbody,iprestr,t0,t1,ithermal,
22     &  nk,set,nset,cyclicsymmetry,ibody,ier)
23!
24!     reading the input deck: *STEADY STATE DYNAMICS
25!
26      implicit none
27!
28      logical nodalset
29!
30      character*1 inpc(*)
31      character*3 harmonic
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,nset,
37     &  ipoinp(2,*),inp(3,*),iperturb(*),isolver,i,ndata,nfour,mcs,
38     &  ipoinpc(0:*),nforc,nload,nbody,iprestr,ithermal(*),j,nk,ipos,
39     &  cyclicsymmetry,ier,ibody(3,*)
40!
41      real*8 fmin,fmax,bias,tmin,tmax,xmodal(*),cs(17,*),t0(*),t1(*)
42!
43      iexpl=0
44      iperturb(2)=0
45      harmonic='YES'
46      if((mcs.ne.0).and.(cs(2,1).ge.0.d0)) then
47         cyclicsymmetry=1
48      endif
49      nodalset=.false.
50!
51      if(istep.lt.1) then
52         write(*,*) '*ERROR reading *STEADY STATE DYNAMICS:'
53         write(*,*) '               *STEADY STATE DYNAMICS'
54         write(*,*) '  can only be used within a STEP'
55         ier=1
56         return
57      endif
58!
59!     default solver
60!
61      solver='                    '
62      if(isolver.eq.0) then
63         solver(1:7)='SPOOLES'
64      elseif(isolver.eq.2) then
65         solver(1:16)='ITERATIVESCALING'
66      elseif(isolver.eq.3) then
67         solver(1:17)='ITERATIVECHOLESKY'
68      elseif(isolver.eq.4) then
69         solver(1:3)='SGI'
70      elseif(isolver.eq.5) then
71         solver(1:5)='TAUCS'
72      elseif(isolver.eq.7) then
73         solver(1:7)='PARDISO'
74      elseif(isolver.eq.8) then
75         solver(1:6)='PASTIX'
76      endif
77!
78      do i=2,n
79         if(textpart(i)(1:7).eq.'SOLVER=') then
80            read(textpart(i)(8:27),'(a20)') solver
81         elseif(textpart(i)(1:9).eq.'HARMONIC=') then
82            read(textpart(i)(10:12),'(a3)') harmonic
83         else
84            write(*,*) '*WARNING reading *STEADY STATE DYNAMICS:'
85            write(*,*) '         parameter not recognized:'
86            write(*,*) '         ',
87     &                 textpart(i)(1:index(textpart(i),' ')-1)
88            call inputwarning(inpc,ipoinpc,iline,
89     &"*STEADY STATE DYNAMICS%")
90         endif
91      enddo
92!
93      if(solver(1:7).eq.'SPOOLES') then
94         isolver=0
95      elseif(solver(1:16).eq.'ITERATIVESCALING') then
96         write(*,*) '*WARNING reading *STEADY STATE DYNAMICS:'
97         write(*,*) '         the iterative scaling'
98         write(*,*) '         procedure is not available for modal'
99         write(*,*) '         dynamic calculations; the default solver'
100         write(*,*) '         is used'
101      elseif(solver(1:17).eq.'ITERATIVECHOLESKY') then
102         write(*,*) '*WARNING reading *STEADY STATE DYNAMICS:'
103         write(*,*) '         the iterative scaling'
104         write(*,*) '         procedure is not available for modal'
105         write(*,*) '         dynamic calculations; the default solver'
106         write(*,*) '         is used'
107      elseif(solver(1:3).eq.'SGI') then
108         isolver=4
109      elseif(solver(1:5).eq.'TAUCS') then
110         isolver=5
111      elseif(solver(1:7).eq.'PARDISO') then
112         isolver=7
113      elseif(solver(1:6).eq.'PASTIX') then
114         isolver=8
115      else
116         write(*,*) '*WARNING reading *STEADY STATE DYNAMICS:'
117         write(*,*) '         unknown solver;'
118         write(*,*) '         the default solver is used'
119      endif
120!
121c      if((isolver.eq.2).or.(isolver.eq.3)) then
122c         write(*,*) '*ERROR reading *STEADY STATE DYNAMICS:'
123c         write(*,*) '       the default solver ',solver
124c         write(*,*) '       cannot be used for modal dynamic'
125c         write(*,*) '       calculations '
126c         ier=1
127c         return
128c      endif
129!
130      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
131     &     ipoinp,inp,ipoinpc)
132      if((istat.lt.0).or.(key.eq.1)) then
133         write(*,*) '*ERROR reading *STEADY STATE DYNAMICS:'
134         write(*,*) '       definition not complete'
135         write(*,*) '       '
136         call inputerror(inpc,ipoinpc,iline,
137     &        "*STEADY STATE DYNAMICS%",ier)
138         return
139      endif
140      read(textpart(1)(1:20),'(f20.0)',iostat=istat) fmin
141      if(istat.gt.0) then
142         call inputerror(inpc,ipoinpc,iline,
143     &        "*STEADY STATE DYNAMICS%",ier)
144         return
145      endif
146      xmodal(3)=fmin
147      read(textpart(2)(1:20),'(f20.0)',iostat=istat) fmax
148      if(istat.gt.0) then
149         call inputerror(inpc,ipoinpc,iline,
150     &        "*STEADY STATE DYNAMICS%",ier)
151         return
152      endif
153      xmodal(4)=fmax
154      read(textpart(3)(1:20),'(i10)',iostat=istat) ndata
155      if(istat.gt.0) then
156         call inputerror(inpc,ipoinpc,iline,
157     &        "*STEADY STATE DYNAMICS%",ier)
158         return
159      endif
160      if(textpart(3)(1:1).eq.' ') then
161         ndata=20
162      elseif(ndata.lt.2) then
163         write(*,*) '*WARNING reading *STEADY STATE DYNAMICS:'
164         write(*,*) '         number of data points:',ndata
165         write(*,*) '         is smaller than 2;'
166         write(*,*) '         2 data points are taken instead'
167         ndata=2
168      endif
169      xmodal(5)=ndata+0.5d0
170      read(textpart(4)(1:20),'(f20.0)',iostat=istat) bias
171      if(istat.gt.0) then
172         call inputerror(inpc,ipoinpc,iline,
173     &        "*STEADY STATE DYNAMICS%",ier)
174         return
175      endif
176      if(bias.lt.1.) bias=3.d0
177      xmodal(6)=bias
178!
179      if(harmonic.eq.'YES') then
180         xmodal(7)=-0.5d0
181      else
182         read(textpart(5)(1:10),'(i10)',iostat=istat) nfour
183         if(istat.gt.0) then
184            call inputerror(inpc,ipoinpc,iline,
185     &           "*STEADY STATE DYNAMICS%",ier)
186            return
187         endif
188         if(nfour.le.0) nfour=20
189         if(n.ge.6) then
190            read(textpart(6)(1:20),'(f20.0)',iostat=istat) tmin
191            if(istat.gt.0) then
192               call inputerror(inpc,ipoinpc,iline,
193     &              "*STEADY STATE DYNAMICS%",ier)
194               return
195            endif
196         else
197            tmin=0.d0
198         endif
199         if(n.ge.7) then
200            read(textpart(7)(1:20),'(f20.0)',iostat=istat) tmax
201            if(istat.gt.0) then
202               call inputerror(inpc,ipoinpc,iline,
203     &              "*STEADY STATE DYNAMICS%",ier)
204               return
205            endif
206         else
207            tmax=1.d0
208         endif
209         xmodal(7)=nfour+0.5d0
210         xmodal(8)=tmin
211         xmodal(9)=tmax
212      endif
213!
214!     removing the present loading except centrigufal loading
215!
216      nforc=0
217      nload=0
218      do i=1,nbody
219         if(ibody(1,i).ne.1) ibody(1,i)=0
220      enddo
221c      nbody=0
222      iprestr=0
223      if((ithermal(1).eq.1).or.(ithermal(1).eq.3)) then
224         do j=1,nk
225            t1(j)=t0(j)
226         enddo
227      endif
228!
229      nmethod=5
230!
231!     correction for cyclic symmetric structures:
232!       if the present step was not preceded by a frequency step
233!       no nodal diameter has been selected. To make sure that
234!       mastructcs is called instead of mastruct a fictitious
235!       minimum nodal diameter is stored
236!
237      if((cyclicsymmetry.eq.1).and.(mcs.ne.0).and.(cs(2,1)<0.d0))
238     &         cs(2,1)=0.d0
239!
240      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
241     &     ipoinp,inp,ipoinpc)
242!
243      return
244      end
245
246