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 
274  id_lagrange = 0
275  if( associated(fstrsolid%CONT_NFORCE) ) fstrsolid%CONT_NFORCE(:) = 0.d0
276  if( associated(fstrsolid%CONT_FRIC) ) fstrsolid%CONT_FRIC(:) = 0.d0
277  if( associated(fstrsolid%EMBED_NFORCE) ) fstrsolid%EMBED_NFORCE(:) = 0.d0
278 
279  do i = 1, fstrsolid%n_contacts
280 
281  grpid = fstrsolid%contacts(i)%group
282  if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
283 
284  algtype = fstrsolid%contacts(i)%algtype
285  nlag = fstr_get_num_lagrange_pernode(algtype)
286 
287  do j = 1, size(fstrsolid%contacts(i)%slave)
288 
289  if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
290 
291  ctsurf = fstrsolid%contacts(i)%states(j)%surface
292  etype = fstrsolid%contacts(i)%master(ctsurf)%etype
293  nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
294  ndlocal(1) = fstrsolid%contacts(i)%slave(j)
295  ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
296  do k = 1, nnode+1
297  nddu((k-1)*3+1:(k-1)*3+3) = fstrsolid%dunode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
298  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)
299  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)
300  enddo
301 
302  do k=1,nlag
303  id_lagrange = id_lagrange + 1
304  lagrange = heclagmat%Lagrange(id_lagrange)
305 
306  fstrsolid%contacts(i)%states(j)%multiplier(k) = heclagmat%Lagrange(id_lagrange)
307 
308  if( algtype == contactsslid .or. algtype == contactfslid ) then
309  ! Obtain contact nodal force vector of contact pair
310  call getcontactnodalforce(etype,nnode,ndcoord,nddu,fstrsolid%contacts(i)%states(j), &
311  fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,ctnforce,cttforce,.true.)
312  ! Update non-eqilibrited force vector
313  call update_ndforce_contact(nnode,ndlocal,id_lagrange,lagrange,ctnforce,cttforce, &
314  & conmat,fstrsolid%CONT_NFORCE,fstrsolid%CONT_FRIC)
315  else if( algtype == contacttied ) then
316  call gettiednodalforce(etype,nnode,k,ndu,fstrsolid%contacts(i)%states(j),lagrange,ctnforce,cttforce)
317  ! Update non-eqilibrited force vector
318  call update_ndforce_contact(nnode,ndlocal,id_lagrange,1.d0,ctnforce,cttforce, &
319  & conmat,fstrsolid%CONT_NFORCE,fstrsolid%CONT_FRIC)
320  endif
321 
322  enddo
323  enddo
324 
325  enddo
326 
327  do i = 1, fstrsolid%n_embeds
328 
329  grpid = fstrsolid%embeds(i)%group
330  if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
331 
332  do j = 1, size(fstrsolid%embeds(i)%slave)
333 
334  if( fstrsolid%embeds(i)%states(j)%state == contactfree ) cycle
335 
336  ctsurf = fstrsolid%embeds(i)%states(j)%surface
337  etype = fstrsolid%embeds(i)%master(ctsurf)%etype
338  nnode = size(fstrsolid%embeds(i)%master(ctsurf)%nodes)
339  ndlocal(1) = fstrsolid%embeds(i)%slave(j)
340  ndlocal(2:nnode+1) = fstrsolid%embeds(i)%master(ctsurf)%nodes(1:nnode)
341  do k = 1, nnode+1
342  nddu((k-1)*3+1:(k-1)*3+3) = fstrsolid%dunode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
343  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)
344  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)
345  enddo
346 
347  do k=1,3
348  id_lagrange = id_lagrange + 1
349  lagrange = heclagmat%Lagrange(id_lagrange)
350 
351  fstrsolid%embeds(i)%states(j)%multiplier(k) = heclagmat%Lagrange(id_lagrange)
352 
353  call gettiednodalforce(etype,nnode,k,ndu,fstrsolid%embeds(i)%states(j),lagrange,ctnforce,cttforce)
354  ! Update non-eqilibrited force vector
355  call update_ndforce_contact(nnode,ndlocal,id_lagrange,-1.d0,ctnforce,cttforce, &
356  & conmat,fstrsolid%EMBED_NFORCE,fstrsolid%EMBED_NFORCE)
357 
358  enddo
359  enddo
360  enddo
361 
362  ! Consider SPC condition
363  call fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, hecmat%B)
364  call fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, conmat%B)
365 
366  end subroutine fstr_update_ndforce_contact
367 
369  subroutine gettiednodalforce(etype,nnode,idof,ndu,ctState,lagrange,ctNForce,ctTForce)
370  integer(kind=kint) :: etype
371  integer(kind=kint) :: nnode
372  integer(kind=kint) :: idof
373  real(kind=kreal) :: ndu((nnode + 1)*3)
374  type(tcontactstate) :: ctstate
375  real(kind=kreal) :: lagrange
376  real(kind=kreal) :: ctnforce((nnode+1)*3+1)
377  real(kind=kreal) :: cttforce((nnode+1)*3+1)
378 
379  integer(kind=kint) :: i, j
380  real(kind=kreal) :: normal(3), shapefunc(nnode)
381  real(kind=kreal) :: ntm((nnode + 1)*3)
382  real(kind=kreal) :: ttm((nnode + 1)*3)
383 
384  ctnforce = 0.0d0
385  cttforce = 0.0d0
386 
387  call getshapefunc( etype, ctstate%lpos, shapefunc )
388 
389  normal(1:3) = ctstate%direction(1:3)
390 
391  ntm(1:3) = -normal(idof)*normal(1:3)
392  do i = 1, nnode
393  ntm(i*3+1:i*3+3) = -shapefunc(i)*ntm(1:3)
394  enddo
395  ttm = 0.d0
396  ttm(idof) = -1.0
397  ttm(1:3) = ttm(1:3)-ntm(1:3)
398  do i = 1, nnode
399  ttm(i*3+1:i*3+3) = -shapefunc(i)*ttm(1:3)
400  enddo
401 
402  do j = 1, (nnode+1)*3
403  ctnforce(j) = lagrange*ntm(j)
404  cttforce(j) = lagrange*ttm(j)
405  enddo
406  j = (nnode+1)*3 + 1
407  ctnforce(j) = dot_product(ntm,ndu)
408  cttforce(j) = dot_product(ttm,ndu)
409  end subroutine
410 
412  subroutine getcontactnodalforce(etype,nnode,ndCoord,ndDu,ctState,tPenalty,fcoeff,lagrange,ctNForce,ctTForce,cflag)
413 
414  type(tcontactstate) :: ctstate
415  integer(kind=kint) :: etype, nnode
416  integer(kind=kint) :: i, j
417  real(kind=kreal) :: fcoeff, tpenalty
418  real(kind=kreal) :: lagrange
419  real(kind=kreal) :: ndcoord((nnode + 1)*3), nddu((nnode + 1)*3)
420  real(kind=kreal) :: normal(3), shapefunc(nnode)
421  real(kind=kreal) :: ntm((nnode + 1)*3)
422  real(kind=kreal) :: tf_trial(3), length_tft, tangent(3), tf_final(3)
423  real(kind=kreal) :: ctnforce((nnode+1)*3+1)
424  real(kind=kreal) :: cttforce((nnode+1)*3+1)
425  logical :: cflag
426 
427  ctnforce = 0.0d0
428  cttforce = 0.0d0
429 
430  call getshapefunc( etype, ctstate%lpos, shapefunc )
431 
432  normal(1:3) = ctstate%direction(1:3)
433 
434  ntm(1:3) = -normal(1:3)
435  do i = 1, nnode
436  ntm(i*3+1:i*3+3) = shapefunc(i)*normal(1:3)
437  enddo
438 
439  do j = 1, (nnode+1)*3
440  ctnforce(j) = lagrange*ntm(j)
441  enddo
442  j = (nnode+1)*3 + 1
443  ctnforce(j) = dot_product(ntm,ndcoord)
444 
445  if(fcoeff /= 0.0d0 .and. lagrange > 0.0d0)then
446 
447  if( cflag ) then !calc tf_final and set it to tangentForce_final
448  call gettrialfricforceandcheckfricstate(nnode,tpenalty,fcoeff,lagrange,normal,shapefunc,ntm,nddu,ctstate)
449 
450  if( ctstate%state == contactstick ) then
451  tf_final(1:3) = ctstate%tangentForce_trial(1:3)
452  elseif( ctstate%state == contactslip ) then
453  tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
454  length_tft = dsqrt(dot_product(tf_trial,tf_trial))
455  tangent(1:3) = tf_trial(1:3)/length_tft
456  tf_final(1:3) = fcoeff*dabs(lagrange)*tangent(1:3)
457  endif
458  ctstate%tangentForce_final(1:3) = tf_final(1:3)
459  else ! just set tangentForce_final to tf_final (used for fstr_ass_load_contact)
460  tf_final(1:3) = ctstate%tangentForce_final(1:3)
461  endif
462 
463  cttforce(1:3) = -tf_final(1:3)
464  do j = 1, nnode
465  cttforce(j*3+1:j*3+3) = shapefunc(j)*tf_final(1:3)
466  enddo
467 
468  endif
469 
470  end subroutine getcontactnodalforce
471 
472 
474  subroutine gettrialfricforceandcheckfricstate(nnode,tPenalty,fcoeff,lagrange,normal,shapefunc,nTm,ndDu,ctstate)
475 
476  type(tcontactstate) :: ctState
477  integer(kind=kint) :: nnode
478  integer(kind=kint) :: i, j
479  real(kind=kreal) :: fcoeff, tpenalty
480  real(kind=kreal) :: lagrange
481  real(kind=kreal) :: nddu((nnode + 1)*3)
482  real(kind=kreal) :: normal(3), shapefunc(nnode)
483  real(kind=kreal) :: ntm((nnode + 1)*3)
484  real(kind=kreal) :: dotp
485  real(kind=kreal) :: relativedisp(3)
486  real(kind=kreal) :: tf_yield
487 
488  relativedisp = 0.0d0
489 
490  dotp = dot_product(ntm,nddu)
491  do i = 1, 3
492  relativedisp(i) = - nddu(i)
493  do j = 1, nnode
494  relativedisp(i) = relativedisp(i) + shapefunc(j)*nddu(j*3+i)
495  enddo
496  relativedisp(i) = relativedisp(i) - dotp*normal(i)
497  ctstate%reldisp(i) = -relativedisp(i)
498  ctstate%tangentForce_trial(i) = ctstate%tangentForce1(i) -tpenalty*relativedisp(i)
499  enddo
500 
501  tf_yield = fcoeff*dabs(lagrange)
502  if(ctstate%state == contactslip) tf_yield =0.99d0*tf_yield
503  if( dsqrt(dot_product(ctstate%tangentForce_trial,ctstate%tangentForce_trial)) <= tf_yield ) then
504  ctstate%state = contactstick
505  else
506  ctstate%state = contactslip
507  endif
508 
509  end subroutine gettrialfricforceandcheckfricstate
510 
511 
514  subroutine update_ndforce_contact(nnode,ndLocal,id_lagrange,lagrange,ctNForce,ctTForce,hecMAT,cont_nforce,cont_fric)
515 
516  type(hecmwst_matrix) :: hecmat
517  integer(kind=kint) :: nnode, ndlocal(nnode + 1)
519  integer(kind=kint) :: id_lagrange
520  real(kind=kreal) :: lagrange
521  integer(kind=kint) :: np, ndof
522  integer (kind=kint) :: i, inod, idx
523  real(kind=kreal) :: ctnforce((nnode+1)*3+1)
524  real(kind=kreal) :: cttforce((nnode+1)*3+1)
525  real(kind=kreal), pointer :: cont_nforce(:)
526  real(kind=kreal), pointer, optional :: cont_fric(:)
527 
528  np = hecmat%NP; ndof = hecmat%NDOF
529 
530  do i = 1, nnode + 1
531  inod = ndlocal(i)
532  idx = (inod-1)*3+1
533  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)
534  cont_nforce(idx:idx+2) = cont_nforce(idx:idx+2) + ctnforce((i-1)*3+1:(i-1)*3+3)
535  if( present(cont_fric) ) cont_fric(idx:idx+2) = cont_fric(idx:idx+2) + cttforce((i-1)*3+1:(i-1)*3+3)
536  enddo
537 
538  hecmat%B(np*ndof+id_lagrange) = ctnforce((nnode+1)*3+1)+cttforce((nnode+1)*3+1)
539 
540  end subroutine update_ndforce_contact
541 
542 
545  subroutine fstr_ass_load_contact(cstep, hecMESH, hecMAT, fstrSOLID, hecLagMAT)
546 
547  type(hecmwst_local_mesh) :: hecmesh
548  type(hecmwst_matrix) :: hecmat
549  type(fstr_solid) :: fstrsolid
550  type(hecmwst_matrix_lagrange) :: heclagmat
551  integer(kind=kint) :: cstep
552  integer(kind=kint) :: np, ndof
553  integer(kind=kint) :: i, j, k, l, id_lagrange, lnod, grpid
554  integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(21)
555  real(kind=kreal) :: ndu(21*3), ndcoord(21*3), lagrange
556  real(kind=kreal) :: normal(3), shapefunc(21)
557  real(kind=kreal) :: ntm(10*3)
558  real(kind=kreal) :: tf_final(3)
559  real(kind=kreal) :: ctforce(21*3 + 1)
560 
561  integer(kind=kint) :: algtype, nlag
562  real(kind=kreal) :: ctnforce(21*3+1)
563  real(kind=kreal) :: cttforce(21*3+1)
564 
565  np = hecmat%NP ; ndof = hecmat%NDOF
566 
567  id_lagrange = 0
568 
569  do i = 1, fstrsolid%n_contacts
570 
571  grpid = fstrsolid%contacts(i)%group
572  if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
573 
574  algtype = fstrsolid%contacts(i)%algtype
575  nlag = fstr_get_num_lagrange_pernode(algtype)
576 
577  do j = 1, size(fstrsolid%contacts(i)%slave)
578 
579  if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
580 
581  ctsurf = fstrsolid%contacts(i)%states(j)%surface
582  etype = fstrsolid%contacts(i)%master(ctsurf)%etype
583  nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
584  ndlocal(1) = fstrsolid%contacts(i)%slave(j)
585  ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
586  do k = 1, nnode+1
587  ndu((k-1)*3+1:(k-1)*3+3) = fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
588  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)
589  enddo
590 
591  do k=1,nlag
592  id_lagrange = id_lagrange + 1
593  lagrange = fstrsolid%contacts(i)%states(j)%multiplier(k)
594 
595  if( algtype == contactsslid .or. algtype == contactfslid ) then
596  ! Obtain contact nodal force vector of contact pair
597  call getcontactnodalforce(etype,nnode,ndcoord,ndu,fstrsolid%contacts(i)%states(j), &
598  fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,ctnforce,cttforce,.false.)
599  ! Update non-eqilibrited force vector
600  call update_ndforce_contact(nnode,ndlocal,id_lagrange,lagrange,ctnforce,cttforce, &
601  & hecmat,fstrsolid%CONT_NFORCE,fstrsolid%CONT_FRIC)
602  else if( algtype == contacttied ) then
603  call gettiednodalforce(etype,nnode,k,ndu,fstrsolid%contacts(i)%states(j),lagrange,ctnforce,cttforce)
604  ! Update non-eqilibrited force vector
605  call update_ndforce_contact(nnode,ndlocal,id_lagrange,-1.d0,ctnforce,cttforce,hecmat,fstrsolid%CONT_NFORCE)
606  endif
607 
608  enddo
609  enddo
610  enddo
611 
612  do i = 1, fstrsolid%n_embeds
613 
614  grpid = fstrsolid%embeds(i)%group
615  if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
616 
617 
618  do j = 1, size(fstrsolid%embeds(i)%slave)
619 
620  if( fstrsolid%embeds(i)%states(j)%state == contactfree ) cycle
621 
622  ctsurf = fstrsolid%embeds(i)%states(j)%surface
623  etype = fstrsolid%embeds(i)%master(ctsurf)%etype
624  nnode = size(fstrsolid%embeds(i)%master(ctsurf)%nodes)
625  ndlocal(1) = fstrsolid%embeds(i)%slave(j)
626  ndlocal(2:nnode+1) = fstrsolid%embeds(i)%master(ctsurf)%nodes(1:nnode)
627  do k = 1, nnode+1
628  ndu((k-1)*3+1:(k-1)*3+3) = fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
629  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)
630  enddo
631 
632  do k=1,3
633  id_lagrange = id_lagrange + 1
634  lagrange = fstrsolid%embeds(i)%states(j)%multiplier(k)
635 
636  call gettiednodalforce(etype,nnode,k,ndu,fstrsolid%embeds(i)%states(j),lagrange,ctnforce,cttforce)
637  ! Update non-eqilibrited force vector
638  call update_ndforce_contact(nnode,ndlocal,id_lagrange,-1.d0,ctnforce,cttforce,hecmat,fstrsolid%EMBED_NFORCE)
639 
640  enddo
641  enddo
642  enddo
643 
644  end subroutine fstr_ass_load_contact
645 
646 end module m_addcontactstiffness
647 
648 
649 
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
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:238
m_fstr::fstr_iscontactactive
logical function fstr_iscontactactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1052
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:515
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:370
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_fstr::fstr_isembedactive
logical function fstr_isembedactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1062
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:546
m_contact_lib::tcontactstate
This structure records contact status.
Definition: contact_lib.f90:27