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 springforc_f2f(xl,vl,imat,elcon,nelcon,
20     &  elas,fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,
21     &  plicon,nplicon,npmat_,senergy,nener,cstr,mi,
22     &  springarea,nmethod,ne0,nstate_,xstateini,
23     &  xstate,reltime,ielas,jfaces,igauss,
24     &  pslavsurf,pmastsurf,clearini,venergy,kscale,
25     &  konl,iout,nelem,smscale,mscalmethod)
26!
27!     calculates the force of the spring (face-to-face penalty)
28!
29      implicit none
30!
31      character*8 lakonl
32!
33      integer i,j,k,imat,ncmat_,ntmat_,nope,iflag,mi(*),
34     &  kode,niso,id,nplicon(0:ntmat_,*),npmat_,nelcon(2,*),nener,
35     &  nmethod,ne0,nstate_,ielas,jfaces,kscale,konl(26),
36     &  igauss,nopes,nopem,nopep,iout,nelem,mscalmethod
37!
38      real*8 xl(3,10),elas(21),t1l,al(3),vl(0:mi(2),19),stickslope,
39     &  pl(3,19),xn(3),alpha,beta,fnl(3,19),tp(3),te(3),ftrial(3),
40     &  t(3),dftrial,elcon(0:ncmat_,ntmat_,*),pproj(3),clear,
41     &  xi,et,elconloc(*),plconloc(82),xk,val,xiso(20),yiso(20),
42     &  plicon(0:2*npmat_,ntmat_,*),um,senergy,cstr(6),dg,venergy,
43     &  dfshear,dfnl,springarea(2),overlap,pres,clearini(3,9,*),
44     &  xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),t1(3),t2(3),
45     &  dt1,dte,alnew(3),reltime,weight,xsj2m(3),xs2m(3,7),shp2m(7,9),
46     &  xsj2s(3),xs2s(3,7),shp2s(7,9),pslavsurf(3,*),pmastsurf(6,*),
47     &  smscale(*)
48!
49      include "gauss.f"
50!
51      iflag=2
52!
53      if(nener.eq.1) venergy=0.d0
54!
55!     # of master nodes
56!
57      nopem=ichar(lakonl(8:8))-48
58!
59!     # of slave nodes
60!
61      nopes=nope-nopem
62!
63!     actual positions of the master nodes belonging to the contact spring
64!
65      do i=1,nopem
66         do j=1,3
67            pl(j,i)=xl(j,i)+vl(j,i)
68         enddo
69      enddo
70!
71!     actual positions of the slave nodes belonging to the contact spring
72!
73      if(nmethod.ne.2) then
74         do i=nopem+1,nope
75            do j=1,3
76               pl(j,i)=xl(j,i)+clearini(j,i-nopem,jfaces)*reltime
77     &                        +vl(j,i)
78            enddo
79         enddo
80      else
81!
82!        for frequency calculations the eigenmodes are freely
83!        scalable, leading to problems with contact finding
84!
85         do i=nopem+1,nope
86            do j=1,3
87               pl(j,i)=xl(j,i)+clearini(j,i-nopem,jfaces)*reltime
88            enddo
89         enddo
90      endif
91!
92!     location of integration point in slave face
93!
94      xi=pslavsurf(1,igauss)
95      et=pslavsurf(2,igauss)
96      weight=pslavsurf(3,igauss)
97!
98      iflag=1
99      if(nopes.eq.8) then
100          call shape8q(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
101      elseif(nopes.eq.4) then
102          call shape4q(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
103      elseif(nopes.eq.6) then
104          call shape6tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
105      else
106          call shape3tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
107      endif
108!
109      nopep=nope+1
110!
111      do k=1,3
112          pl(k,nopep)=0.d0
113          vl(k,nopep)=0.d0
114          do j=1,nopes
115              pl(k,nopep)=pl(k,nopep)+shp2s(4,j)*pl(k,nopem+j)
116              vl(k,nopep)=vl(k,nopep)+shp2s(4,j)*vl(k,nopem+j)
117          enddo
118      enddo
119!
120!     corresponding location in the master face
121!
122      xi=pmastsurf(1,igauss)
123      et=pmastsurf(2,igauss)
124!
125!     determining the jacobian vector on the surface
126!
127      iflag=2
128      if(nopem.eq.8) then
129         call shape8q(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
130      elseif(nopem.eq.4) then
131         call shape4q(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
132      elseif(nopem.eq.6) then
133         call shape6tri(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
134      else
135         call shape3tri(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
136      endif
137!
138      do i=1,3
139         pproj(i)=0.d0
140         do j=1,nopem
141            pproj(i)=pproj(i)+shp2m(4,j)*pl(i,j)
142         enddo
143      enddo
144!
145!     distance vector between both
146!
147      do i=1,3
148         al(i)=pl(i,nopep)-pproj(i)
149      enddo
150!
151!     normal on the master face
152!
153      xn(1)=pmastsurf(4,igauss)
154      xn(2)=pmastsurf(5,igauss)
155      xn(3)=pmastsurf(6,igauss)
156!
157!     distance from surface along normal (= clearance)
158!
159      clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
160!
161!     check for a reduction of the initial penetration, if any
162!
163      if(nmethod.eq.1) then
164         clear=clear-springarea(2)*(1.d0-reltime)
165      endif
166      cstr(1)=clear
167!
168!
169      if(int(elcon(3,1,imat)).eq.1) then
170!
171!        exponential overclosure
172!
173         if(dabs(elcon(2,1,imat)).lt.1.d-30) then
174            elas(1)=0.d0
175            beta=1.d0
176         else
177!
178            alpha=elcon(2,1,imat)*springarea(1)
179            beta=elcon(1,1,imat)
180            if(-beta*clear.gt.23.d0-dlog(alpha)) then
181               beta=(dlog(alpha)-23.d0)/clear
182            endif
183            elas(1)=dexp(-beta*clear+dlog(alpha))
184         endif
185         if(nener.eq.1) then
186            senergy=elas(1)/beta;
187         endif
188      elseif((int(elcon(3,1,imat)).eq.2).or.
189     &       (int(elcon(3,1,imat)).eq.4)) then
190!
191!        linear overclosure/tied overclosure
192!
193         elas(1)=-springarea(1)*elcon(2,1,imat)*clear/kscale
194!
195!     spring scaling for explicit dynamics
196!
197         if((mscalmethod.eq.2).or.(mscalmethod.eq.3)) then
198            elas(1)=elas(1)*smscale(nelem)
199         endif
200!
201         if(nener.eq.1) then
202            senergy=-elas(1)*clear/2.d0;
203         endif
204      elseif(int(elcon(3,1,imat)).eq.3) then
205!
206!        tabular overclosure
207!
208!        interpolating the material data
209!
210         call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
211     &     elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
212         overlap=-clear
213         niso=int(plconloc(81))
214         do i=1,niso
215            xiso(i)=plconloc(2*i-1)
216            yiso(i)=plconloc(2*i)
217         enddo
218         call ident(xiso,overlap,niso,id)
219         if(id.eq.0) then
220            pres=yiso(1)
221            if(nener.eq.1) then
222               senergy=yiso(1)*overlap;
223            endif
224         elseif(id.eq.niso) then
225            pres=yiso(niso)
226            if(nener.eq.1) then
227               senergy=yiso(1)*xiso(1)
228               do i=2,niso
229                  senergy=senergy+(xiso(i)-xiso(i-1))*
230     &                            (yiso(i)+yiso(i-1))/2.d0
231               enddo
232               senergy=senergy+(pres-xiso(niso))*yiso(niso)
233            endif
234         else
235            xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
236            pres=yiso(id)+xk*(overlap-xiso(id))
237            if(nener.eq.1) then
238               senergy=yiso(1)*xiso(1)
239               do i=2, id
240                  senergy=senergy+(xiso(i)-xiso(i-1))*
241     &                 (yiso(i)+yiso(i-1))/2.d0
242               enddo
243               senergy=senergy+(overlap-xiso(id))*(pres+yiso(id))/2.d0
244            endif
245         endif
246         elas(1)=springarea(1)*pres
247         if(nener.eq.1) senergy=springarea(1)*senergy
248      endif
249!
250!     forces in the nodes of the contact element
251!
252      do i=1,3
253         fnl(i,nopep)=-elas(1)*xn(i)
254      enddo
255      cstr(4)=elas(1)/springarea(1)
256!
257!     Coulomb friction for static calculations
258!
259      if((int(elcon(3,1,imat)).eq.4).and.(ncmat_.lt.7)) then
260         write(*,*) '*ERROR in springforc_f2f: for contact'
261         write(*,*) '       with pressure-overclosure=tied'
262         write(*,*) '       a stick slope using the *FRICTION'
263         write(*,*) '       card is mandatory'
264         call exit(201)
265      endif
266!
267      if(ncmat_.ge.7) then
268!
269!        tied contact
270!
271         if(int(elcon(3,1,imat)).eq.4) then
272            um=1.d30
273         else
274            um=elcon(6,1,imat)
275         endif
276         stickslope=elcon(7,1,imat)/kscale
277!
278!     spring scaling for explicit dynamics
279!
280         if((mscalmethod.eq.2).or.(mscalmethod.eq.3)) then
281            stickslope=stickslope*smscale(nelem)
282         endif
283!
284         if(um.gt.0.d0) then
285            if(1.d0 - dabs(xn(1)).lt.1.5231d-6) then
286!
287!     calculating the local directions on master surface
288!
289               t1(1)=-xn(3)*xn(1)
290               t1(2)=-xn(3)*xn(2)
291               t1(3)=1.d0-xn(3)*xn(3)
292            else
293               t1(1)=1.d0-xn(1)*xn(1)
294               t1(2)=-xn(1)*xn(2)
295               t1(3)=-xn(1)*xn(3)
296            endif
297            dt1=dsqrt(t1(1)*t1(1)+t1(2)*t1(2)+t1(3)*t1(3))
298            do i=1,3
299               t1(i)=t1(i)/dt1
300            enddo
301            t2(1)=xn(2)*t1(3)-xn(3)*t1(2)
302            t2(2)=xn(3)*t1(1)-xn(1)*t1(3)
303            t2(3)=xn(1)*t1(2)-xn(2)*t1(1)
304!
305!     linear stiffness of the shear stress versus
306!     slip curve
307!
308            xk=stickslope*springarea(1)
309!
310!     calculating the relative displacement between the slave node
311!     and its projection on the master surface
312!
313            do i=1,3
314               alnew(i)=vl(i,nopep)
315               do j=1,nopem
316                  alnew(i)=alnew(i)-shp2m(4,j)*vl(i,j)
317               enddo
318            enddo
319!
320!     calculating the difference in relative displacement since
321!     the start of the increment = lamda^*
322!
323            do i=1,3
324               al(i)=alnew(i)-xstateini(3+i,1,ne0+igauss)
325            enddo
326!
327!     ||lambda^*||
328!
329            val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
330!
331!     update the relative tangential displacement
332!
333            do i=1,3
334               t(i)=xstateini(6+i,1,ne0+igauss)+al(i)-val*xn(i)
335            enddo
336!
337!     store the actual relative displacement and
338!     the actual relative tangential displacement
339!
340            do i=1,3
341               xstate(3+i,1,ne0+igauss)=alnew(i)
342               xstate(6+i,1,ne0+igauss)=t(i)
343            enddo
344!
345!     size of normal force
346!
347            dfnl=dsqrt(fnl(1,nopep)**2+fnl(2,nopep)**2+
348     &           fnl(3,nopep)**2)
349!
350!     maximum size of shear force
351!
352            if(int(elcon(3,1,imat)).eq.4) then
353               dfshear=1.d30
354            else
355               dfshear=um*dfnl
356            endif
357!
358!     plastic and elastic slip
359!
360            do i=1,3
361               tp(i)=xstateini(i,1,ne0+igauss)
362               te(i)=t(i)-tp(i)
363            enddo
364            dte=dsqrt(te(1)*te(1)+te(2)*te(2)+te(3)*te(3))
365!
366!     trial force
367!
368            do i=1,3
369               ftrial(i)=xk*te(i)
370            enddo
371            dftrial=xk*dte
372c            dftrial=dsqrt(ftrial(1)**2+ftrial(2)**2+ftrial(3)**2)
373!
374!     check whether stick or slip
375!
376            if((dftrial.lt.dfshear).or.(dftrial.le.0.d0).or.
377     &           (ielas.eq.1)) then
378!
379!     stick
380!
381               do i=1,3
382                  fnl(i,nopep)=fnl(i,nopep)+ftrial(i)
383                  xstate(i,1,ne0+igauss)=tp(i)
384               enddo
385               cstr(5)=(ftrial(1)*t1(1)+ftrial(2)*t1(2)+
386     &              ftrial(3)*t1(3))/springarea(1)
387               cstr(6)=(ftrial(1)*t2(1)+ftrial(2)*t2(2)+
388     &              ftrial(3)*t2(3))/springarea(1)
389!
390!              shear elastic energy
391!
392               if(nener.eq.1) senergy=senergy+dftrial*dte
393            else
394!
395!     slip
396!
397               dg=(dftrial-dfshear)/xk
398               do i=1,3
399                  ftrial(i)=te(i)/dte
400                  fnl(i,nopep)=fnl(i,nopep)+dfshear*ftrial(i)
401                  xstate(i,1,ne0+igauss)=tp(i)+dg*ftrial(i)
402               enddo
403               cstr(5)=(dfshear*ftrial(1)*t1(1)+
404     &              dfshear*ftrial(2)*t1(2)+
405     &              dfshear*ftrial(3)*t1(3))/springarea(1)
406               cstr(6)=(dfshear*ftrial(1)*t2(1)+
407     &              dfshear*ftrial(2)*t2(2)+
408     &              dfshear*ftrial(3)*t2(3))/springarea(1)
409!
410!              shear elastic and viscous energy
411!
412               if(nener.eq.1) then
413                  senergy=senergy+dfshear*dfshear/xk
414                  venergy=dg*dfshear
415               endif
416!
417            endif
418         endif
419!
420!     storing the tangential displacements
421!
422         cstr(2)=t(1)*t1(1)+t(2)*t1(2)+t(3)*t1(3)
423         cstr(3)=t(1)*t2(1)+t(2)*t2(2)+t(3)*t2(3)
424      endif
425!
426!     force in the master nodes
427!
428      do i=1,3
429         do j=1,nopem
430            fnl(i,j)=-shp2m(4,j)*fnl(i,nopep)
431         enddo
432      enddo
433!
434!     forces in the nodes of the slave face
435!
436      do i=1,3
437          do j=1,nopes
438              fnl(i,nopem+j)=shp2s(4,j)*fnl(i,nopep)
439          enddo
440      enddo
441!
442!     write statements for Malte Krack
443!
444c      if(iout.gt.0) then
445c         write(*,*) 'contact element: ',nelem
446c         write(*,*) 'undeformed location of the integration point'
447c         write(*,*) ((pl(j,nopep)-vl(j,nopep)),j=1,3)
448c         write(*,*) 'deformed location of the integration point'
449c         write(*,*) (pl(j,nopep),j=1,3)
450c         write(*,*) 'nodes and shape values'
451c         do j=1,nopem
452c            write(*,*) konl(j),-shp2m(4,j)
453c         enddo
454c         do j=1,nopes
455c            write(*,*) konl(nopem+j),shp2s(4,j)
456c         enddo
457c         write(*,*) 'contact force'
458c         write(*,*) fnl(1,nopep),fnl(2,nopep),fnl(3,nopep)
459c         write(*,*) 'slave area'
460c         write(*,*) springarea(1)
461c      endif
462!
463      return
464      end
465
466