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 calcstabletimeinccont(ne,lakon,kon,ipkon,mi,
20     &  ielmat,elcon,mortar,adb,alpha,nactdof,springarea,
21     &  ne0,ntmat_,ncmat_,dtcont,smscale,dtset,mscalmethod)
22!
23!     Calculates the critical time increment (CTI) for contact
24!     spring elements based on the Courant
25!     Criterion for Explicit Dynamics calculations.
26!
27!     mscalmethod<0: no explicit dynamics
28!                 0: explicit dynamics, no scaling
29!                 1: explicit dynamics, volumetric scaling active
30!                 2: explicit dynamics, contact scaling active
31!                 3: explicit dynamics, volumetric and contact scaling
32!
33      implicit none
34!
35      character*8 lakon(*),lakonl
36!
37      integer j,ne,nope,kon(*),ipkon(*),indexe,nelem,
38     &  ncmat_,ntmat_,mi(*),ielmat(mi(3),*),imat,mortar,
39     &  indexn,nactdof(0:mi(2),*),ne0,nopem,mscalmethod,
40     &  icount
41!
42      real*8 elcon(0:ncmat_,ntmat_,*),safefac,xk,adb(*),alpha(*),
43     &  critom,damping,springms,springmm,springfac,dtcont,xmacont,
44     &  springarea(2,*),areaslav,bet,gam,smscale(*),dtset
45!
46      dtcont=1.d30
47      safefac=0.80d0
48!
49      xmacont=0.0d0
50!
51      damping=0.d0
52      icount=0
53!
54      bet=(1.d0-alpha(1))*(1.d0-alpha(1))/4.d0
55      gam=0.5d0-alpha(1)
56!
57!     Omega Critical
58!     Om_cr=dt*freq_max
59!
60      critom=dsqrt(damping*damping*(1.d0+2.d0*alpha(1)*(1.d0-gam))
61     &     *(1.d0+2.d0*alpha(1)*(1.d0-gam))
62     &    +   2.d0*(gam+2.d0*alpha(1)*(gam-bet)))
63      critom=0.98d0*(-damping*(1.d0+2.d0*alpha(1)*(1.d0-gam))+critom)
64     &     /(gam+2.d0*alpha(1)*(gam-bet)) !eq 25 miranda
65!
66!     ** DO per element
67!
68      do nelem=ne0+1,ne
69         indexe=ipkon(nelem)
70         lakonl=lakon(nelem)
71         imat=ielmat(1,nelem)
72!
73               xk=elcon(2,1,imat)
74!
75               springmm=0.0d0
76               springms=0.0d0
77!
78               if(mortar.eq.0) then
79!
80!                 node-to-face
81!
82!                 nope is the total number of nodes:
83!                 master nodes+1 slave node
84!                 (notice -47 instead of -48)
85!
86                  nope=ichar(lakonl(8:8))-47
87                  springfac=1.d0
88!
89                  do indexn=1,nope
90                     xmacont=0.0d0
91!
92                     do j=1,3
93                       if(nactdof(j,kon(indexe+indexn)).gt.0)then
94                           xmacont=max(xmacont,
95     &                         adb(nactdof(j,kon(indexe+indexn))))
96                        endif
97                     enddo
98!
99                     if(indexn.eq.nope)then
100                        springms=springms+xmacont
101!
102!                       mean mass at master nodes
103!
104                        springmm=springmm/(nope-1.0d0)
105                     else
106                        springmm=springmm+xmacont
107                     endif
108                  enddo
109!
110                  areaslav=springarea(1,kon(indexe+nope+1))
111                  xk=xk*areaslav
112!
113!                 checking whether mass or slave side is constrained
114!
115                  if((springmm.le.0.d0).and.(springms.le.0.d0)) then
116                     cycle
117                  elseif(springmm.le.0.d0) then
118                     springmm=springms
119                  elseif(springms.le.0.d0) then
120                     springms=springmm
121                  endif
122!
123               elseif(mortar.eq.1) then
124!
125!                 face-to-face
126!
127                  nopem=ichar(lakonl(8:8))-48
128                  springfac=0.1d0
129!
130                  do indexn=1,kon(indexe)
131!
132                     xmacont=0.0d0
133                     do j=1,3
134                        if(nactdof(j,kon(indexe+indexn)).gt.0)then
135                           xmacont=max(xmacont,
136     &                          adb(nactdof(j,kon(indexe+indexn))))
137                        endif
138                     enddo
139                     if(indexn.gt.nopem)then !slave
140                        springms=springms+xmacont
141                     else !master
142                        springmm=springmm+xmacont
143                     endif
144                  enddo
145!
146!           mean mass at slave and master nodes
147!
148                  springms=springms/(kon(indexe)-nopem)
149                  springmm=springmm/nopem
150!
151                  areaslav=springarea(1,kon(1+indexe+kon(indexe)))
152                  xk=xk*areaslav
153!
154!                 checking whether mass or slave side is constrained
155!
156                  if((springmm.le.0.d0).and.(springms.le.0.d0)) then
157                     cycle
158                  elseif(springmm.le.0.d0) then
159                     springmm=springms
160                  elseif(springms.le.0.d0) then
161                     springms=springmm
162                  endif
163               endif
164!
165!              linear pressure-overclosure
166!
167               if(int(elcon(3,1,imat)).eq.2) then
168!
169                  springmm=springmm/2.0d0
170                  springms=springms/2.0d0
171!
172!                 time increment needed by contact spring
173!
174                  smscale(nelem)=springfac*critom*
175     &                 dsqrt(springmm*springms/
176     &                       ((springmm+springms)*xk))
177!
178                  dtcont =
179     &                 min(dtcont,smscale(nelem))
180!
181               else
182                  write(*,*) '*ERROR in calcstabletimeinccont:'
183                  write(*,*) '       in explicit dynamic calculations'
184                  write(*,*) '       only linear pressure-overclosure'
185                  write(*,*) '       is allowed'
186                  call exit(201)
187               endif
188      enddo
189!     ** ENDDO per element
190!
191!     check whether spring scaling is necessary in order to meet the
192!     smallest allowed time increment specified by the user
193!
194      if(dtcont.lt.dtset/safefac) then
195         dtset=dtset/safefac
196         if(mscalmethod.eq.1) then
197            mscalmethod=3
198         else
199            mscalmethod=2
200         endif
201!
202         do nelem=ne0+1,ne
203            if(smscale(nelem).ge.dtset) then
204               smscale(nelem)=1.d0
205            else
206               smscale(nelem)=(smscale(nelem)/dtset)**2
207               icount=icount+1
208            endif
209         enddo
210!
211         write(*,*) 'Selective Spring Scaling is active'
212         write(*,*)
213         write(*,*) 'Scaling factor of time increment:',dtset/dtcont
214         write(*,*) 'Reduction of spring stiffness by (maximum) =',
215     &        (dtset/dtcont)**2
216         write(*,*) 'In total ',icount,'elements were scaled of'
217         write(*,*) '         ',ne-ne0,' contact elements'
218         write(*,*)
219!
220         dtset=dtset*safefac
221!
222      else
223         if(mscalmethod.eq.3) then
224            mscalmethod=1
225         elseif(mscalmethod.eq.2) then
226            mscalmethod=0
227         endif
228      endif
229!
230      dtcont=dtcont*safefac
231!
232      return
233      end
234