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 springdamp_n2f(xl,elas,voldl,s,imat,elcon,
20     &  ncmat_,ntmat_,nope,iperturb,springarea,nmethod,mi,
21     &  reltime,nasym)
22!
23!     calculates the contact damping matrix (node-to-face penalty)
24!     (User's manual -> theory -> boundary conditions ->
25!      node-to-face penalty contact)
26!
27      implicit none
28!
29      integer i,j,imat,ncmat_,ntmat_,k,l,nope,nterms,iflag,i1,
30     &  iperturb(*),nmethod,mi(*),nasym
31!
32      real*8 xl(3,10),elas(21),ratio(9),pproj(3),val,shp2(7,9),
33     &  al(3),s(60,60),voldl(0:mi(2),10),pl(3,10),xn(3),dm,
34     &  c1,c2,c3,c4,elcon(0:ncmat_,ntmat_,*),xm(3),xmu(3,3,10),
35     &  dxmu(3,10),dval(3,10),fpu(3,3,10),xi,et,xs2(3,7),xk,
36     &  a11,a12,a22,b1(3,10),b2(3,10),dal(3,3,10),qxxy(3),fnl(3),
37     &  qxyy(3),dxi(3,10),det(3,10),determinant,c11,c12,c22,
38     &  qxyx(3),qyxy(3),springarea(2),dist,tu(3,3,10),
39     &  clear,reltime,alnew(3)
40!
41      data iflag /4/
42!
43!     actual positions of the nodes belonging to the contact spring
44!     (otherwise no contact force)
45!
46      do i=1,nope
47         do j=1,3
48            pl(j,i)=xl(j,i)+voldl(j,i)
49         enddo
50      enddo
51!
52!     contact springs
53!
54      nterms=nope-1
55!
56!     vector al connects the actual position of the slave node
57!     with its projection on the master face = vec_r (User's
58!     manual -> theory -> boundary conditions -> node-to-face
59!     penalty contact)
60!
61      do i=1,3
62         pproj(i)=pl(i,nope)
63      enddo
64      call attach_2d(pl,pproj,nterms,ratio,dist,xi,et)
65      do i=1,3
66         al(i)=pl(i,nope)-pproj(i)
67      enddo
68!
69!     determining the jacobian vector on the surface
70!
71      if(nterms.eq.8) then
72         call shape8q(xi,et,pl,xm,xs2,shp2,iflag)
73      elseif(nterms.eq.4) then
74         call shape4q(xi,et,pl,xm,xs2,shp2,iflag)
75      elseif(nterms.eq.6) then
76         call shape6tri(xi,et,pl,xm,xs2,shp2,iflag)
77      else
78         call shape3tri(xi,et,pl,xm,xs2,shp2,iflag)
79      endif
80!
81!     d xi / d vec_u_j
82!     d eta / d vec_u_j
83!
84!     dxi and det are determined from the orthogonality
85!     condition
86!
87      a11=-(xs2(1,1)*xs2(1,1)+xs2(2,1)*xs2(2,1)+xs2(3,1)*xs2(3,1))
88     &    +al(1)*xs2(1,5)+al(2)*xs2(2,5)+al(3)*xs2(3,5)
89      a12=-(xs2(1,1)*xs2(1,2)+xs2(2,1)*xs2(2,2)+xs2(3,1)*xs2(3,2))
90     &    +al(1)*xs2(1,6)+al(2)*xs2(2,6)+al(3)*xs2(3,6)
91      a22=-(xs2(1,2)*xs2(1,2)+xs2(2,2)*xs2(2,2)+xs2(3,2)*xs2(3,2))
92     &    +al(1)*xs2(1,7)+al(2)*xs2(2,7)+al(3)*xs2(3,7)
93!
94      do i=1,3
95         do j=1,nterms
96            b1(i,j)=shp2(4,j)*xs2(i,1)-shp2(1,j)*al(i)
97            b2(i,j)=shp2(4,j)*xs2(i,2)-shp2(2,j)*al(i)
98         enddo
99         b1(i,nope)=-xs2(i,1)
100         b2(i,nope)=-xs2(i,2)
101      enddo
102!
103      determinant=a11*a22-a12*a12
104      c11=a22/determinant
105      c12=-a12/determinant
106      c22=a11/determinant
107!
108      do i=1,3
109         do j=1,nope
110            dxi(i,j)=c11*b1(i,j)+c12*b2(i,j)
111            det(i,j)=c12*b1(i,j)+c22*b2(i,j)
112         enddo
113      enddo
114!
115!     d vec_r / d vec_u_k
116!
117      do i=1,nope
118         do j=1,3
119            do k=1,3
120               dal(j,k,i)=-xs2(j,1)*dxi(k,i)-xs2(j,2)*det(k,i)
121            enddo
122         enddo
123      enddo
124      do i=1,nterms
125         do j=1,3
126            dal(j,j,i)=dal(j,j,i)-shp2(4,i)
127         enddo
128      enddo
129      do j=1,3
130         dal(j,j,nope)=dal(j,j,nope)+1.d0
131      enddo
132!
133!     d2q/dxx x dq/dy
134!
135      qxxy(1)=xs2(2,5)*xs2(3,2)-xs2(3,5)*xs2(2,2)
136      qxxy(2)=xs2(3,5)*xs2(1,2)-xs2(1,5)*xs2(3,2)
137      qxxy(3)=xs2(1,5)*xs2(2,2)-xs2(2,5)*xs2(1,2)
138!
139!     dq/dx x d2q/dyy
140!
141      qxyy(1)=xs2(2,1)*xs2(3,7)-xs2(3,1)*xs2(2,7)
142      qxyy(2)=xs2(3,1)*xs2(1,7)-xs2(1,1)*xs2(3,7)
143      qxyy(3)=xs2(1,1)*xs2(2,7)-xs2(2,1)*xs2(1,7)
144!
145!     Modified by Stefan Sicklinger
146!
147!     dq/dx x d2q/dxy
148!
149      qxyx(1)=xs2(2,1)*xs2(3,6)-xs2(3,1)*xs2(2,6)
150      qxyx(2)=xs2(3,1)*xs2(1,6)-xs2(1,1)*xs2(3,6)
151      qxyx(3)=xs2(1,1)*xs2(2,6)-xs2(2,1)*xs2(1,6)
152!
153!     d2q/dxy x dq/dy
154!
155      qyxy(1)=xs2(2,6)*xs2(3,2)-xs2(3,6)*xs2(2,2)
156      qyxy(2)=xs2(3,6)*xs2(1,2)-xs2(1,6)*xs2(3,2)
157      qyxy(3)=xs2(1,6)*xs2(2,2)-xs2(2,6)*xs2(1,2)
158!
159!
160!     End modifications
161!
162!     normal on the surface
163!
164      dm=dsqrt(xm(1)*xm(1)+xm(2)*xm(2)+xm(3)*xm(3))
165      do i=1,3
166         xn(i)=xm(i)/dm
167      enddo
168!
169!     distance from surface along normal (= clearance)
170!
171      clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
172      if(nmethod.eq.1) then
173         clear=clear-springarea(2)*(1.d0-reltime)
174      endif
175!
176!     representative area: usually the slave surface stored in
177!     springarea; however, if no area was assigned because the
178!     node does not belong to any element, the master surface
179!     is used
180!
181      if(springarea(1).le.0.d0) then
182         if(nterms.eq.3) then
183            springarea(1)=dm/2.d0
184         else
185            springarea(1)=dm*4.d0
186         endif
187      endif
188!
189      elas(2)=-springarea(1)*elcon(5,1,imat)
190      elas(1)=elas(2)*clear
191!
192!     contact force
193!
194      do i=1,3
195         fnl(i)=-elas(1)*xn(i)
196      enddo
197!
198!     derivatives of the jacobian vector w.r.t. the displacement
199!     vectors (d m / d u_k)
200!
201      do k=1,nterms
202         xmu(1,1,k)=0.d0
203         xmu(2,2,k)=0.d0
204         xmu(3,3,k)=0.d0
205         xmu(1,2,k)=shp2(1,k)*xs2(3,2)-shp2(2,k)*xs2(3,1)
206         xmu(2,3,k)=shp2(1,k)*xs2(1,2)-shp2(2,k)*xs2(1,1)
207         xmu(3,1,k)=shp2(1,k)*xs2(2,2)-shp2(2,k)*xs2(2,1)
208         xmu(1,3,k)=-xmu(3,1,k)
209         xmu(2,1,k)=-xmu(1,2,k)
210         xmu(3,2,k)=-xmu(2,3,k)
211      enddo
212      do i=1,3
213         do j=1,3
214            xmu(i,j,nope)=0.d0
215         enddo
216      enddo
217!
218!     correction due to change of xi and eta
219!
220      do k=1,nope
221         do j=1,3
222            do i=1,3
223!
224!     modified by Stefan Sicklinger
225!
226               xmu(i,j,k)=xmu(i,j,k)+(qxxy(i)+qxyx(i))*dxi(j,k)
227     &                                +(qxyy(i)+qyxy(i))*det(j,k)
228            enddo
229         enddo
230      enddo
231!
232!     derivatives of the size of the jacobian vector w.r.t. the
233!     displacement vectors (d ||m||/d u_k)
234!
235      do k=1,nope
236         do i=1,3
237            dxmu(i,k)=xn(1)*xmu(1,i,k)+xn(2)*xmu(2,i,k)+
238     &           xn(3)*xmu(3,i,k)
239         enddo
240!
241!        auxiliary variable: (d clear d u_k)*||m||
242!
243         do i=1,3
244            dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k)+
245     &               al(3)*xmu(3,i,k)-clear*dxmu(i,k)+
246     &               xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k)+xm(3)*dal(3,i,k)
247         enddo
248!
249      enddo
250!
251      c1=1.d0/dm
252      c2=c1*c1
253      c3=elas(2)*c2
254      c4=elas(1)*c1
255!
256!     derivatives of the forces w.r.t. the displacement vectors
257!
258      do k=1,nope
259         do j=1,3
260            do i=1,3
261               fpu(i,j,k)=-c3*xm(i)*dval(j,k)
262     &                    +c4*(xn(i)*dxmu(j,k)-xmu(i,j,k))
263            enddo
264         enddo
265      enddo
266!
267!     tangential damping
268!
269      if(elcon(8,1,imat).gt.0.d0) then
270!
271!     stiffness of shear stress versus slip curve
272!
273         xk=elcon(8,1,imat)*elcon(5,1,imat)*springarea(1)
274!
275!     calculating the relative displacement between the slave node
276!     and its projection on the master surface
277!
278         do i=1,3
279            alnew(i)=voldl(i,nope)
280            do j=1,nterms
281               alnew(i)=alnew(i)-ratio(j)*voldl(i,j)
282            enddo
283         enddo
284!
285!     calculating the difference in relative displacement since
286!     the start of the increment = vec_s
287!
288         do i=1,3
289            al(i)=alnew(i)
290         enddo
291!
292!     s=||vec_s||
293!
294         val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
295!
296!     d vec_s/ d vec_u_k
297!     notice: xi & et are const.
298!
299         do k=1,nope
300            do i=1,3
301               do j=1,3
302                  dal(i,j,k)=0.d0
303               enddo
304            enddo
305         enddo
306!
307         do i=1,nterms
308            do j=1,3
309               dal(j,j,i)=-shp2(4,i)
310            enddo
311         enddo
312!
313         do j=1,3
314            dal(j,j,nope)=1.d0
315         enddo
316!
317!     d s/ dvec_u_k.||m||
318!
319         do k=1,nope
320            do i=1,3
321               dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k)
322     &              +al(3)*xmu(3,i,k)-val*dxmu(i,k)
323     &              +xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k)
324     &              +xm(3)*dal(3,i,k)
325            enddo
326         enddo
327!
328!     d vec_t/d vec_u_k
329!
330         do k=1,nope
331            do j=1,3
332               do i=1,3
333                  tu(i,j,k)=dal(i,j,k)
334     &                 -c1*(xn(i)*(dval(j,k)-val*dxmu(j,k))
335     &                 +val*xmu(i,j,k))
336               enddo
337            enddo
338         enddo
339!
340!     damping matrix
341!
342         do k=1,nope
343            do j=1,3
344               do i=1,3
345                  fpu(i,j,k)=fpu(i,j,k)+xk*tu(i,j,k)
346               enddo
347            enddo
348         enddo
349      endif
350!
351!     determining the damping matrix contributions
352!
353!     complete field shp2
354!
355      shp2(1,nope)=0.d0
356      shp2(2,nope)=0.d0
357      shp2(4,nope)=-1.d0
358!
359      do k=1,nope
360         do l=1,nope
361            do i=1,3
362               i1=i+(k-1)*3
363               do j=1,3
364                  s(i1,j+(l-1)*3)=-shp2(4,k)*fpu(i,j,l)
365     &                            -shp2(1,k)*fnl(i)*dxi(j,l)
366     &                            -shp2(2,k)*fnl(i)*det(j,l)
367               enddo
368            enddo
369         enddo
370      enddo
371!
372!     symmetrizing the matrix
373!     this is done in the absence of friction or for modal dynamic
374!     calculations
375!
376      if((nasym.eq.0).or.((nmethod.eq.4).and.(iperturb(1).le.1))) then
377         do j=1,3*nope
378            do i=1,j-1
379               s(i,j)=(s(i,j)+s(j,i))/2.d0
380            enddo
381         enddo
382      endif
383!
384      return
385      end
386
387