FrontISTR  5.7.1
Large-scale structural analysis program with finit element method
fstr_AddContactStiff.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
12 
14 
15  use m_fstr
16  use elementinfo
17  use m_contact_lib
20  use m_fstr_residual
21 
22  implicit none
23 
24  public :: fstr_addcontactstiffness
26  public :: update_ndforce_contact
27  public :: fstr_ass_load_contact
28 
29  private :: getcontactstiffness
30  private :: getcontactnodalforce
31  private :: gettrialfricforceandcheckfricstate
32 
33 contains
34 
37  subroutine fstr_addcontactstiffness(cstep,iter,hecMAT,hecLagMAT,fstrSOLID)
38 
39  integer(kind=kint) :: cstep
40  integer(kind=kint) :: iter
41  type(hecmwst_matrix) :: hecmat
42  type(fstr_solid) :: fstrsolid
43  type(hecmwst_matrix_lagrange) :: heclagmat
44  integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(21)
45  integer(kind=kint) :: i, j, k, nlag, id_lagrange, grpid, algtype
46  real(kind=kreal) :: lagrange
47  real(kind=kreal) :: stiffness(21*3 + 1, 21*3 + 1)
48 
49  heclagmat%AL_lagrange = 0.0d0
50  heclagmat%AU_lagrange = 0.0d0
51 
52  id_lagrange = 0
53 
54  do i = 1, fstrsolid%n_contacts
55 
56  grpid = fstrsolid%contacts(i)%group
57  if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
58 
59  algtype = fstrsolid%contacts(i)%algtype
60  nlag = fstr_get_num_lagrange_pernode(algtype)
61 
62  do j = 1, size(fstrsolid%contacts(i)%slave)
63 
64  if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
65 
66  ctsurf = fstrsolid%contacts(i)%states(j)%surface
67  etype = fstrsolid%contacts(i)%master(ctsurf)%etype
68  nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
69  ndlocal(1) = fstrsolid%contacts(i)%slave(j)
70  ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
71 
72  do k=1,nlag
73  id_lagrange = id_lagrange + 1
74  lagrange = heclagmat%Lagrange(id_lagrange)
75 
76  if( algtype == contactsslid .or. algtype == contactfslid ) then
77  ! Obtain contact stiffness matrix of contact pair
78  call getcontactstiffness(iter,etype,nnode,fstrsolid%contacts(i)%states(j), &
79  fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,stiffness)
80  else if( algtype == contacttied ) then
81  call gettiedstiffness(etype,nnode,k,fstrsolid%contacts(i)%states(j),lagrange,stiffness)
82  endif
83 
84  ! Assemble contact stiffness matrix of contact pair into global stiffness matrix
85  call hecmw_mat_ass_contactlag(nnode,ndlocal,id_lagrange,fstrsolid%contacts(i)%fcoeff,stiffness,hecmat,heclagmat)
86  enddo
87  enddo
88  enddo
89 
90  do i = 1, fstrsolid%n_embeds
91 
92  grpid = fstrsolid%embeds(i)%group
93  if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
94 
95  do j = 1, size(fstrsolid%embeds(i)%slave)
96 
97  if( fstrsolid%embeds(i)%states(j)%state == contactfree ) cycle
98 
99  ctsurf = fstrsolid%embeds(i)%states(j)%surface
100  etype = fstrsolid%embeds(i)%master(ctsurf)%etype
101  nnode = size(fstrsolid%embeds(i)%master(ctsurf)%nodes)
102  ndlocal(1) = fstrsolid%embeds(i)%slave(j)
103  ndlocal(2:nnode+1) = fstrsolid%embeds(i)%master(ctsurf)%nodes(1:nnode)
104 
105  do k=1,3
106  id_lagrange = id_lagrange + 1
107  lagrange = heclagmat%Lagrange(id_lagrange)
108 
109  call gettiedstiffness(etype,nnode,k,fstrsolid%embeds(i)%states(j),lagrange,stiffness)
110 
111  ! Assemble contact stiffness matrix of contact pair into global stiffness matrix
112  call hecmw_mat_ass_contactlag(nnode,ndlocal,id_lagrange,0.d0,stiffness,hecmat,heclagmat)
113  enddo
114  enddo
115  enddo
116 
117  end subroutine fstr_addcontactstiffness
118 
120  subroutine gettiedstiffness(etype,nnode,idof,ctState,lagrange,stiffness)
121  integer(kind=kint) :: etype
122  integer(kind=kint) :: nnode
123  integer(kind=kint) :: idof
124  type(tcontactstate) :: ctState
125  real(kind=kreal) :: lagrange
126  real(kind=kreal) :: stiffness(:,:)
127 
128  integer(kind=kint) :: i, j, k, l
129  real(kind=kreal) :: shapefunc(nnode)
130  real(kind=kreal) :: ntm((nnode+1)*3)
131 
132  stiffness = 0.0d0
133 
134  call getshapefunc( etype, ctstate%lpos(:), shapefunc )
135 
136  ntm = 0.d0
137  ntm(idof) = 1.d0
138  do i = 1, nnode
139  ntm(i*3+idof) = -shapefunc(i)
140  enddo
141 
142  i = (nnode+1)*3 + 1
143  do j = 1, (nnode+1)*3
144  stiffness(i,j) = ntm(j); stiffness(j,i) = ntm(j)
145  enddo
146 
147  end subroutine gettiedstiffness
148 
150  subroutine getcontactstiffness(iter,etype,nnode,ctState,tPenalty,fcoeff,lagrange,stiffness)
151 
152  type(tcontactstate) :: ctstate
153  integer(kind=kint) :: iter
154  integer(kind=kint) :: etype, nnode
155  integer(kind=kint) :: i, j, k, l
156  real(kind=kreal) :: normal(3), shapefunc(nnode)
157  real(kind=kreal) :: ntm((nnode + 1)*3)
158  real(kind=kreal) :: fcoeff, tpenalty
159  real(kind=kreal) :: lagrange
160  real(kind=kreal) :: tf_trial(3), length_tft
161  real(kind=kreal) :: tangent(3), ttm((nnode + 1)*3)
162  real(kind=kreal) :: stiffness(:,:)
163 
164  stiffness = 0.0d0
165 
166  call getshapefunc( etype, ctstate%lpos(:), shapefunc )
167 
168  normal(1:3) = ctstate%direction(1:3)
169 
170  ntm(1:3) = normal(1:3)
171  do i = 1, nnode
172  ntm(i*3+1:i*3+3) = -shapefunc(i)*normal(1:3)
173  enddo
174 
175  i = (nnode+1)*3 + 1
176  do j = 1, (nnode+1)*3
177  stiffness(i,j) = ntm(j); stiffness(j,i) = ntm(j)
178  enddo
179 
180 
181  if( fcoeff /= 0.0d0 ) then
182  if( lagrange>0.0d0 .or. iter==1 ) then
183 
184  do i = 1, nnode+1
185  do j = 1, i
186  do k = 1, 3
187  do l = 1, 3
188  stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) - tpenalty*ntm((i-1)*3+k)*ntm((j-1)*3+l)
189  if( k==l ) then
190  if(i==1 .and. j==1)then
191  stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) + tpenalty
192  elseif(i>1 .and. j==1)then
193  stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) - tpenalty*shapefunc(i-1)
194  elseif(i>1 .and. j>1)then
195  stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) + tpenalty*shapefunc(i-1)*shapefunc(j-1)
196  endif
197  endif
198  if(i==j) cycle
199  stiffness((j-1)*3+l,(i-1)*3+k) = stiffness((i-1)*3+k,(j-1)*3+l)
200  enddo
201  enddo
202  enddo
203  enddo
204 
205  if( ctstate%state == contactslip ) then
206 
207  tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
208  length_tft = dsqrt(dot_product(tf_trial,tf_trial))
209  tangent(1:3) = tf_trial(1:3)/length_tft
210 
211  ttm(1:3) = -tangent(1:3)
212  do i = 1, nnode
213  ttm(i*3+1:i*3+3) = shapefunc(i)*tangent(1:3)
214  enddo
215 
216  do i = 1, nnode+1
217  do j = 1, nnode+1
218  do k = 1, 3
219  do l = 1, 3
220  stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) &
221  + tpenalty*(-ttm((i-1)*3+k)*ttm((j-1)*3+l) &
222  +ttm((i-1)*3+k)*ntm((j-1)*3+l)*dot_product(tangent,normal))
223  enddo
224  enddo
225  enddo
226  enddo
227  stiffness(1:(nnode+1)*3,1:(nnode+1)*3) = (fcoeff*lagrange/length_tft)*stiffness(1:(nnode+1)*3,1:(nnode+1)*3)
228 
229  !
230  ! do i = 1, (nnode + 1)*3
231  ! do j = 1, i
232  ! do k = 1, 3
233  ! do l = 1, k
234  ! stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) - mut*tTm((i-1)*3+k)*tTm((j-1)*3+l))
235  ! if( k/=l )stiffness((j-1)*3+l,(i-1)*3+k) = stiffness((i-1)*3+k,(j-1)*3+l)
236  ! enddo !- l
237  ! enddo !- k
238  ! enddo !- j
239  ! enddo !- i
240  ! stiffness(1:(nnode+1)*3,1:(nnode+1)*3) = (fcoeff*dabs(lagrange)/length_tft)*stiffness(1:(nnode+1)*3,1:(nnode+1)*3)
241 
242  ! j = (nnode+1)*3 + 1
243  ! do i = 1, (nnode+1)*3
244  ! stiffness(i,j) = stiffness(i,j) - fcoeff*tTm(i)
245  ! enddo
246  stiffness(1:(nnode+1)*3,(nnode+1)*3+1) = stiffness(1:(nnode+1)*3,(nnode+1)*3+1) - fcoeff*ttm(1:(nnode+1)*3)
247 
248  endif
249  endif
250  endif
251 
252  end subroutine getcontactstiffness
253 
256  subroutine fstr_update_ndforce_contact(cstep,hecMESH,hecMAT,hecLagMAT,fstrSOLID,conMAT)
257 
258  type(hecmwst_local_mesh) :: hecmesh
259  type(hecmwst_matrix) :: hecmat
260  type(fstr_solid) :: fstrsolid
261  type(hecmwst_matrix_lagrange) :: heclagmat
262  type(hecmwst_matrix) :: conmat
263  integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(21)
264  integer(kind=kint) :: i, j, k, nlag, id_lagrange, algtype
265  real(kind=kreal) :: ndcoord(21*3)
266  real(kind=kreal) :: ndu(21*3), nddu(21*3)
267  real(kind=kreal) :: lagrange
268  real(kind=kreal) :: ctnforce(21*3+1)
269  real(kind=kreal) :: cttforce(21*3+1)
270 
271  integer(kind=kint) :: cstep
272  integer(kind=kint) :: grpid
273  logical :: if_flag
274  real(kind=kreal) :: ctime, etime
275  integer(kind=kint) :: if_type
276 
277  id_lagrange = 0
278  if( associated(fstrsolid%CONT_NFORCE) ) fstrsolid%CONT_NFORCE(:) = 0.d0
279  if( associated(fstrsolid%CONT_FRIC) ) fstrsolid%CONT_FRIC(:) = 0.d0
280  if( associated(fstrsolid%EMBED_NFORCE) ) fstrsolid%EMBED_NFORCE(:) = 0.d0
281 
282  do i = 1, fstrsolid%n_contacts
283 
284  grpid = fstrsolid%contacts(i)%group
285  if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
286 
287  algtype = fstrsolid%contacts(i)%algtype
288  nlag = fstr_get_num_lagrange_pernode(algtype)
289  if_flag = (fstrsolid%contacts(i)%if_type /= 0)
290  if(if_flag)then
291  ctime = fstrsolid%contacts(i)%ctime
292  etime = fstrsolid%contacts(i)%if_etime
293  if_type = fstrsolid%contacts(i)%if_type
294  end if
295 
296  do j = 1, size(fstrsolid%contacts(i)%slave)
297 
298  if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
299  if(if_flag) call set_shrink_factor(ctime, fstrsolid%contacts(i)%states(j), etime, if_type)
300 
301  ctsurf = fstrsolid%contacts(i)%states(j)%surface
302  etype = fstrsolid%contacts(i)%master(ctsurf)%etype
303  nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
304  ndlocal(1) = fstrsolid%contacts(i)%slave(j)
305  ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
306  do k = 1, nnode+1
307  nddu((k-1)*3+1:(k-1)*3+3) = fstrsolid%dunode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
308  ndu((k-1)*3+1:(k-1)*3+3) = fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) + nddu((k-1)*3+1:(k-1)*3+3)
309  ndcoord((k-1)*3+1:(k-1)*3+3) = hecmesh%node((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) + ndu((k-1)*3+1:(k-1)*3+3)
310  enddo
311 
312  do k=1,nlag
313  id_lagrange = id_lagrange + 1
314  lagrange = heclagmat%Lagrange(id_lagrange)
315 
316  fstrsolid%contacts(i)%states(j)%multiplier(k) = heclagmat%Lagrange(id_lagrange)
317 
318  if( algtype == contactsslid .or. algtype == contactfslid ) then
319  ! Obtain contact nodal force vector of contact pair
320  if(if_flag) call get_shrink_dummy_surf(fstrsolid%contacts(i)%states(j),ndcoord, nnode)
321 
322  call getcontactnodalforce(etype,nnode,ndcoord,nddu,fstrsolid%contacts(i)%states(j), &
323  fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,ctnforce,cttforce,.true.)
324  ! Update non-eqilibrited force vector
325  call update_ndforce_contact(nnode,ndlocal,id_lagrange,lagrange,ctnforce,cttforce, &
326  & conmat,fstrsolid%CONT_NFORCE,fstrsolid%CONT_FRIC)
327  else if( algtype == contacttied ) then
328  call gettiednodalforce(etype,nnode,k,ndu,fstrsolid%contacts(i)%states(j),lagrange,ctnforce,cttforce)
329  ! Update non-eqilibrited force vector
330  call update_ndforce_contact(nnode,ndlocal,id_lagrange,1.d0,ctnforce,cttforce, &
331  & conmat,fstrsolid%CONT_NFORCE,fstrsolid%CONT_FRIC)
332  endif
333 
334  enddo
335  enddo
336 
337  enddo
338 
339  do i = 1, fstrsolid%n_embeds
340 
341  grpid = fstrsolid%embeds(i)%group
342  if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
343 
344  do j = 1, size(fstrsolid%embeds(i)%slave)
345 
346  if( fstrsolid%embeds(i)%states(j)%state == contactfree ) cycle
347 
348  ctsurf = fstrsolid%embeds(i)%states(j)%surface
349  etype = fstrsolid%embeds(i)%master(ctsurf)%etype
350  nnode = size(fstrsolid%embeds(i)%master(ctsurf)%nodes)
351  ndlocal(1) = fstrsolid%embeds(i)%slave(j)
352  ndlocal(2:nnode+1) = fstrsolid%embeds(i)%master(ctsurf)%nodes(1:nnode)
353  do k = 1, nnode+1
354  nddu((k-1)*3+1:(k-1)*3+3) = fstrsolid%dunode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
355  ndu((k-1)*3+1:(k-1)*3+3) = fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) + nddu((k-1)*3+1:(k-1)*3+3)
356  ndcoord((k-1)*3+1:(k-1)*3+3) = hecmesh%node((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) + ndu((k-1)*3+1:(k-1)*3+3)
357  enddo
358 
359  do k=1,3
360  id_lagrange = id_lagrange + 1
361  lagrange = heclagmat%Lagrange(id_lagrange)
362 
363  fstrsolid%embeds(i)%states(j)%multiplier(k) = heclagmat%Lagrange(id_lagrange)
364 
365  call gettiednodalforce(etype,nnode,k,ndu,fstrsolid%embeds(i)%states(j),lagrange,ctnforce,cttforce)
366  ! Update non-eqilibrited force vector
367  call update_ndforce_contact(nnode,ndlocal,id_lagrange,-1.d0,ctnforce,cttforce, &
368  & conmat,fstrsolid%EMBED_NFORCE,fstrsolid%EMBED_NFORCE)
369 
370  enddo
371  enddo
372  enddo
373 
374  ! Consider SPC condition
375  call fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, hecmat%B)
376  call fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, conmat%B)
377 
378  end subroutine fstr_update_ndforce_contact
379 
381  subroutine gettiednodalforce(etype,nnode,idof,ndu,ctState,lagrange,ctNForce,ctTForce)
382  integer(kind=kint) :: etype
383  integer(kind=kint) :: nnode
384  integer(kind=kint) :: idof
385  real(kind=kreal) :: ndu((nnode + 1)*3)
386  type(tcontactstate) :: ctstate
387  real(kind=kreal) :: lagrange
388  real(kind=kreal) :: ctnforce((nnode+1)*3+1)
389  real(kind=kreal) :: cttforce((nnode+1)*3+1)
390 
391  integer(kind=kint) :: i, j
392  real(kind=kreal) :: normal(3), shapefunc(nnode)
393  real(kind=kreal) :: ntm((nnode + 1)*3)
394  real(kind=kreal) :: ttm((nnode + 1)*3)
395 
396  ctnforce = 0.0d0
397  cttforce = 0.0d0
398 
399  call getshapefunc( etype, ctstate%lpos, shapefunc )
400 
401  normal(1:3) = ctstate%direction(1:3)
402 
403  ntm(1:3) = -normal(idof)*normal(1:3)
404  do i = 1, nnode
405  ntm(i*3+1:i*3+3) = -shapefunc(i)*ntm(1:3)
406  enddo
407  ttm = 0.d0
408  ttm(idof) = -1.0
409  ttm(1:3) = ttm(1:3)-ntm(1:3)
410  do i = 1, nnode
411  ttm(i*3+1:i*3+3) = -shapefunc(i)*ttm(1:3)
412  enddo
413 
414  do j = 1, (nnode+1)*3
415  ctnforce(j) = lagrange*ntm(j)
416  cttforce(j) = lagrange*ttm(j)
417  enddo
418  j = (nnode+1)*3 + 1
419  ctnforce(j) = dot_product(ntm,ndu)
420  cttforce(j) = dot_product(ttm,ndu)
421  end subroutine
422 
424  subroutine getcontactnodalforce(etype,nnode,ndCoord,ndDu,ctState,tPenalty,fcoeff,lagrange,ctNForce,ctTForce,cflag)
425 
426  type(tcontactstate) :: ctstate
427  integer(kind=kint) :: etype, nnode
428  integer(kind=kint) :: i, j
429  real(kind=kreal) :: fcoeff, tpenalty
430  real(kind=kreal) :: lagrange
431  real(kind=kreal) :: ndcoord((nnode + 1)*3), nddu((nnode + 1)*3)
432  real(kind=kreal) :: normal(3), shapefunc(nnode)
433  real(kind=kreal) :: ntm((nnode + 1)*3)
434  real(kind=kreal) :: tf_trial(3), length_tft, tangent(3), tf_final(3)
435  real(kind=kreal) :: ctnforce((nnode+1)*3+1)
436  real(kind=kreal) :: cttforce((nnode+1)*3+1)
437  logical :: cflag
438 
439  ctnforce = 0.0d0
440  cttforce = 0.0d0
441 
442  call getshapefunc( etype, ctstate%lpos, shapefunc )
443 
444  normal(1:3) = ctstate%direction(1:3)
445 
446  ntm(1:3) = -normal(1:3)
447  do i = 1, nnode
448  ntm(i*3+1:i*3+3) = shapefunc(i)*normal(1:3)
449  enddo
450 
451  do j = 1, (nnode+1)*3
452  ctnforce(j) = lagrange*ntm(j)
453  enddo
454  j = (nnode+1)*3 + 1
455  ctnforce(j) = dot_product(ntm,ndcoord)
456 
457  if(fcoeff /= 0.0d0 .and. lagrange > 0.0d0)then
458 
459  if( cflag ) then !calc tf_final and set it to tangentForce_final
460  call gettrialfricforceandcheckfricstate(nnode,tpenalty,fcoeff,lagrange,normal,shapefunc,ntm,nddu,ctstate)
461 
462  if( ctstate%state == contactstick ) then
463  tf_final(1:3) = ctstate%tangentForce_trial(1:3)
464  elseif( ctstate%state == contactslip ) then
465  tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
466  length_tft = dsqrt(dot_product(tf_trial,tf_trial))
467  tangent(1:3) = tf_trial(1:3)/length_tft
468  tf_final(1:3) = fcoeff*dabs(lagrange)*tangent(1:3)
469  endif
470  ctstate%tangentForce_final(1:3) = tf_final(1:3)
471  else ! just set tangentForce_final to tf_final (used for fstr_ass_load_contact)
472  tf_final(1:3) = ctstate%tangentForce_final(1:3)
473  endif
474 
475  cttforce(1:3) = -tf_final(1:3)
476  do j = 1, nnode
477  cttforce(j*3+1:j*3+3) = shapefunc(j)*tf_final(1:3)
478  enddo
479 
480  endif
481 
482  end subroutine getcontactnodalforce
483 
484 
486  subroutine gettrialfricforceandcheckfricstate(nnode,tPenalty,fcoeff,lagrange,normal,shapefunc,nTm,ndDu,ctstate)
487 
488  type(tcontactstate) :: ctState
489  integer(kind=kint) :: nnode
490  integer(kind=kint) :: i, j
491  real(kind=kreal) :: fcoeff, tpenalty
492  real(kind=kreal) :: lagrange
493  real(kind=kreal) :: nddu((nnode + 1)*3)
494  real(kind=kreal) :: normal(3), shapefunc(nnode)
495  real(kind=kreal) :: ntm((nnode + 1)*3)
496  real(kind=kreal) :: dotp
497  real(kind=kreal) :: relativedisp(3)
498  real(kind=kreal) :: tf_yield
499 
500  relativedisp = 0.0d0
501 
502  dotp = dot_product(ntm,nddu)
503  do i = 1, 3
504  relativedisp(i) = - nddu(i)
505  do j = 1, nnode
506  relativedisp(i) = relativedisp(i) + shapefunc(j)*nddu(j*3+i)
507  enddo
508  relativedisp(i) = relativedisp(i) - dotp*normal(i)
509  ctstate%reldisp(i) = -relativedisp(i)
510  ctstate%tangentForce_trial(i) = ctstate%tangentForce1(i) -tpenalty*relativedisp(i)
511  enddo
512 
513  tf_yield = fcoeff*dabs(lagrange)
514  if(ctstate%state == contactslip) tf_yield =0.99d0*tf_yield
515  if( dsqrt(dot_product(ctstate%tangentForce_trial,ctstate%tangentForce_trial)) <= tf_yield ) then
516  ctstate%state = contactstick
517  else
518  ctstate%state = contactslip
519  endif
520 
521  end subroutine gettrialfricforceandcheckfricstate
522 
523 
526  subroutine update_ndforce_contact(nnode,ndLocal,id_lagrange,lagrange,ctNForce,ctTForce,hecMAT,cont_nforce,cont_fric)
527 
528  type(hecmwst_matrix) :: hecmat
529  integer(kind=kint) :: nnode, ndlocal(nnode + 1)
531  integer(kind=kint) :: id_lagrange
532  real(kind=kreal) :: lagrange
533  integer(kind=kint) :: np, ndof
534  integer (kind=kint) :: i, inod, idx
535  real(kind=kreal) :: ctnforce((nnode+1)*3+1)
536  real(kind=kreal) :: cttforce((nnode+1)*3+1)
537  real(kind=kreal), pointer :: cont_nforce(:)
538  real(kind=kreal), pointer, optional :: cont_fric(:)
539 
540  np = hecmat%NP; ndof = hecmat%NDOF
541 
542  do i = 1, nnode + 1
543  inod = ndlocal(i)
544  idx = (inod-1)*3+1
545  hecmat%B(idx:idx+2) = hecmat%B(idx:idx+2) + ctnforce((i-1)*3+1:(i-1)*3+3) + cttforce((i-1)*3+1:(i-1)*3+3)
546  cont_nforce(idx:idx+2) = cont_nforce(idx:idx+2) + ctnforce((i-1)*3+1:(i-1)*3+3)
547  if( present(cont_fric) ) cont_fric(idx:idx+2) = cont_fric(idx:idx+2) + cttforce((i-1)*3+1:(i-1)*3+3)
548  enddo
549 
550  hecmat%B(np*ndof+id_lagrange) = ctnforce((nnode+1)*3+1)+cttforce((nnode+1)*3+1)
551 
552  end subroutine update_ndforce_contact
553 
554 
557  subroutine fstr_ass_load_contact(cstep, hecMESH, hecMAT, fstrSOLID, hecLagMAT)
558 
559  type(hecmwst_local_mesh) :: hecmesh
560  type(hecmwst_matrix) :: hecmat
561  type(fstr_solid) :: fstrsolid
562  type(hecmwst_matrix_lagrange) :: heclagmat
563  integer(kind=kint) :: cstep
564  integer(kind=kint) :: np, ndof
565  integer(kind=kint) :: i, j, k, l, id_lagrange, lnod, grpid
566  integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(21)
567  real(kind=kreal) :: ndu(21*3), ndcoord(21*3), lagrange
568  real(kind=kreal) :: normal(3), shapefunc(21)
569  real(kind=kreal) :: ntm(10*3)
570  real(kind=kreal) :: tf_final(3)
571  real(kind=kreal) :: ctforce(21*3 + 1)
572 
573  integer(kind=kint) :: algtype, nlag
574  real(kind=kreal) :: ctnforce(21*3+1)
575  real(kind=kreal) :: cttforce(21*3+1)
576  logical :: if_flag
577  real(kind=kreal) :: ctime, etime
578  integer(kind=kint) :: if_type
579 
580  np = hecmat%NP ; ndof = hecmat%NDOF
581 
582  id_lagrange = 0
583 
584  do i = 1, fstrsolid%n_contacts
585 
586  grpid = fstrsolid%contacts(i)%group
587  if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
588 
589  algtype = fstrsolid%contacts(i)%algtype
590  nlag = fstr_get_num_lagrange_pernode(algtype)
591 
592  if_flag = (fstrsolid%contacts(i)%if_type /= 0)
593  if(if_flag)then
594  ctime = fstrsolid%contacts(i)%ctime
595  etime = fstrsolid%contacts(i)%if_etime
596  if_type = fstrsolid%contacts(i)%if_type
597  end if
598 
599  do j = 1, size(fstrsolid%contacts(i)%slave)
600 
601  if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
602  if(if_flag) call set_shrink_factor(ctime, fstrsolid%contacts(i)%states(j), etime, if_type)
603 
604  ctsurf = fstrsolid%contacts(i)%states(j)%surface
605  etype = fstrsolid%contacts(i)%master(ctsurf)%etype
606  nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
607  ndlocal(1) = fstrsolid%contacts(i)%slave(j)
608  ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
609  do k = 1, nnode+1
610  ndu((k-1)*3+1:(k-1)*3+3) = fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
611  ndcoord((k-1)*3+1:(k-1)*3+3) = hecmesh%node((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) + ndu((k-1)*3+1:(k-1)*3+3)
612  enddo
613 
614  do k=1,nlag
615  id_lagrange = id_lagrange + 1
616  lagrange = fstrsolid%contacts(i)%states(j)%multiplier(k)
617  if(if_flag) call get_shrink_dummy_surf(fstrsolid%contacts(i)%states(j),ndcoord, nnode)
618  if( algtype == contactsslid .or. algtype == contactfslid ) then
619  ! Obtain contact nodal force vector of contact pair
620  call getcontactnodalforce(etype,nnode,ndcoord,ndu,fstrsolid%contacts(i)%states(j), &
621  fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,ctnforce,cttforce,.false.)
622  ! Update non-eqilibrited force vector
623  call update_ndforce_contact(nnode,ndlocal,id_lagrange,lagrange,ctnforce,cttforce, &
624  & hecmat,fstrsolid%CONT_NFORCE,fstrsolid%CONT_FRIC)
625  else if( algtype == contacttied ) then
626  call gettiednodalforce(etype,nnode,k,ndu,fstrsolid%contacts(i)%states(j),lagrange,ctnforce,cttforce)
627  ! Update non-eqilibrited force vector
628  call update_ndforce_contact(nnode,ndlocal,id_lagrange,-1.d0,ctnforce,cttforce,hecmat,fstrsolid%CONT_NFORCE)
629  endif
630 
631  enddo
632  enddo
633  enddo
634 
635  do i = 1, fstrsolid%n_embeds
636 
637  grpid = fstrsolid%embeds(i)%group
638  if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
639 
640 
641  do j = 1, size(fstrsolid%embeds(i)%slave)
642 
643  if( fstrsolid%embeds(i)%states(j)%state == contactfree ) cycle
644 
645  ctsurf = fstrsolid%embeds(i)%states(j)%surface
646  etype = fstrsolid%embeds(i)%master(ctsurf)%etype
647  nnode = size(fstrsolid%embeds(i)%master(ctsurf)%nodes)
648  ndlocal(1) = fstrsolid%embeds(i)%slave(j)
649  ndlocal(2:nnode+1) = fstrsolid%embeds(i)%master(ctsurf)%nodes(1:nnode)
650  do k = 1, nnode+1
651  ndu((k-1)*3+1:(k-1)*3+3) = fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
652  ndcoord((k-1)*3+1:(k-1)*3+3) = hecmesh%node((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) + ndu((k-1)*3+1:(k-1)*3+3)
653  enddo
654 
655  do k=1,3
656  id_lagrange = id_lagrange + 1
657  lagrange = fstrsolid%embeds(i)%states(j)%multiplier(k)
658 
659  call gettiednodalforce(etype,nnode,k,ndu,fstrsolid%embeds(i)%states(j),lagrange,ctnforce,cttforce)
660  ! Update non-eqilibrited force vector
661  call update_ndforce_contact(nnode,ndlocal,id_lagrange,-1.d0,ctnforce,cttforce,hecmat,fstrsolid%EMBED_NFORCE)
662 
663  enddo
664  enddo
665  enddo
666 
667  end subroutine fstr_ass_load_contact
668 
669  subroutine get_shrink_dummy_surf(cstate, coords, nnode)
670  use m_fstr_timeinc
671  integer(kind=kint) :: nnode, i
672  type(tcontactstate) :: cstate
673  real(kind=kreal) :: coords(:), ctime
674  real(kind=kreal) :: factor, d(3)
675 
676  d = - cstate%shrink_factor * cstate%direction(1:3)
677 
678  do i = 1, nnode
679  coords(1+i*3:(i+1)*3) = coords(1+i*3:(i+1)*3) + d(1:3)
680  enddo
681 
682  end subroutine get_shrink_dummy_surf
683 
684 end module m_addcontactstiffness
685 
686 
687 
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:647
This module provides functions of reconstructing.
integer(kind=kint) function, public fstr_get_num_lagrange_pernode(algtype)
subroutine, public hecmw_mat_ass_contactlag(nnode, ndLocal, id_lagrange, fcoeff, stiffness, hecMAT, hecLagMAT)
This subroutine assembles contact stiffness matrix of a contact pair into global stiffness matrix.
This module provides functions: 1) obtain contact stiffness matrix of each contact pair and assemble ...
subroutine gettiednodalforce(etype, nnode, idof, ndu, ctState, lagrange, ctNForce, ctTForce)
This subroutine obtains contact nodal force vector of contact pair.
subroutine, public fstr_ass_load_contact(cstep, hecMESH, hecMAT, fstrSOLID, hecLagMAT)
This subroutine adds initial contact force vector to the right-hand side vector \at the beginning of ...
subroutine, public update_ndforce_contact(nnode, ndLocal, id_lagrange, lagrange, ctNForce, ctTForce, hecMAT, cont_nforce, cont_fric)
This subroutine assembles contact nodal force vector into right-hand side vector to update non-equili...
subroutine get_shrink_dummy_surf(cstate, coords, nnode)
subroutine, public fstr_update_ndforce_contact(cstep, hecMESH, hecMAT, hecLagMAT, fstrSOLID, conMAT)
This subroutine obtains contact nodal force vector of each contact pair and assembles it into right-h...
subroutine, public fstr_addcontactstiffness(cstep, iter, hecMAT, hecLagMAT, fstrSOLID)
This subroutine obtains contact stiffness matrix of each contact pair and assembles it into global st...
subroutine gettiedstiffness(etype, nnode, idof, ctState, lagrange, stiffness)
This subroutine obtains contact stiffness matrix of contact pair.
This module provide functions of contact stiffness calculation.
Definition: contact_lib.f90:6
integer, parameter contacttied
contact type or algorithm definition
Definition: contact_lib.f90:21
integer, parameter contactfslid
Definition: contact_lib.f90:24
subroutine set_shrink_factor(ctime, cstate, etime, if_type)
integer, parameter contactslip
Definition: contact_lib.f90:18
integer, parameter contactsslid
Definition: contact_lib.f90:23
integer, parameter contactfree
contact state definition
Definition: contact_lib.f90:16
integer, parameter contactstick
Definition: contact_lib.f90:17
This module provides function to calculate residual of nodal force.
subroutine, public fstr_update_ndforce_spc(cstep, hecMESH, fstrSOLID, B)
This module provides functions to deal with time and increment of stress analysis.
This module defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
logical function fstr_isembedactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1063
real(kind=kreal) etime
Definition: m_fstr.f90:140
logical function fstr_iscontactactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1053
This structure records contact status.
Definition: contact_lib.f90:30