FrontISTR  5.7.0
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 
13 module m_addcontactstiffness
14 
15  use m_fstr
16  use elementinfo
17  use m_contact_lib
18  use fstr_matrix_con_contact
19  use hecmw_matrix_ass
20  use m_fstr_residual
21 
22  implicit none
23 
24  public :: fstr_addcontactstiffness
25  public :: fstr_update_ndforce_contact
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 
m_contact_lib::contactsslid
integer, parameter contactsslid
Definition: contact_lib.f90:23
m_addcontactstiffness
This module provides functions: 1) obtain contact stiffness matrix of each contact pair and assemble ...
Definition: fstr_AddContactStiff.f90:13
m_fstr_timeinc
This module provides functions to deal with time and increment of stress analysis.
Definition: fstr_Ctrl_TimeInc.f90:7
fstr_matrix_con_contact::fstr_get_num_lagrange_pernode
integer(kind=kint) function, public fstr_get_num_lagrange_pernode(algtype)
Definition: fstr_mat_con_contact.f90:32
m_contact_lib::contacttied
integer, parameter contacttied
contact type or algorithm definition
Definition: contact_lib.f90:21
m_fstr::fstr_solid
Definition: m_fstr.f90:239
m_fstr::fstr_iscontactactive
logical function fstr_iscontactactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1053
m_addcontactstiffness::fstr_update_ndforce_contact
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...
Definition: fstr_AddContactStiff.f90:257
m_addcontactstiffness::gettiedstiffness
subroutine gettiedstiffness(etype, nnode, idof, ctState, lagrange, stiffness)
This subroutine obtains contact stiffness matrix of contact pair.
Definition: fstr_AddContactStiff.f90:121
m_fstr_residual
This module provides function to calculate residual of nodal force.
Definition: fstr_Residual.f90:7
m_contact_lib::contactslip
integer, parameter contactslip
Definition: contact_lib.f90:18
m_contact_lib::contactfree
integer, parameter contactfree
contact state definition
Definition: contact_lib.f90:16
m_fstr_residual::fstr_update_ndforce_spc
subroutine, public fstr_update_ndforce_spc(cstep, hecMESH, fstrSOLID, B)
Definition: fstr_Residual.f90:98
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
m_contact_lib::contactstick
integer, parameter contactstick
Definition: contact_lib.f90:17
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
hecmw_matrix_ass
Definition: hecmw_mat_ass.f90:6
m_contact_lib
This module provide functions of contact stiffness calculation.
Definition: contact_lib.f90:6
fstr_matrix_con_contact
This module provides functions of reconstructing.
Definition: fstr_mat_con_contact.f90:9
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:647
hecmw_matrix_ass::hecmw_mat_ass_contactlag
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.
Definition: hecmw_mat_ass.f90:493
m_addcontactstiffness::update_ndforce_contact
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...
Definition: fstr_AddContactStiff.f90:527
m_addcontactstiffness::gettiednodalforce
subroutine gettiednodalforce(etype, nnode, idof, ndu, ctState, lagrange, ctNForce, ctTForce)
This subroutine obtains contact nodal force vector of contact pair.
Definition: fstr_AddContactStiff.f90:382
m_contact_lib::set_shrink_factor
subroutine set_shrink_factor(ctime, cstate, etime, if_type)
Definition: contact_lib.f90:506
m_addcontactstiffness::fstr_addcontactstiffness
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...
Definition: fstr_AddContactStiff.f90:38
m_addcontactstiffness::get_shrink_dummy_surf
subroutine get_shrink_dummy_surf(cstate, coords, nnode)
Definition: fstr_AddContactStiff.f90:670
m_fstr::etime
real(kind=kreal) etime
Definition: m_fstr.f90:140
m_fstr::fstr_isembedactive
logical function fstr_isembedactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1063
m_contact_lib::contactfslid
integer, parameter contactfslid
Definition: contact_lib.f90:24
m_addcontactstiffness::fstr_ass_load_contact
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 ...
Definition: fstr_AddContactStiff.f90:558
m_contact_lib::tcontactstate
This structure records contact status.
Definition: contact_lib.f90:30