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 creeps(inpc,textpart,nelcon,nmat,ntmat_,npmat_,
20     &        plicon,nplicon,elcon,iplas,iperturb,nstate_,ncmat_,
21     &        matname,irstrt,istep,istat,n,iline,ipol,inl,ipoinp,inp,
22     &        ipoinpc,ianisoplas,ier)
23!
24!     reading the input deck: *CREEP
25!
26      implicit none
27!
28      logical iso
29!
30      character*1 inpc(*)
31      character*80 matname(*)
32      character*132 textpart(16)
33!
34      integer nelcon(2,*),nmat,ntmat_,ntmat,istep,npmat_,nstate_,
35     &  n,key,i,j,iplas,iperturb(*),istat,nplicon(0:ntmat_,*),ncmat_,
36     &  k,id,irstrt(*),iline,ipol,inl,ipoinp(2,*),inp(3,*),ipoinpc(0:*),
37     &  ianisoplas,ier
38!
39      real*8 temperature,elcon(0:ncmat_,ntmat_,*),t1l,
40     &  plicon(0:2*npmat_,ntmat_,*)
41!
42      iso=.true.
43      ntmat=0
44!
45      if((istep.gt.0).and.(irstrt(1).ge.0)) then
46         write(*,*) '*ERROR reading *CREEP: *CREEP should be placed'
47         write(*,*) '  before all step definitions'
48         ier=1
49         return
50      endif
51!
52      if(nmat.eq.0) then
53         write(*,*) '*ERROR reading *CREEP: *CREEP should be preceded'
54         write(*,*) '  by a *MATERIAL card'
55         ier=1
56         return
57      endif
58!
59!     check for anisotropic creep: assumes a ucreep routine
60!
61!     following if corresponds to "elastic isotropic with or
62!     without plasticity"
63!
64      if((nelcon(1,nmat).ne.2).and.(nelcon(1,nmat).ne.-51)) then
65!
66!        following if corresponds to "elastic anisotropic with or
67!        without plasticity"
68!
69         if((nelcon(1,nmat).ne.9).and.(nelcon(1,nmat).ne.-114)) then
70            write(*,*) '*ERROR reading *CREEP: *CREEP should be'
71            write(*,*) '       preceded by an *ELASTIC,TYPE=ISO card,'
72            write(*,*) '       or an *ELASTIC,TYPE=ORTHO card'
73            ier=1
74            return
75         endif
76!
77         ianisoplas=1
78!
79         if(nelcon(1,nmat).ne.-114) then
80!
81!           viscoplastic material with zero yield surface and
82!           without hardening: no plasticity
83!
84            iperturb(1)=3
85            nelcon(1,nmat)=-114
86            do i=2,n
87               if(textpart(i)(1:8).eq.'LAW=USER') then
88                  nelcon(1,nmat)=-109
89                  exit
90               endif
91            enddo
92            if(nelcon(1,nmat).eq.-109) then
93!
94!              elastic orthotropic
95!              no plasticity
96!              user creep: -109
97!
98!              7 state variables: cf. Section 6.8.13
99!              "Elastic anisotropy with isotropic creep defined by
100!               a creep user subroutine" in the User's Manual
101!
102               nstate_=max(nstate_,7)
103               if(matname(nmat)(70:80).ne.'           ') then
104                  write(*,*) '*ERROR reading *CREEP: the material name'
105                  write(*,*) '       for an elastically anisotropic'
106                  write(*,*) '       material with isotropic creep must'
107                  write(*,*) '       not exceed 69 characters'
108                  ier=1
109                  return
110               else
111                  do i=80,12,-1
112                     matname(nmat)(i:i)=matname(nmat)(i-11:i-11)
113                  enddo
114                  matname(nmat)(1:11)='ANISO_CREEP'
115               endif
116               call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
117     &              ipoinp,inp,ipoinpc)
118               return
119            else
120!
121!              elastic orthotropic
122!              no plasticity
123!              Norton creep: -114
124!
125!              14 state variables: cf. Section 6.8.12
126!              "Elastic anisotropy with isotropic viscoplasticity"
127!              in the User's Manual
128!
129               nstate_=max(nstate_,14)
130               do i=1,nelcon(2,nmat)
131                  elcon(10,i,nmat)=0.d0
132                  elcon(11,i,nmat)=0.d0
133                  elcon(12,i,nmat)=0.d0
134               enddo
135               if(matname(nmat)(71:80).ne.'          ') then
136                  write(*,*) '*ERROR reading *CREEP: the material name'
137                  write(*,*) '       for an elastically anisotropic'
138                  write(*,*) '       material with Norton creep'
139                  write(*,*) '       must not exceed 70 characters'
140                  ier=1
141                  return
142               else
143                  do i=80,11,-1
144                     matname(nmat)(i:i)=matname(nmat)(i-10:i-10)
145                  enddo
146                  matname(nmat)(1:10)='ANISO_PLAS'
147               endif
148            endif
149         else
150!
151!              elastic orthotropic
152!              plasticity
153!              Norton creep: -114 (user creep is not allowed)
154!
155!              14 state variables: cf. Section 6.8.12
156!              "Elastic anisotropy with isotropic viscoplasticity"
157!              in the User's Manual
158!
159            do i=2,n
160               if(textpart(i)(1:8).eq.'LAW=USER') then
161                  write(*,*) '*ERROR reading *CREEP: for an elastically'
162                  write(*,*) '       anisotropic material with von'
163                  write(*,*) '       Mises plasticity only Norton creep'
164                  write(*,*) '       is allowed (no user subroutine)'
165                  ier=1
166                  return
167               endif
168            enddo
169         endif
170      endif
171!
172!     if the *CREEP card is not preceded by a *PLASTIC card, a zero
173!     yield surface is assumed
174!
175      if(nelcon(1,nmat).ne.-114) then
176!
177!        elastic isotropic
178!        plasticity or no plasticity
179!        creep (Norton or user): -52
180!
181         if(nelcon(1,nmat).ne.-51) then
182!
183!           elastic isotropic
184!           no plasticity -> zero yield plasticity
185!           creep (Norton or user)
186!
187            nplicon(0,nmat)=1
188            nplicon(1,nmat)=2
189            plicon(0,1,nmat)=0.d0
190            plicon(1,1,nmat)=0.d0
191            plicon(2,1,nmat)=0.d0
192            plicon(3,1,nmat)=0.d0
193            plicon(4,1,nmat)=10.d10
194         endif
195!
196         iperturb(1)=3
197         iplas=1
198         nelcon(1,nmat)=-52
199c         nstate_=max(nstate_,13)
200!
201!        13 state variables: cf. Sections 6.8.6 and 6.8.7
202!        "Incremental (visco)plasticity: multiplicative decomposition"
203!        and "Incremental (visco)plasticity: additive decomposition"
204!        in the User's Manual
205!
206!        one extra internal variable to store the creep strain rate
207!        (needed for the treatment of cetol)
208!
209         nstate_=max(nstate_,14)
210!
211         do i=2,n
212            if(textpart(i)(1:8).eq.'LAW=USER') then
213               do j=1,nelcon(2,nmat)
214                  elcon(3,j,nmat)=-1.d0
215               enddo
216               call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
217     &              ipoinp,inp,ipoinpc)
218               return
219            elseif(textpart(i)(1:10).eq.'LAW=NORTON') then
220!
221!              default; nothing to do
222!
223            else
224               write(*,*)
225     &              '*WARNING reading *CREEP: parameter not recognized:'
226               write(*,*) '         ',
227     &              textpart(i)(1:index(textpart(i),' ')-1)
228               call inputwarning(inpc,ipoinpc,iline,
229     &"*CREEP%")
230            endif
231         enddo
232!
233!        before interpolation: data are stored in positions 6-9:
234!        A,n,m,temperature
235!        after interpolation: data are stored in positions 3-5:
236!        A,n,m
237!
238         do
239            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
240     &           ipoinp,inp,ipoinpc)
241            if((istat.lt.0).or.(key.eq.1)) exit
242            ntmat=ntmat+1
243            if(ntmat.gt.ntmat_) then
244               write(*,*) '*ERROR reading *CREEP: increase ntmat_'
245               ier=1
246               return
247            endif
248            do i=1,3
249               read(textpart(i)(1:20),'(f20.0)',iostat=istat)
250     &               elcon(i+5,ntmat,nmat)
251               if(istat.gt.0) then
252                  call inputerror(inpc,ipoinpc,iline,
253     &                 "*CREEP%",ier)
254                  return
255               endif
256            enddo
257            if(elcon(6,ntmat,nmat).le.0.d0) then
258               write(*,*) '*ERROR reading *CREEP: parameter A'
259               write(*,*) '       in the Norton law is nonpositive'
260               ier=1
261               return
262            endif
263            if(elcon(7,ntmat,nmat).le.0.d0) then
264               write(*,*) '*ERROR reading *CREEP: parameter n'
265               write(*,*) '       in the Norton law is nonpositive'
266               ier=1
267               return
268            endif
269            if(textpart(4)(1:1).ne.' ') then
270               read(textpart(4)(1:20),'(f20.0)',iostat=istat)
271     &               temperature
272               if(istat.gt.0) then
273                  call inputerror(inpc,ipoinpc,iline,
274     &                 "*CREEP%",ier)
275                  return
276               endif
277            else
278               temperature=0.d0
279            endif
280            elcon(9,ntmat,nmat)=temperature
281         enddo
282!
283         if(ntmat.eq.0) then
284            write(*,*) '*ERROR reading *CREEP: Norton law assumed,'
285            write(*,*) '       yet no constants given'
286            ier=1
287            return
288         endif
289!
290!        interpolating the creep data at the elastic temperature
291!        data points
292!
293         write(*,*) '*INFO: interpolating the creep data at the'
294         write(*,*) '       elastic temperature data points;'
295         write(*,*) '       please note that it is preferable'
296         write(*,*) '       to use exactly the same temperature'
297         write(*,*) '       data points for the elastic and creep'
298         write(*,*) '       data (if not already done so)'
299         write(*,*)
300         write(*,*) 'interpolated creep data'
301         write(*,*) 'temperature    A     n     m'
302!
303         do i=1,nelcon(2,nmat)
304            t1l=elcon(0,i,nmat)
305            call ident2(elcon(9,1,nmat),t1l,ntmat,ncmat_+1,id)
306            if(ntmat.eq.0) then
307               continue
308            elseif((ntmat.eq.1).or.(id.eq.0)) then
309               elcon(3,i,nmat)=elcon(6,1,nmat)
310               elcon(4,i,nmat)=elcon(7,1,nmat)
311               elcon(5,i,nmat)=elcon(8,1,nmat)
312            elseif(id.eq.ntmat) then
313               elcon(3,i,nmat)=elcon(6,id,nmat)
314               elcon(4,i,nmat)=elcon(7,id,nmat)
315               elcon(5,i,nmat)=elcon(8,id,nmat)
316            else
317               do k=3,5
318                  elcon(k,i,nmat)=elcon(k+3,id,nmat)+
319     &               (elcon(k+3,id+1,nmat)-elcon(k+3,id,nmat))*
320     &               (t1l-elcon(9,id,nmat))/
321     &               (elcon(9,id+1,nmat)-elcon(9,id,nmat))
322               enddo
323            endif
324            write(*,*) t1l,(elcon(k,i,nmat),k=3,5)
325         enddo
326         write(*,*)
327!
328      else
329!
330!        elastically anisotropic material with isotropic viscoplasticity
331!        (i.e. isotropic plasticity with Norton creep)
332!
333         do
334            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
335     &           ipoinp,inp,ipoinpc)
336            if((istat.lt.0).or.(key.eq.1)) exit
337            ntmat=ntmat+1
338            if(ntmat.gt.ntmat_) then
339               write(*,*) '*ERROR reading *CREEP: increase ntmat_'
340               ier=1
341               return
342            endif
343            do i=1,3
344               read(textpart(i)(1:20),'(f20.0)',iostat=istat)
345     &               elcon(i+15,ntmat,nmat)
346               if(istat.gt.0) then
347                  call inputerror(inpc,ipoinpc,iline,
348     &                 "*CREEP%",ier)
349                  return
350               endif
351            enddo
352            if(textpart(3)(1:1).ne.' ') then
353               read(textpart(3)(1:20),'(f20.0)',iostat=istat)
354     &                  temperature
355               if(istat.gt.0) then
356                  call inputerror(inpc,ipoinpc,iline,
357     &                 "*CREEP%",ier)
358                  return
359               endif
360            else
361               temperature=0.d0
362            endif
363            elcon(19,ntmat,nmat)=temperature
364         enddo
365!
366!        interpolating the creep data at the elastic temperature
367!        data points
368!
369!        before interpolation: data are stored in positions 16-19:
370!        A,n,m,temperature
371!        after interpolation: data are stored in positions 13-15:
372!        A,n,m
373!
374         write(*,*) '*INFO: interpolating the creep data at the'
375         write(*,*) '       elastic temperature data points;'
376         write(*,*) '       please note that it is preferable'
377         write(*,*) '       to use exactly the same temperature'
378         write(*,*) '       data points for the elastic and creep'
379         write(*,*) '       data (if not already done so)'
380         write(*,*)
381         write(*,*) 'interpolated creep data'
382         write(*,*) 'temperature    A     n     m'
383!
384         if(ntmat.eq.0) then
385            write(*,*) '*ERROR reading *CREEP: Norton law assumed,'
386            write(*,*) '       yet no constants given'
387            ier=1
388            return
389         endif
390!
391         do i=1,nelcon(2,nmat)
392            t1l=elcon(0,i,nmat)
393            call ident2(elcon(19,1,nmat),t1l,ntmat,ncmat_+1,id)
394            if(ntmat.eq.0) then
395               continue
396            elseif((ntmat.eq.1).or.(id.eq.0)) then
397               elcon(13,i,nmat)=elcon(16,1,nmat)
398               elcon(14,i,nmat)=elcon(17,1,nmat)
399               elcon(15,i,nmat)=elcon(18,1,nmat)
400            elseif(id.eq.ntmat) then
401               elcon(13,i,nmat)=elcon(16,id,nmat)
402               elcon(14,i,nmat)=elcon(17,id,nmat)
403               elcon(15,i,nmat)=elcon(18,id,nmat)
404            else
405               do k=13,15
406                  elcon(k,i,nmat)=elcon(k+3,id,nmat)+
407     &               (elcon(k+3,id+1,nmat)-elcon(k+3,id,nmat))*
408     &               (t1l-elcon(19,id,nmat))/
409     &               (elcon(19,id+1,nmat)-elcon(19,id,nmat))
410               enddo
411            endif
412            write(*,*) t1l,(elcon(k,i,nmat),k=13,15)
413         enddo
414         write(*,*)
415!
416      endif
417!
418      return
419      end
420
421