FrontISTR  5.2.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 
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  public :: fstr_mat_ass_bc_contact
29 
30  private :: getcontactstiffness
31  private :: fstr_mat_ass_contact
32  private :: getcontactnodalforce
33  private :: gettrialfricforceandcheckfricstate
34 
35 contains
36 
39  subroutine fstr_addcontactstiffness(cstep,iter,hecMAT,fstrMAT,fstrSOLID)
40 
41  integer(kind=kint) :: cstep
42  integer(kind=kint) :: iter
43  type(hecmwst_matrix) :: hecmat
44  type(fstr_solid) :: fstrsolid
45  type(fstrst_matrix_contact_lagrange) :: fstrmat
46  integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(9)
47  integer(kind=kint) :: i, j, id_lagrange, grpid
48  real(kind=kreal) :: lagrange
49  real(kind=kreal) :: stiffness(9*3 + 1, 9*3 + 1)
50 
51  fstrmat%AL_lagrange = 0.0d0
52  fstrmat%AU_lagrange = 0.0d0
53 
54  id_lagrange = 0
55 
56  do i = 1, size(fstrsolid%contacts)
57 
58  grpid = fstrsolid%contacts(i)%group
59  if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
60 
61  do j = 1, size(fstrsolid%contacts(i)%slave)
62 
63  if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
64 
65  id_lagrange = id_lagrange + 1
66  lagrange = fstrmat%Lagrange(id_lagrange)
67 
68  ctsurf = fstrsolid%contacts(i)%states(j)%surface
69  etype = fstrsolid%contacts(i)%master(ctsurf)%etype
70  nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
71  ndlocal(1) = fstrsolid%contacts(i)%slave(j)
72  ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
73 
74  ! Obtain contact stiffness matrix of contact pair
75  call getcontactstiffness(iter,etype,nnode,fstrsolid%contacts(i)%states(j), &
76  fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,stiffness)
77 
78  ! Assemble contact stiffness matrix of contact pair into global stiffness matrix
79  call fstr_mat_ass_contact(nnode,ndlocal,id_lagrange,fstrsolid%contacts(i)%fcoeff,stiffness,hecmat,fstrmat)
80 
81  enddo
82 
83  enddo
84 
85 
86  end subroutine fstr_addcontactstiffness
87 
88 
90  subroutine getcontactstiffness(iter,etype,nnode,ctState,tPenalty,fcoeff,lagrange,stiffness)
91 
92  type(tcontactstate) :: ctstate
93  integer(kind=kint) :: iter
94  integer(kind=kint) :: etype, nnode
95  integer(kind=kint) :: i, j, k, l
96  real(kind=kreal) :: normal(3), shapefunc(nnode)
97  real(kind=kreal) :: ntm((nnode + 1)*3)
98  real(kind=kreal) :: fcoeff, tpenalty
99  real(kind=kreal) :: lagrange
100  real(kind=kreal) :: tf_trial(3), length_tft
101  real(kind=kreal) :: tangent(3), ttm((nnode + 1)*3)
102  real(kind=kreal) :: stiffness(9*3 + 1, 9*3 + 1)
103 
104  stiffness = 0.0d0
105 
106  call getshapefunc( etype, ctstate%lpos(:), shapefunc )
107 
108  normal(1:3) = ctstate%direction(1:3)
109 
110  ntm(1:3) = normal(1:3)
111  do i = 1, nnode
112  ntm(i*3+1:i*3+3) = -shapefunc(i)*normal(1:3)
113  enddo
114 
115  i = (nnode+1)*3 + 1
116  do j = 1, (nnode+1)*3
117  stiffness(i,j) = ntm(j); stiffness(j,i) = ntm(j)
118  enddo
119 
120 
121  if( fcoeff /= 0.0d0 ) then
122  if( lagrange>0.0d0 .or. iter==1 ) then
123 
124  do i = 1, nnode+1
125  do j = 1, i
126  do k = 1, 3
127  do l = 1, 3
128  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)
129  if( k==l ) then
130  if(i==1 .and. j==1)then
131  stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) + tpenalty
132  elseif(i>1 .and. j==1)then
133  stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) - tpenalty*shapefunc(i-1)
134  elseif(i>1 .and. j>1)then
135  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)
136  endif
137  endif
138  if(i==j) cycle
139  stiffness((j-1)*3+l,(i-1)*3+k) = stiffness((i-1)*3+k,(j-1)*3+l)
140  enddo
141  enddo
142  enddo
143  enddo
144 
145  if( ctstate%state == contactslip ) then
146 
147  tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
148  length_tft = dsqrt(dot_product(tf_trial,tf_trial))
149  tangent(1:3) = tf_trial(1:3)/length_tft
150 
151  ttm(1:3) = -tangent(1:3)
152  do i = 1, nnode
153  ttm(i*3+1:i*3+3) = shapefunc(i)*tangent(1:3)
154  enddo
155 
156  do i = 1, nnode+1
157  do j = 1, nnode+1
158  do k = 1, 3
159  do l = 1, 3
160  stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) &
161  + tpenalty*(-ttm((i-1)*3+k)*ttm((j-1)*3+l) &
162  +ttm((i-1)*3+k)*ntm((j-1)*3+l)*dot_product(tangent,normal))
163  enddo
164  enddo
165  enddo
166  enddo
167  stiffness(1:(nnode+1)*3,1:(nnode+1)*3) = (fcoeff*lagrange/length_tft)*stiffness(1:(nnode+1)*3,1:(nnode+1)*3)
168 
169  !
170  ! do i = 1, (nnode + 1)*3
171  ! do j = 1, i
172  ! do k = 1, 3
173  ! do l = 1, k
174  ! 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))
175  ! if( k/=l )stiffness((j-1)*3+l,(i-1)*3+k) = stiffness((i-1)*3+k,(j-1)*3+l)
176  ! enddo !- l
177  ! enddo !- k
178  ! enddo !- j
179  ! enddo !- i
180  ! stiffness(1:(nnode+1)*3,1:(nnode+1)*3) = (fcoeff*dabs(lagrange)/length_tft)*stiffness(1:(nnode+1)*3,1:(nnode+1)*3)
181 
182  ! j = (nnode+1)*3 + 1
183  ! do i = 1, (nnode+1)*3
184  ! stiffness(i,j) = stiffness(i,j) - fcoeff*tTm(i)
185  ! enddo
186  stiffness(1:(nnode+1)*3,(nnode+1)*3+1) = stiffness(1:(nnode+1)*3,(nnode+1)*3+1) - fcoeff*ttm(1:(nnode+1)*3)
187 
188  endif
189  endif
190  endif
191 
192  end subroutine getcontactstiffness
193 
194 
196  subroutine fstr_mat_ass_contact(nnode,ndLocal,id_lagrange,fcoeff,stiffness,hecMAT,fstrMAT)
197 
198  type(hecmwst_matrix) :: hecmat
199  type(fstrst_matrix_contact_lagrange) :: fstrmat
200  integer(kind=kint) :: nnode, ndlocal(nnode + 1), id_lagrange
203  integer(kind=kint) :: i, j, inod, jnod, l
204  integer(kind=kint) :: isl, iel, idxl_base, kl, idxl, isu, ieu, idxu_base, ku, idxu
205  real(kind=kreal) :: fcoeff
206  real(kind=kreal) :: stiffness(9*3 + 1, 9*3 + 1)
207  real(kind=kreal) :: a(3, 3)
208 
209  i = nnode + 1 + 1
210  inod = id_lagrange
211  isl = fstrmat%indexL_lagrange(inod-1)+1
212  iel = fstrmat%indexL_lagrange(inod)
213 
214  do j = 1, nnode + 1
215  jnod = ndlocal(j)
216  isu = fstrmat%indexU_lagrange(jnod-1)+1
217  ieu = fstrmat%indexU_lagrange(jnod)
218 
219  kl = hecmw_array_search_i(fstrmat%itemL_lagrange,isl,iel,jnod)
220  if( kl<isl .or. kl>iel ) then
221  write(*,*) '###ERROR### : cannot find connectivity (Lagrange1)'
222  stop
223  endif
224  ku = hecmw_array_search_i(fstrmat%itemU_lagrange,isu,ieu,inod)
225  if( ku<isu .or. ku>ieu ) then
226  write(*,*) '###ERROR### : cannot find connectivity (Lagrange2)'
227  stop
228  endif
229 
230  idxl_base = (kl-1)*3
231  idxu_base = (ku-1)*3
232 
233  do l = 1, 3
234  idxl = idxl_base + l
235  fstrmat%AL_lagrange(idxl) = fstrmat%AL_lagrange(idxl) + stiffness((i-1)*3+1,(j-1)*3+l)
236  idxu = idxu_base + l
237  fstrmat%AU_lagrange(idxu) = fstrmat%AU_lagrange(idxu) + stiffness((j-1)*3+l,(i-1)*3+1)
238  enddo
239  enddo
240 
241 
242  if(fcoeff /= 0.0d0)then
243 
244  do i = 1, nnode + 1
245  inod = ndlocal(i)
246  do j = 1, nnode + 1
247  jnod = ndlocal(j)
248  call stf_get_block(stiffness(1:(nnode+1)*3,1:(nnode+1)*3), 3, i, j, a)
249  call hecmw_mat_add_node(hecmat, inod, jnod, a)
250  enddo
251  enddo
252 
253  endif
254 
255  end subroutine fstr_mat_ass_contact
256 
257 
260  subroutine fstr_update_ndforce_contact(cstep,hecMESH,hecMAT,fstrMAT,fstrSOLID,conMAT)
261 
262  type(hecmwst_local_mesh) :: hecmesh
263  type(hecmwst_matrix) :: hecmat
264  type(fstr_solid) :: fstrsolid
265  type(fstrst_matrix_contact_lagrange) :: fstrmat
266  type(hecmwst_matrix), optional :: conmat
267  integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(9)
268  integer(kind=kint) :: i, j, k, id_lagrange
269  real(kind=kreal) :: ndcoord(9*3), nddu(9*3)
270  real(kind=kreal) :: lagrange
271  real(kind=kreal) :: ctnforce(9*3+1)
272  real(kind=kreal) :: cttforce(9*3+1)
273 
274  integer(kind=kint) :: cstep
275  integer(kind=kint) :: grpid
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 
281  do i = 1, size(fstrsolid%contacts)
282 
283  grpid = fstrsolid%contacts(i)%group
284  if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
285 
286  do j = 1, size(fstrsolid%contacts(i)%slave)
287 
288  if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
289 
290  id_lagrange = id_lagrange + 1
291  lagrange = fstrmat%Lagrange(id_lagrange)
292 
293  fstrsolid%contacts(i)%states(j)%multiplier(1) = fstrmat%Lagrange(id_lagrange)
294 
295  ctsurf = fstrsolid%contacts(i)%states(j)%surface
296  etype = fstrsolid%contacts(i)%master(ctsurf)%etype
297  nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
298  ndlocal(1) = fstrsolid%contacts(i)%slave(j)
299  ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
300  do k = 1, nnode+1
301  ndcoord((k-1)*3+1:(k-1)*3+3) = hecmesh%node((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) &
302  + fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) &
303  + fstrsolid%dunode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
304  nddu((k-1)*3+1:(k-1)*3+3) = fstrsolid%dunode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
305  enddo
306  ! Obtain contact nodal force vector of contact pair
307  call getcontactnodalforce(etype,nnode,ndcoord,nddu,fstrsolid%contacts(i)%states(j), &
308  fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,ctnforce,cttforce)
309  ! Update non-eqilibrited force vector
310  if(present(conmat)) then
311  call update_ndforce_contact(nnode,ndlocal,id_lagrange,lagrange,ctnforce,cttforce,fstrsolid,conmat)
312  else
313  call update_ndforce_contact(nnode,ndlocal,id_lagrange,lagrange,ctnforce,cttforce,fstrsolid,hecmat)
314  endif
315 
316  enddo
317 
318  enddo
319 
320  ! Consider SPC condition
321  call fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, hecmat%B)
322  if(present(conmat)) call fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, conmat%B)
323 
324  end subroutine fstr_update_ndforce_contact
325 
327  subroutine getcontactnodalforce(etype,nnode,ndCoord,ndDu,ctState,tPenalty,fcoeff,lagrange,ctNForce,ctTForce)
328 
329  type(tcontactstate) :: ctstate
330  integer(kind=kint) :: etype, nnode
331  integer(kind=kint) :: i, j
332  real(kind=kreal) :: fcoeff, tpenalty
333  real(kind=kreal) :: lagrange
334  real(kind=kreal) :: ndcoord((nnode + 1)*3), nddu((nnode + 1)*3)
335  real(kind=kreal) :: normal(3), shapefunc(nnode)
336  real(kind=kreal) :: ntm((nnode + 1)*3)
337  real(kind=kreal) :: tf_trial(3), length_tft, tangent(3), tf_final(3)
338  real(kind=kreal) :: ctnforce((nnode+1)*3+1)
339  real(kind=kreal) :: cttforce((nnode+1)*3+1)
340 
341  ctnforce = 0.0d0
342  cttforce = 0.0d0
343 
344  call getshapefunc( etype, ctstate%lpos(:), shapefunc )
345 
346  normal(1:3) = ctstate%direction(1:3)
347 
348  ntm(1:3) = -normal(1:3)
349  do i = 1, nnode
350  ntm(i*3+1:i*3+3) = shapefunc(i)*normal(1:3)
351  enddo
352 
353  do j = 1, (nnode+1)*3
354  ctnforce(j) = lagrange*ntm(j)
355  enddo
356  j = (nnode+1)*3 + 1
357  ctnforce(j) = dot_product(ntm,ndcoord)
358 
359  if(fcoeff /= 0.0d0 .and. lagrange > 0.0d0)then
360 
361  call gettrialfricforceandcheckfricstate(nnode,tpenalty,fcoeff,lagrange,normal,shapefunc,ntm,nddu,ctstate)
362 
363  if( ctstate%state == contactstick ) then
364  tf_final(1:3) = ctstate%tangentForce_trial(1:3)
365  elseif( ctstate%state == contactslip ) then
366  tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
367  length_tft = dsqrt(dot_product(tf_trial,tf_trial))
368  tangent(1:3) = tf_trial(1:3)/length_tft
369  tf_final(1:3) = fcoeff*dabs(lagrange)*tangent(1:3)
370  endif
371 
372  cttforce(1:3) = -tf_final(1:3)
373  do j = 1, nnode
374  cttforce(j*3+1:j*3+3) = shapefunc(j)*tf_final(1:3)
375  enddo
376 
377  ctstate%tangentForce_final(1:3) = tf_final(1:3)
378 
379  endif
380 
381  end subroutine getcontactnodalforce
382 
383 
385  subroutine gettrialfricforceandcheckfricstate(nnode,tPenalty,fcoeff,lagrange,normal,shapefunc,nTm,ndDu,ctstate)
386 
387  type(tcontactstate) :: ctstate
388  integer(kind=kint) :: nnode
389  integer(kind=kint) :: i, j
390  real(kind=kreal) :: fcoeff, tpenalty
391  real(kind=kreal) :: lagrange
392  real(kind=kreal) :: nddu((nnode + 1)*3)
393  real(kind=kreal) :: normal(3), shapefunc(nnode)
394  real(kind=kreal) :: ntm((nnode + 1)*3)
395  real(kind=kreal) :: dotp
396  real(kind=kreal) :: relativedisp(3)
397  real(kind=kreal) :: tf_yield
398 
399  relativedisp = 0.0d0
400 
401  dotp = dot_product(ntm,nddu)
402  do i = 1, 3
403  relativedisp(i) = - nddu(i)
404  do j = 1, nnode
405  relativedisp(i) = relativedisp(i) + shapefunc(j)*nddu(j*3+i)
406  enddo
407  relativedisp(i) = relativedisp(i) - dotp*normal(i)
408  ctstate%reldisp(i) = -relativedisp(i)
409  ctstate%tangentForce_trial(i) = ctstate%tangentForce1(i) -tpenalty*relativedisp(i)
410  enddo
411 
412  tf_yield = fcoeff*dabs(lagrange)
413  if(ctstate%state == contactslip) tf_yield =0.99d0*tf_yield
414  if( dsqrt(dot_product(ctstate%tangentForce_trial,ctstate%tangentForce_trial)) <= tf_yield ) then
415  ctstate%state = contactstick
416  else
417  ctstate%state = contactslip
418  endif
419 
420  end subroutine gettrialfricforceandcheckfricstate
421 
422 
425  subroutine update_ndforce_contact(nnode,ndLocal,id_lagrange,lagrange,ctNForce,ctTForce,fstrSOLID,hecMAT)
426 
427  type(fstr_solid) :: fstrsolid
428  type(hecmwst_matrix) :: hecmat
429  integer(kind=kint) :: nnode, ndlocal(nnode + 1)
431  integer(kind=kint) :: id_lagrange
432  real(kind=kreal) :: lagrange
433  integer(kind=kint) :: np, ndof
434  integer (kind=kint) :: i, inod, idx
435  real(kind=kreal) :: ctnforce((nnode+1)*3+1)
436  real(kind=kreal) :: cttforce((nnode+1)*3+1)
437 
438  np = hecmat%NP; ndof = hecmat%NDOF
439 
440  do i = 1, nnode + 1
441  inod = ndlocal(i)
442  idx = (inod-1)*3+1
443  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)
444  if( lagrange < 0.d0 ) cycle
445  fstrsolid%CONT_NFORCE(idx:idx+2) = fstrsolid%CONT_NFORCE(idx:idx+2) + ctnforce((i-1)*3+1:(i-1)*3+3)
446  fstrsolid%CONT_FRIC(idx:idx+2) = fstrsolid%CONT_FRIC(idx:idx+2) + cttforce((i-1)*3+1:(i-1)*3+3)
447  enddo
448 
449  hecmat%B(np*ndof+id_lagrange) = ctnforce((nnode+1)*3+1)+cttforce((nnode+1)*3+1)
450 
451  end subroutine update_ndforce_contact
452 
453 
456  subroutine fstr_ass_load_contact(cstep, hecMESH, hecMAT, fstrSOLID, fstrMAT)
457 
458  type(hecmwst_local_mesh) :: hecmesh
459  type(hecmwst_matrix) :: hecmat
460  type(fstr_solid) :: fstrsolid
461  type(fstrst_matrix_contact_lagrange) :: fstrmat
462  integer(kind=kint) :: cstep
463  integer(kind=kint) :: np, ndof
464  integer(kind=kint) :: i, j, k, l, id_lagrange, lnod, grpid
465  integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(9)
466  real(kind=kreal) :: ndcoord(9*3), lagrange
467  real(kind=kreal) :: normal(3), shapefunc(9)
468  real(kind=kreal) :: ntm(10*3)
469  real(kind=kreal) :: tf_final(3)
470  real(kind=kreal) :: ctforce(9*3 + 1)
471 
472  np = hecmat%NP ; ndof = hecmat%NDOF
473 
474  id_lagrange = 0
475 
476  do i = 1, size(fstrsolid%contacts)
477 
478  grpid = fstrsolid%contacts(i)%group
479  if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
480 
481  do j = 1, size(fstrsolid%contacts(i)%slave)
482 
483  if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
484 
485  id_lagrange = id_lagrange + 1
486  lagrange = fstrsolid%contacts(i)%states(j)%multiplier(1)
487 
488  ctsurf = fstrsolid%contacts(i)%states(j)%surface
489  etype = fstrsolid%contacts(i)%master(ctsurf)%etype
490  nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
491  ndlocal(1) = fstrsolid%contacts(i)%slave(j)
492  ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
493  do k = 1, nnode+1
494  ndcoord((k-1)*3+1:(k-1)*3+3) = hecmesh%node((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) &
495  + fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
496  enddo
497 
498  ctforce = 0.0d0
499 
500  call getshapefunc( etype, fstrsolid%contacts(i)%states(j)%lpos(:), shapefunc )
501  normal(1:3) = fstrsolid%contacts(i)%states(j)%direction(1:3)
502  ntm(1:3) = -normal(1:3)
503  do k = 1, nnode
504  ntm(k*3+1:k*3+3) = shapefunc(k)*normal(1:3)
505  enddo
506  do l = 1, (nnode+1)*3
507  ctforce(l) = lagrange*ntm(l)
508  enddo
509  l = (nnode+1)*3 + 1
510  ctforce(l) = dot_product(ntm(1:(nnode+1)*3),ndcoord(1:(nnode+1)*3))
511 
512  if( fstrsolid%contacts(i)%fcoeff/=0.0d0 .and. lagrange>0.0d0 )then
513  tf_final(1:3) = fstrsolid%contacts(i)%states(j)%tangentForce_final(1:3)
514  ctforce(1:3) = ctforce(1:3) - tf_final(1:3)
515  do l = 1, nnode
516  ctforce(l*3+1:l*3+3) = ctforce(l*3+1:l*3+3) + shapefunc(l)*tf_final(1:3)
517  enddo
518  endif
519 
520  do l = 1, nnode + 1
521  lnod = ndlocal(l)
522  hecmat%B((lnod-1)*3+1:(lnod-1)*3+3) = hecmat%B((lnod-1)*3+1:(lnod-1)*3+3) + ctforce((l-1)*3+1:(l-1)*3+3)
523  enddo
524  hecmat%B(np*ndof+id_lagrange) = ctforce((nnode+1)*3+1)
525 
526  enddo
527 
528  enddo
529 
530  end subroutine fstr_ass_load_contact
531 
532 
535  subroutine fstr_mat_ass_bc_contact(hecMAT,fstrMAT,inode,idof,RHS)
536 
537  type(hecmwst_matrix) :: hecmat
538  type(fstrst_matrix_contact_lagrange) :: fstrmat
539  integer(kind=kint) :: inode, idof
540  integer(kind=kint) :: isu, ieu, isl, iel, i, l, k
541  real(kind=kreal) :: rhs
542 
543  isu = fstrmat%indexU_lagrange(inode-1)+1
544  ieu = fstrmat%indexU_lagrange(inode)
545  do i = isu, ieu
546  fstrmat%AU_lagrange((i-1)*3+idof) = 0.0d0
547  l = fstrmat%itemU_lagrange(i)
548  isl = fstrmat%indexL_lagrange(l-1)+1
549  iel = fstrmat%indexL_lagrange(l)
550  k = hecmw_array_search_i(fstrmat%itemL_lagrange,isl,iel,inode)
551  if(k < isl .or. k > iel) cycle
552  hecmat%B(hecmat%NP*hecmat%NDOF+l) = hecmat%B(hecmat%NP*hecmat%NDOF+l) - fstrmat%AL_lagrange((k-1)*3+idof)*rhs
553  fstrmat%AL_lagrange((k-1)*3+idof) = 0.0d0
554  enddo
555 
556  end subroutine fstr_mat_ass_bc_contact
557 
558 end module m_addcontactstiffness
559 
560 
561 
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 coodinate system.
Definition: element.f90:640
This module provides functions of reconstructing.
integer(kind=kint) function, public hecmw_array_search_i(array, is, iE, ival)
subroutine, public stf_get_block(stiffness, ndof, inod, jnod, a)
subroutine, public hecmw_mat_add_node(hecMAT, inod, jnod, a)
This module provides functions: 1) obtain contact stiffness matrix of each contact pair and assemble ...
subroutine, public fstr_update_ndforce_contact(cstep, hecMESH, hecMAT, fstrMAT, fstrSOLID, conMAT)
This subroutine obtains contact nodal force vector of each contact pair and assembles it into right-h...
subroutine, public fstr_mat_ass_bc_contact(hecMAT, fstrMAT, inode, idof, RHS)
Modify Lagrange multiplier-related part of stiffness matrix and right-hand side vector for dealing wi...
subroutine, public fstr_ass_load_contact(cstep, hecMESH, hecMAT, fstrSOLID, fstrMAT)
This subroutine adds initial contact force vector to the right-hand side vector \at the beginning of ...
subroutine, public fstr_addcontactstiffness(cstep, iter, hecMAT, fstrMAT, fstrSOLID)
This subroutine obtains contact stiffness matrix of each contact pair and assembles it into global st...
subroutine, public update_ndforce_contact(nnode, ndLocal, id_lagrange, lagrange, ctNForce, ctTForce, fstrSOLID, hecMAT)
This subroutine assembles contact nodal force vector into right-hand side vector to update non-equili...
This module provide functions of contact stiffness calculation.
Definition: contact_lib.f90:6
integer, parameter contactslip
Definition: contact_lib.f90:18
integer, parameter contactfree
contact state definition
Definition: contact_lib.f90:16
integer, parameter contactstick
Definition: contact_lib.f90:17
This module provides function to calcualte residual of nodal force.
subroutine, public fstr_update_ndforce_spc(cstep, hecMESH, fstrSOLID, B)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
logical function fstr_iscontactactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1021
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...
This structure records contact status.
Definition: contact_lib.f90:27