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 couptempdisps(inpc,textpart,nmethod,
20     &  iperturb,isolver,
21     &  istep,istat,n,tinc,tper,tmin,tmax,idrct,ithermal,iline,ipol,
22     &  inl,ipoinp,inp,ipoinpc,alpha,ctrl,iexpl,tincf,ttime,nener,
23     &  ier)
24!
25!     reading the input deck: *COUPLED TEMPERATURE-DISPLACEMENT
26!
27!     isolver=0: SPOOLES
28!             2: iterative solver with diagonal scaling
29!             3: iterative solver with Cholesky preconditioning
30!             4: sgi solver
31!             5: TAUCS
32!             7: pardiso
33!             8: pastix
34!
35!      iexpl==0:  structure:implicit, fluid:incompressible
36!      iexpl==1:  structure:implicit, fluid:compressible
37!
38      implicit none
39!
40      logical timereset
41!
42      character*1 inpc(*)
43      character*20 solver
44      character*132 textpart(16)
45!
46      integer nmethod,iperturb(*),isolver,istep,istat,n,key,i,idrct,
47     &  ithermal(*),iline,ipol,inl,ipoinp(2,*),inp(3,*),ipoinpc(0:*),
48     &  iexpl,nener,ier
49!
50      real*8 tinc,tper,tmin,tmax,alpha(*),ctrl(*),tincf,ttime
51!
52      idrct=0
53      alpha(1)=-0.05d0
54      tmin=0.d0
55      tmax=0.d0
56c      tincf=1.d-2
57      tincf=-1.d0
58      nmethod=4
59      timereset=.false.
60!
61      if(iperturb(1).eq.0) then
62         iperturb(1)=2
63      elseif((iperturb(1).eq.1).and.(istep.gt.1)) then
64         write(*,*) '*ERROR reading *COUPLED TEMPERATURE-DISPLACEMENT:'
65         write(*,*) '       perturbation analysis is not provided in a'
66         write(*,*) '       *COUPLED TEMPERATURE-DISPLACEMENT step.'
67         ier=1
68         return
69      endif
70!
71      if(istep.lt.1) then
72         write(*,*) '*ERROR reading *COUPLED TEMPERATURE-DISPLACEMENT:'
73         write(*,*) '       *COUPLED TEMPERATURE-DISPLACEMENT can only '
74         write(*,*) '       be used within a STEP'
75         ier=1
76         return
77      endif
78!
79!     default solver
80!
81      solver='                    '
82      if(isolver.eq.0) then
83         solver(1:7)='SPOOLES'
84      elseif(isolver.eq.2) then
85         solver(1:16)='ITERATIVESCALING'
86      elseif(isolver.eq.3) then
87         solver(1:17)='ITERATIVECHOLESKY'
88      elseif(isolver.eq.4) then
89         solver(1:3)='SGI'
90      elseif(isolver.eq.5) then
91         solver(1:5)='TAUCS'
92      elseif(isolver.eq.7) then
93         solver(1:7)='PARDISO'
94      elseif(isolver.eq.8) then
95         solver(1:6)='PASTIX'
96      endif
97!
98      do i=2,n
99         if(textpart(i)(1:6).eq.'ALPHA=') then
100            read(textpart(i)(7:26),'(f20.0)',iostat=istat) alpha(1)
101            if(istat.gt.0) then
102               call inputerror(inpc,ipoinpc,iline,
103     &              "*COUPLED TEMPERATURE-DISPLACEMENT%",ier)
104               return
105            endif
106            if(alpha(1).lt.-1.d0/3.d0) then
107               write(*,*) '*WARNING in dynamics: alpha is smaller'
108               write(*,*) '  than -1/3 and is reset to -1/3'
109               alpha(1)=-1.d0/3.d0
110            elseif(alpha(1).gt.0.d0) then
111               write(*,*) '*WARNING in dynamics: alpha is greater'
112               write(*,*) '  than 0 and is reset to 0'
113               alpha(1)=0.d0
114            endif
115         elseif(textpart(i)(1:7).eq.'SOLVER=') then
116            read(textpart(i)(8:27),'(a20)') solver
117         elseif(textpart(i)(1:12).eq.'COMPRESSIBLE') then
118            iexpl=1
119         elseif((textpart(i)(1:6).eq.'DIRECT').and.
120     &          (textpart(i)(1:9).ne.'DIRECT=NO')) then
121            idrct=1
122         elseif(textpart(i)(1:11).eq.'STEADYSTATE') then
123            nmethod=1
124         elseif(textpart(i)(1:7).eq.'DELTMX=') then
125            read(textpart(i)(8:27),'(f20.0)',iostat=istat) ctrl(27)
126         elseif(textpart(i)(1:9).eq.'TIMERESET') then
127            timereset=.true.
128         elseif(textpart(i)(1:17).eq.'TOTALTIMEATSTART=') then
129            read(textpart(i)(18:37),'(f20.0)',iostat=istat) ttime
130         else
131            write(*,*)
132     &        '*WARNING reading *COUPLED TEMPERATURE-DISPLACEMENT:'
133            write(*,*) '         parameter not recognized:'
134            write(*,*) '         ',
135     &                 textpart(i)(1:index(textpart(i),' ')-1)
136            call inputwarning(inpc,ipoinpc,iline,
137     &"*COUPLED TEMPERATURE-DISPLACEMENT%")
138         endif
139      enddo
140      if(nmethod.eq.1) ctrl(27)=1.d30
141!
142      if((ithermal(1).eq.0).and.(nmethod.ne.1).and.
143     &   (nmethod.ne.2).and.(iperturb(1).ne.0)) then
144         write(*,*) '*ERROR reading *COUPLED TEMPERATURE-DISPLACEMENT:'
145         write(*,*) '       please define initial '
146         write(*,*) '       conditions for the temperature'
147         ier=1
148         return
149      else
150         ithermal(1)=3
151      endif
152!
153      if(solver(1:7).eq.'SPOOLES') then
154         isolver=0
155      elseif(solver(1:16).eq.'ITERATIVESCALING') then
156         isolver=2
157      elseif(solver(1:17).eq.'ITERATIVECHOLESKY') then
158         isolver=3
159      elseif(solver(1:3).eq.'SGI') then
160         isolver=4
161      elseif(solver(1:5).eq.'TAUCS') then
162         isolver=5
163      elseif(solver(1:7).eq.'PARDISO') then
164         isolver=7
165      elseif(solver(1:6).eq.'PASTIX') then
166         isolver=8
167      else
168        write(*,*) '*WARNING reading *COUPLED TEMPERATURE-DISPLACEMENT:'
169         write(*,*) '         unknown solver;'
170         write(*,*) '         the default solver is used'
171      endif
172!
173      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
174     &     ipoinp,inp,ipoinpc)
175      if((istat.lt.0).or.(key.eq.1)) then
176         if(iperturb(1).ge.2) then
177            write(*,*)
178     &          '*WARNING reading *COUPLED TEMPERATURE-DISPLACEMENT:'
179            write(*,*)
180     &          '         a nonlinear geometricanalysis is requested'
181            write(*,*) '         but no time increment nor step is speci
182     &fied'
183            write(*,*) '         the defaults (1,1) are used'
184            tinc=1.d0
185            tper=1.d0
186            tmin=1.d-5
187            tmax=1.d+30
188            tincf=-1.d0
189c            tincf=1.d-2
190         endif
191         if(timereset)ttime=ttime-tper
192         return
193      endif
194!
195      read(textpart(1)(1:20),'(f20.0)',iostat=istat) tinc
196      if(istat.gt.0) then
197         call inputerror(inpc,ipoinpc,iline,
198     &        "*COUPLED TEMPERATURE-DISPLACEMENT%",ier)
199         return
200      endif
201      read(textpart(2)(1:20),'(f20.0)',iostat=istat) tper
202      if(istat.gt.0) then
203         call inputerror(inpc,ipoinpc,iline,
204     &        "*COUPLED TEMPERATURE-DISPLACEMENT%",ier)
205         return
206      endif
207      read(textpart(3)(1:20),'(f20.0)',iostat=istat) tmin
208      if(istat.gt.0) then
209         call inputerror(inpc,ipoinpc,iline,
210     &        "*COUPLED TEMPERATURE-DISPLACEMENT%",ier)
211         return
212      endif
213      read(textpart(4)(1:20),'(f20.0)',iostat=istat) tmax
214      if(istat.gt.0) then
215         call inputerror(inpc,ipoinpc,iline,
216     &        "*COUPLED TEMPERATURE-DISPLACEMENT%",ier)
217         return
218      endif
219      read(textpart(5)(1:20),'(f20.0)',iostat=istat) tincf
220      if(istat.gt.0) then
221         call inputerror(inpc,ipoinpc,iline,
222     &        "*COUPLED TEMPERATURE-DISPLACEMENT%",ier)
223         return
224      endif
225!
226      if(tinc.le.0.d0) then
227         write(*,*) '*ERROR reading *COUPLED TEMPERATURE-DISPLACEMENT:'
228         write(*,*) '       initial increment size is negative'
229      endif
230      if(tper.le.0.d0) then
231         write(*,*) '*ERROR reading *COUPLED TEMPERATURE-DISPLACEMENT:'
232         write(*,*) '       step size is negative'
233      endif
234      if(tinc.gt.tper) then
235         write(*,*) '*ERROR reading *COUPLED TEMPERATURE-DISPLACEMENT:'
236         write(*,*) '       initial increment size exceeds step size'
237      endif
238!
239      if(idrct.ne.1) then
240         if(dabs(tmin).lt.1.d-6*tper) then
241            tmin=min(tinc,1.d-6*tper)
242         endif
243         if(dabs(tmax).lt.1.d-10) then
244            tmax=1.d+30
245         endif
246         if(tinc.gt.dabs(tmax)) then
247            write(*,*)
248     &          '*WARNING reading *COUPLED TEMPERATURE-DISPLACEMENT:'
249            write(*,*) '         the initial increment ',tinc
250            write(*,*) '         exceeds the maximum increment ',
251     &          tmax
252            write(*,*) '         the initial increment is reduced'
253            write(*,*) '         to the maximum value'
254            tinc=dabs(tmax)
255         endif
256      endif
257!
258      if(timereset)ttime=ttime-tper
259!
260      if(nmethod.eq.4) nener=1
261!
262      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
263     &     ipoinp,inp,ipoinpc)
264!
265      return
266      end
267
268
269
270
271
272
273
274
275