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 springstiff_n2f(xl,elas,konl,voldl,s,imat,elcon,nelcon,
20     &  ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,plicon,
21     &  nplicon,npmat_,iperturb,springarea,nmethod,mi,ne0,
22     &  nstate_,xstateini,xstate,reltime,nasym,ielorien,orab,
23     &  norien,nelem)
24!
25!     calculates the stiffness of a spring (node-to-face penalty)
26!     (User's manual -> theory -> boundary conditions ->
27!      node-to-face penalty contact)
28!
29      implicit none
30!
31      character*8 lakonl
32!
33      integer konl(20),i,j,imat,ncmat_,ntmat_,k,l,nope,nterms,iflag,
34     &  i1,kode,niso,id,nplicon(0:ntmat_,*),npmat_,nelcon(2,*),
35     &  iperturb(*),nmethod,mi(*),ne0,nstate_,nasym,ielorien(mi(3),*),
36     &  norien,nelem,idof,idof1,idof2,iorien
37!
38      real*8 xl(3,10),elas(21),ratio(9),pproj(3),val,shp2(7,9),
39     &  al(3),s(60,60),voldl(0:mi(2),10),pl(3,10),xn(3),dm,
40     &  c1,c2,c3,c4,alpha,beta,elcon(0:ncmat_,ntmat_,*),xm(3),
41     &  xmu(3,3,10),dxmu(3,10),dval(3,10),fpu(3,3,10),xi,et,
42     &  xs2(3,7),t1l,elconloc(21),plconloc(802),xk,fk,
43     &  xiso(200),yiso(200),dd0,plicon(0:2*npmat_,ntmat_,*),
44     &  a11,a12,a22,b1(3,10),b2(3,10),dal(3,3,10),qxxy(3),fnl(3),
45     &  qxyy(3),dxi(3,10),det(3,10),determinant,c11,c12,c22,
46     &  qxyx(3),qyxy(3),springarea(2),dd,dist,t(3),tu(3,3,10),
47     &  xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
48     &  um,eps,pi,dftdt(3,3),tp(3),te(3),ftrial(3),clear,
49     &  dftrial,dfnl,dfshear,dg,dte,alnew(3),dfn(3,10),reltime,
50     &  overlap,pres,dpresdoverlap,overclosure,orab(7,*),
51     &  xn1(3),xn2(3),a(3,3)
52!
53!
54!
55      iflag=4
56!
57!     actual positions of the nodes belonging to the contact spring
58!     (otherwise no contact force)
59!
60      do i=1,nope
61         do j=1,3
62            pl(j,i)=xl(j,i)+voldl(j,i)
63         enddo
64      enddo
65!
66      if(lakonl(7:7).eq.'A') then
67         if(lakonl(4:6).eq.'RNG') then
68!
69!           SPRINGA-element
70!
71            dd0=dsqrt((xl(1,2)-xl(1,1))**2
72     &           +(xl(2,2)-xl(2,1))**2
73     &           +(xl(3,2)-xl(3,1))**2)
74            dd=dsqrt((pl(1,2)-pl(1,1))**2
75     &           +(pl(2,2)-pl(2,1))**2
76     &           +(pl(3,2)-pl(3,1))**2)
77            if(dd.le.0.d0) then
78               write(*,*) '*ERROR in springstiff_n2f: spring connects'
79               write(*,*) '       two nodes at the same location:'
80               write(*,*) '       spring length is zero'
81               call exit(201)
82            endif
83            do i=1,3
84               xn(i)=(pl(i,2)-pl(i,1))/dd
85            enddo
86            val=dd-dd0
87!
88!     interpolating the material data
89!
90            call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
91     &           elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
92!
93!     calculating the spring force and the spring constant
94!
95            if(kode.eq.2)then
96               xk=elconloc(1)
97               fk=xk*val
98            else
99               niso=int(plconloc(801))
100               do i=1,niso
101                  xiso(i)=plconloc(2*i-1)
102                  yiso(i)=plconloc(2*i)
103               enddo
104               call ident(xiso,val,niso,id)
105               if(id.eq.0) then
106                  xk=0.d0
107                  fk=yiso(1)
108               elseif(id.eq.niso) then
109                  xk=0.d0
110                  fk=yiso(niso)
111               else
112                  xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
113                  fk=yiso(id)+xk*(val-xiso(id))
114               endif
115            endif
116!
117            c1=fk/dd
118            c2=xk-c1
119            do i=1,3
120               do j=1,3
121                  s(i,j)=c2*xn(i)*xn(j)
122               enddo
123               s(i,i)=s(i,i)+c1
124            enddo
125         else
126!
127!           GAP-element
128!           interpolating the material data
129!
130            call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
131     &           elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
132!
133            dd=elconloc(1)
134            xn(1)=elconloc(2)
135            xn(2)=elconloc(3)
136            xn(3)=elconloc(4)
137            xk=elconloc(5)
138            pi=4.d0*datan(1.d0)
139            eps=pi*elconloc(6)/xk
140            overclosure=-dd-xn(1)*(voldl(1,2)-voldl(1,1))
141     &                     -xn(2)*(voldl(2,2)-voldl(2,1))
142     &                     -xn(3)*(voldl(3,2)-voldl(3,1))
143            fk=-xk*overclosure*(0.5d0+datan(overclosure/eps)/pi)
144            c2=xk*((0.5d0+datan(overclosure/eps)/pi)
145     &        +overclosure/(pi*eps*(1.d0+(overclosure/eps)**2)))
146            do i=1,3
147               do j=1,3
148                  s(i,j)=c2*xn(i)*xn(j)
149               enddo
150            enddo
151         endif
152!
153         do i=1,3
154            do j=1,3
155               s(i+3,j)=-s(i,j)
156               s(i,j+3)=-s(i,j)
157               s(i+3,j+3)=s(i,j)
158            enddo
159         enddo
160         return
161      elseif(lakonl(7:7).eq.'1') then
162!
163!        spring1-element
164!        determine the direction of action of the spring
165!
166         idof=nint(elcon(3,1,imat))
167!
168         if(norien.gt.0) then
169            iorien=max(0,ielorien(1,nelem))
170         else
171            iorien=0
172         endif
173!
174         if(iorien.eq.0) then
175            do i=1,3
176               xn(i)=0.d0
177            enddo
178            xn(idof)=1.d0
179         else
180            call transformatrix(orab(1,iorien),xl(1,1),a)
181            do i=1,3
182               xn(i)=a(i,idof)
183            enddo
184         endif
185!
186!     interpolating the material data
187!
188         call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
189     &           elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
190!
191!     calculating the spring constant
192!
193         if(kode.eq.2)then
194            xk=elconloc(1)
195         else
196            niso=int(plconloc(801))
197            do i=1,niso
198               xiso(i)=plconloc(2*i-1)
199               yiso(i)=plconloc(2*i)
200            enddo
201            call ident(xiso,val,niso,id)
202            if(id.eq.0) then
203               xk=0.d0
204            elseif(id.eq.niso) then
205               xk=0.d0
206            else
207               xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
208            endif
209         endif
210!
211!        assembling the stiffness matrix
212!
213         do i=1,3
214            do j=1,3
215               s(i,j)=xk*xn(i)*xn(j)
216            enddo
217         enddo
218         return
219      elseif(lakonl(7:7).eq.'2') then
220!
221!        spring2-element
222!        determine the direction of action of the spring
223!
224         idof1=nint(elcon(3,1,imat))
225         idof2=nint(elcon(4,1,imat))
226!
227         if(norien.gt.0) then
228            iorien=max(0,ielorien(1,nelem))
229         else
230            iorien=0
231         endif
232!
233         if(iorien.eq.0) then
234            do i=1,3
235               xn1(i)=0.d0
236               xn2(i)=0.d0
237            enddo
238            xn1(idof1)=1.d0
239            xn2(idof2)=1.d0
240         else
241            call transformatrix(orab(1,iorien),xl(1,1),a)
242            do i=1,3
243               xn1(i)=a(i,idof1)
244            enddo
245            call transformatrix(orab(1,iorien),xl(1,2),a)
246            do i=1,3
247               xn2(i)=a(i,idof2)
248            enddo
249         endif
250!
251!     interpolating the material data
252!
253         call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
254     &           elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
255!
256!        calculating the spring constant
257!
258         if(kode.eq.2)then
259            xk=elconloc(1)
260         else
261            niso=int(plconloc(801))
262            do i=1,niso
263               xiso(i)=plconloc(2*i-1)
264               yiso(i)=plconloc(2*i)
265            enddo
266            call ident(xiso,val,niso,id)
267            if(id.eq.0) then
268               xk=0.d0
269            elseif(id.eq.niso) then
270               xk=0.d0
271            else
272               xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
273            endif
274         endif
275!
276!        assembling the stiffness matrix
277!
278         do i=1,3
279            do j=1,3
280               s(i,j)=xk*xn1(i)*xn1(j)
281               s(i,j+3)=-xk*xn1(i)*xn2(j)
282               s(i+3,j)=-xk*xn2(i)*xn1(j)
283               s(i+3,j+3)=xk*xn2(i)*xn2(j)
284            enddo
285         enddo
286         return
287!
288      endif
289!
290!     contact springs
291!
292      nterms=nope-1
293!
294!     vector al connects the actual position of the slave node
295!     with its projection on the master face = vec_r (User's
296!     manual -> theory -> boundary conditions -> node-to-face
297!     penalty contact)
298!
299      do i=1,3
300         pproj(i)=pl(i,nope)
301      enddo
302      call attach_2d(pl,pproj,nterms,ratio,dist,xi,et)
303      do i=1,3
304         al(i)=pl(i,nope)-pproj(i)
305      enddo
306!
307!     determining the jacobian vector on the surface
308!
309      if(nterms.eq.8) then
310         call shape8q(xi,et,pl,xm,xs2,shp2,iflag)
311      elseif(nterms.eq.4) then
312         call shape4q(xi,et,pl,xm,xs2,shp2,iflag)
313      elseif(nterms.eq.6) then
314         call shape6tri(xi,et,pl,xm,xs2,shp2,iflag)
315      else
316         call shape3tri(xi,et,pl,xm,xs2,shp2,iflag)
317      endif
318!
319!     d xi / d vec_u_j
320!     d eta / d vec_u_j
321!
322!     dxi and det are determined from the orthogonality
323!     condition
324!
325      a11=-(xs2(1,1)*xs2(1,1)+xs2(2,1)*xs2(2,1)+xs2(3,1)*xs2(3,1))
326     &    +al(1)*xs2(1,5)+al(2)*xs2(2,5)+al(3)*xs2(3,5)
327      a12=-(xs2(1,1)*xs2(1,2)+xs2(2,1)*xs2(2,2)+xs2(3,1)*xs2(3,2))
328     &    +al(1)*xs2(1,6)+al(2)*xs2(2,6)+al(3)*xs2(3,6)
329      a22=-(xs2(1,2)*xs2(1,2)+xs2(2,2)*xs2(2,2)+xs2(3,2)*xs2(3,2))
330     &    +al(1)*xs2(1,7)+al(2)*xs2(2,7)+al(3)*xs2(3,7)
331!
332      do i=1,3
333         do j=1,nterms
334            b1(i,j)=shp2(4,j)*xs2(i,1)-shp2(1,j)*al(i)
335            b2(i,j)=shp2(4,j)*xs2(i,2)-shp2(2,j)*al(i)
336         enddo
337         b1(i,nope)=-xs2(i,1)
338         b2(i,nope)=-xs2(i,2)
339      enddo
340!
341      determinant=a11*a22-a12*a12
342      c11=a22/determinant
343      c12=-a12/determinant
344      c22=a11/determinant
345!
346      do i=1,3
347         do j=1,nope
348            dxi(i,j)=c11*b1(i,j)+c12*b2(i,j)
349            det(i,j)=c12*b1(i,j)+c22*b2(i,j)
350         enddo
351      enddo
352!
353!     d vec_r / d vec_u_k
354!
355      do i=1,nope
356         do j=1,3
357            do k=1,3
358               dal(j,k,i)=-xs2(j,1)*dxi(k,i)-xs2(j,2)*det(k,i)
359            enddo
360         enddo
361      enddo
362      do i=1,nterms
363         do j=1,3
364            dal(j,j,i)=dal(j,j,i)-shp2(4,i)
365         enddo
366      enddo
367      do j=1,3
368         dal(j,j,nope)=dal(j,j,nope)+1.d0
369      enddo
370!
371!     d2q/dxx x dq/dy
372!
373      qxxy(1)=xs2(2,5)*xs2(3,2)-xs2(3,5)*xs2(2,2)
374      qxxy(2)=xs2(3,5)*xs2(1,2)-xs2(1,5)*xs2(3,2)
375      qxxy(3)=xs2(1,5)*xs2(2,2)-xs2(2,5)*xs2(1,2)
376!
377!     dq/dx x d2q/dyy
378!
379      qxyy(1)=xs2(2,1)*xs2(3,7)-xs2(3,1)*xs2(2,7)
380      qxyy(2)=xs2(3,1)*xs2(1,7)-xs2(1,1)*xs2(3,7)
381      qxyy(3)=xs2(1,1)*xs2(2,7)-xs2(2,1)*xs2(1,7)
382!
383!     Modified by Stefan Sicklinger
384!
385!     dq/dx x d2q/dxy
386!
387      qxyx(1)=xs2(2,1)*xs2(3,6)-xs2(3,1)*xs2(2,6)
388      qxyx(2)=xs2(3,1)*xs2(1,6)-xs2(1,1)*xs2(3,6)
389      qxyx(3)=xs2(1,1)*xs2(2,6)-xs2(2,1)*xs2(1,6)
390!
391!     d2q/dxy x dq/dy
392!
393      qyxy(1)=xs2(2,6)*xs2(3,2)-xs2(3,6)*xs2(2,2)
394      qyxy(2)=xs2(3,6)*xs2(1,2)-xs2(1,6)*xs2(3,2)
395      qyxy(3)=xs2(1,6)*xs2(2,2)-xs2(2,6)*xs2(1,2)
396!
397!
398!     End modifications
399!
400!     normal on the surface
401!
402      dm=dsqrt(xm(1)*xm(1)+xm(2)*xm(2)+xm(3)*xm(3))
403      do i=1,3
404         xn(i)=xm(i)/dm
405      enddo
406!
407!     distance from surface along normal (= clearance)
408!
409      clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
410      if(nmethod.eq.1) then
411         clear=clear-springarea(2)*(1.d0-reltime)
412      endif
413!
414!     representative area: usually the slave surface stored in
415!     springarea; however, if no area was assigned because the
416!     node does not belong to any element, the master surface
417!     is used
418!
419      if(springarea(1).le.0.d0) then
420         if(nterms.eq.3) then
421            springarea(1)=dm/2.d0
422         else
423            springarea(1)=dm*4.d0
424         endif
425      endif
426!
427!     alpha and beta, taking the representative area into account
428!     (conversion of pressure into force)
429!
430      if(int(elcon(3,1,imat)).eq.1) then
431!
432!        exponential overclosure
433!
434         if(dabs(elcon(2,1,imat)).lt.1.d-30) then
435            elas(1)=0.d0
436            elas(2)=0.d0
437         else
438            alpha=elcon(2,1,imat)*springarea(1)
439            beta=elcon(1,1,imat)
440            if(-beta*clear.gt.23.d0-dlog(alpha)) then
441               beta=(dlog(alpha)-23.d0)/clear
442            endif
443            elas(1)=dexp(-beta*clear+dlog(alpha))
444            elas(2)=-beta*elas(1)
445         endif
446      elseif(int(elcon(3,1,imat)).eq.2) then
447!
448!        linear overclosure
449!
450         pi=4.d0*datan(1.d0)
451         eps=elcon(1,1,imat)*pi/elcon(2,1,imat)
452         elas(1)=(-springarea(1)*elcon(2,1,imat)*clear*
453     &            (0.5d0+datan(-clear/eps)/pi))
454         elas(2)=-springarea(1)*elcon(2,1,imat)*
455     &            ((0.5d0+datan(-clear/eps)/pi)-
456     &             clear/(pi*eps*(1.d0+(clear/eps)**2)))
457      elseif(int(elcon(3,1,imat)).eq.3) then
458!
459!        tabular overclosure
460!
461!        interpolating the material data
462!
463         call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
464     &     elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
465         overlap=-clear
466         niso=int(plconloc(801))
467         do i=1,niso
468            xiso(i)=plconloc(2*i-1)
469            yiso(i)=plconloc(2*i)
470         enddo
471         call ident(xiso,overlap,niso,id)
472         if(id.eq.0) then
473            dpresdoverlap=0.d0
474            pres=yiso(1)
475         elseif(id.eq.niso) then
476            dpresdoverlap=0.d0
477            pres=yiso(niso)
478         else
479            dpresdoverlap=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
480            pres=yiso(id)+dpresdoverlap*(overlap-xiso(id))
481         endif
482         elas(1)=springarea(1)*pres
483         elas(2)=-springarea(1)*dpresdoverlap
484      endif
485!
486!     contact force
487!
488      do i=1,3
489         fnl(i)=-elas(1)*xn(i)
490      enddo
491!
492!     derivatives of the jacobian vector w.r.t. the displacement
493!     vectors (d m / d u_k)
494!
495      do k=1,nterms
496         xmu(1,1,k)=0.d0
497         xmu(2,2,k)=0.d0
498         xmu(3,3,k)=0.d0
499         xmu(1,2,k)=shp2(1,k)*xs2(3,2)-shp2(2,k)*xs2(3,1)
500         xmu(2,3,k)=shp2(1,k)*xs2(1,2)-shp2(2,k)*xs2(1,1)
501         xmu(3,1,k)=shp2(1,k)*xs2(2,2)-shp2(2,k)*xs2(2,1)
502         xmu(1,3,k)=-xmu(3,1,k)
503         xmu(2,1,k)=-xmu(1,2,k)
504         xmu(3,2,k)=-xmu(2,3,k)
505      enddo
506      do i=1,3
507         do j=1,3
508            xmu(i,j,nope)=0.d0
509         enddo
510      enddo
511!
512!     correction due to change of xi and eta
513!
514      do k=1,nope
515         do j=1,3
516            do i=1,3
517!
518!     modified by Stefan Sicklinger
519!
520               xmu(i,j,k)=xmu(i,j,k)+(qxxy(i)+qxyx(i))*dxi(j,k)
521     &                                +(qxyy(i)+qyxy(i))*det(j,k)
522            enddo
523         enddo
524      enddo
525!
526!     derivatives of the size of the jacobian vector w.r.t. the
527!     displacement vectors (d ||m||/d u_k)
528!
529      do k=1,nope
530         do i=1,3
531            dxmu(i,k)=xn(1)*xmu(1,i,k)+xn(2)*xmu(2,i,k)+
532     &           xn(3)*xmu(3,i,k)
533         enddo
534!
535!        auxiliary variable: (d clear d u_k)*||m||
536!
537         do i=1,3
538            dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k)+
539     &               al(3)*xmu(3,i,k)-clear*dxmu(i,k)+
540     &               xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k)+xm(3)*dal(3,i,k)
541         enddo
542!
543      enddo
544!
545      c1=1.d0/dm
546      c2=c1*c1
547      c3=elas(2)*c2
548      c4=elas(1)*c1
549!
550!     derivatives of the forces w.r.t. the displacement vectors
551!
552      do k=1,nope
553         do j=1,3
554            do i=1,3
555               fpu(i,j,k)=-c3*xm(i)*dval(j,k)
556     &                    +c4*(xn(i)*dxmu(j,k)-xmu(i,j,k))
557            enddo
558         enddo
559      enddo
560!
561!     Coulomb friction for static calculations
562!
563      if(ncmat_.ge.7) then
564         um=elcon(6,1,imat)
565         if(um.gt.0.d0) then
566!
567!     stiffness of shear stress versus slip curve
568!
569            xk=elcon(7,1,imat)*springarea(1)
570!
571!     calculating the relative displacement between the slave node
572!     and its projection on the master surface
573!
574            do i=1,3
575               alnew(i)=voldl(i,nope)
576               do j=1,nterms
577                  alnew(i)=alnew(i)-ratio(j)*voldl(i,j)
578               enddo
579            enddo
580!
581!     calculating the difference in relative displacement since
582!     the start of the increment = vec_s
583!
584            do i=1,3
585               al(i)=alnew(i)-xstateini(3+i,1,ne0+konl(nope+1))
586            enddo
587!
588!     s=||vec_s||
589!
590            val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
591!
592!     update the relative tangential displacement vec_t
593!
594            do i=1,3
595               t(i)=xstateini(6+i,1,ne0+konl(nope+1))+al(i)-val*xn(i)
596            enddo
597!
598!     store the actual relative displacement and
599!     the actual relative tangential displacement
600!
601            do i=1,3
602               xstate(3+i,1,ne0+konl(nope+1))=alnew(i)
603               xstate(6+i,1,ne0+konl(nope+1))=t(i)
604            enddo
605!
606!     d vec_s/ d vec_u_k
607!     notice: xi & et are const.
608!
609            do k=1,nope
610               do i=1,3
611                  do j=1,3
612                     dal(i,j,k)=0.d0
613                  enddo
614               enddo
615            enddo
616!
617            do i=1,nterms
618               do j=1,3
619                  dal(j,j,i)=-shp2(4,i)
620               enddo
621            enddo
622!
623            do j=1,3
624               dal(j,j,nope)=1.d0
625            enddo
626!
627!     d s/ dvec_u_k.||m||
628!
629            do k=1,nope
630               do i=1,3
631                  dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k)
632     &                 +al(3)*xmu(3,i,k)-val*dxmu(i,k)
633     &                 +xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k)
634     &                 +xm(3)*dal(3,i,k)
635               enddo
636            enddo
637!
638!     d vec_t/d vec_u_k
639!
640            do k=1,nope
641               do j=1,3
642                  do i=1,3
643                     tu(i,j,k)=dal(i,j,k)
644     &                    -c1*(xn(i)*(dval(j,k)-val*dxmu(j,k))
645     &                    +val*xmu(i,j,k))
646                  enddo
647               enddo
648            enddo
649!
650!     size of normal force
651!
652            dfnl=dsqrt(fnl(1)**2+fnl(2)**2+fnl(3)**2)
653!
654!     maximum size of shear force
655!
656            dfshear=um*dfnl
657!
658!     plastic and elastic slip
659!
660            do i=1,3
661               tp(i)=xstateini(i,1,ne0+konl(nope+1))
662               te(i)=t(i)-tp(i)
663            enddo
664!
665!     the force due to normal contact is in -xn
666!     direction (internal force) -> minus signs
667!
668            do k=1,nope
669               do i=1,3
670                  dfn(i,k)=-xn(1)*fpu(1,i,k)-xn(2)*fpu(2,i,k)-
671     &                 xn(3)*fpu(3,i,k)
672               enddo
673            enddo
674!
675            dte=dsqrt(te(1)*te(1)+te(2)*te(2)+te(3)*te(3))
676!
677!     trial force
678!
679            do i=1,3
680               ftrial(i)=xk*te(i)
681            enddo
682            dftrial=dsqrt(ftrial(1)**2+ftrial(2)**2+ftrial(3)**2)
683!
684!     check whether stick or slip
685!
686            if((dftrial.lt.dfshear) .or. (dftrial.le.0.d0)) then
687!
688!     stick force
689!
690               do i=1,3
691                  fnl(i)=fnl(i)+ftrial(i)
692                  xstate(i,1,ne0+konl(nope+1))=tp(i)
693               enddo
694!
695!     stick stiffness
696!
697               do k=1,nope
698                  do j=1,3
699                     do i=1,3
700                        fpu(i,j,k)=fpu(i,j,k)+xk*tu(i,j,k)
701                     enddo
702                  enddo
703               enddo
704            else
705!
706!     slip force
707!
708               dg=(dftrial-dfshear)/xk
709               do i=1,3
710                  ftrial(i)=te(i)/dte
711                  fnl(i)=fnl(i)+dfshear*ftrial(i)
712                  xstate(i,1,ne0+konl(nope+1))=tp(i)+dg*ftrial(i)
713               enddo
714!
715!     slip stiffness
716!
717               c1=xk*dfshear/dftrial
718               do i=1,3
719                  do j=1,3
720                     dftdt(i,j)=-c1*ftrial(i)*ftrial(j)
721                  enddo
722                  dftdt(i,i)=dftdt(i,i)+c1
723               enddo
724!
725               do k=1,nope
726                  do j=1,3
727                     do i=1,3
728                        do l=1,3
729                           fpu(i,j,k)=fpu(i,j,k)+dftdt(i,l)*tu(l,j,k)
730                        enddo
731                        if((nmethod.ne.4).or.(iperturb(1).gt.1)) then
732                           fpu(i,j,k)=fpu(i,j,k)+um*ftrial(i)*dfn(j,k)
733                        endif
734                     enddo
735                  enddo
736               enddo
737            endif
738         endif
739      endif
740!
741!     determining the stiffness matrix contributions
742!
743!     complete field shp2
744!
745      shp2(1,nope)=0.d0
746      shp2(2,nope)=0.d0
747      shp2(4,nope)=-1.d0
748!
749      do k=1,nope
750         do l=1,nope
751            do i=1,3
752               i1=i+(k-1)*3
753               do j=1,3
754                  s(i1,j+(l-1)*3)=-shp2(4,k)*fpu(i,j,l)
755     &                            -shp2(1,k)*fnl(i)*dxi(j,l)
756     &                            -shp2(2,k)*fnl(i)*det(j,l)
757               enddo
758            enddo
759         enddo
760      enddo
761!
762!     symmetrizing the matrix
763!     this is done in the absence of friction or for modal dynamic
764!     calculations
765!
766      if((nasym.eq.0).or.((nmethod.eq.4).and.(iperturb(1).le.1))) then
767         do j=1,3*nope
768            do i=1,j-1
769               s(i,j)=(s(i,j)+s(j,i))/2.d0
770            enddo
771         enddo
772      endif
773!
774      return
775      end
776
777