FrontISTR  5.8.0
Large-scale structural analysis program with finit element method
contact_lib.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 !-------------------------------------------------------------------------------
6 module m_contact_lib
7  use elementinfo
8  implicit none
9 
10  integer, parameter, private :: kreal = kind(0.0d0)
11  integer, parameter, private :: l_max_surface_node = 20
12  integer, parameter, private :: l_max_elem_node = 100
13 
14  integer, parameter :: contactunknown = -1
16  integer, parameter :: contactfree = -1
17  integer, parameter :: contactstick = 1
18  integer, parameter :: contactslip = 2
19 
21  integer, parameter :: contacttied = 1
22  integer, parameter :: contactglued = 2
23  integer, parameter :: contactsslid = 3
24  integer, parameter :: contactfslid = 4
25 
27  integer, parameter :: c_if_slave = 1
28  integer, parameter :: c_if_master = 2
30  type tcontactstate
31  integer :: state
32  integer :: surface
33  real(kind=kreal) :: distance
34  real(kind=kreal) :: wkdist
35  real(kind=kreal) :: lpos(3)
36  real(kind=kreal) :: gpos(3)
37  real(kind=kreal) :: direction(3)
38  real(kind=kreal) :: multiplier(3)
40  real(kind=kreal) :: tangentforce(3)
41  real(kind=kreal) :: tangentforce1(3)
42  real(kind=kreal) :: tangentforce_trial(3)
43  real(kind=kreal) :: tangentforce_final(3)
44  real(kind=kreal) :: reldisp(3)
45  !
46  real(kind=kreal) :: shrink_factor
47  real(kind=kreal) :: time_factor
48  real(kind=kreal) :: init_pos
49  real(kind=kreal) :: end_pos
50  integer :: interference_flag
51  end type
52 
53 contains
54 
56  subroutine contact_state_init(cstate)
57  type(tcontactstate), intent(inout) :: cstate
58  cstate%state = -1
59  cstate%surface = -1
60  cstate%distance = 0.0d0
61  cstate%wkdist = 0.0d0
62  cstate%lpos(:) = 0.0d0
63  cstate%gpos(:) = 0.0d0
64  cstate%direction(:) = 0.0d0
65  cstate%multiplier(:) = 0.0d0
66  cstate%tangentForce(:) = 0.0d0
67  cstate%tangentForce1(:) = 0.0d0
68  cstate%tangentForce_trial(:) = 0.0d0
69  cstate%tangentForce_final(:) = 0.0d0
70  cstate%reldisp(:) = 0.0d0
71  cstate%shrink_factor = 0.0d0
72  cstate%time_factor = 0.0d0
73  cstate%init_pos = 0.0d0
74  cstate%end_pos = 0.0d0
75  cstate%interference_flag = 0
76  end subroutine
77 
79  subroutine contact_state_copy(cstate1, cstate2)
80  type(tcontactstate), intent(in) :: cstate1
81  type(tcontactstate), intent(inout) :: cstate2
82  cstate2 = cstate1
83  end subroutine
84 
86  subroutine print_contact_state(fnum, cstate)
87  integer, intent(in) :: fnum
88  type(tcontactstate), intent(in) :: cstate
89  write(fnum, *) "--Contact state=",cstate%state
90  write(fnum, *) cstate%surface, cstate%distance
91  write(fnum, *) cstate%lpos
92  write(fnum, *) cstate%direction
93  write(fnum, *) cstate%multiplier
94  end subroutine
95 
97  subroutine contact2mpcval( cstate, etype, nnode, mpcval )
98  type(tcontactstate), intent(in) :: cstate
99  integer, intent(in) :: etype
100  integer, intent(in) :: nnode
101  real(kind=kreal), intent(out) :: mpcval(nnode*3 + 4)
102 
103  integer :: i,j
104  real(kind=kreal) :: shapefunc(nnode)
105 
106  call getshapefunc( etype, cstate%lpos(1:2), shapefunc )
107  mpcval(1:3) = cstate%direction(1:3)
108  do i=1,nnode
109  do j=1,3
110  mpcval( i*3+j ) = -cstate%direction(j)*shapefunc(i)
111  enddo
112  enddo
113  mpcval( 3*nnode+4 )=cstate%distance
114  end subroutine
115 
117  subroutine contact2stiff( flag, cstate, etype, nnode, ele, mu, mut, &
118  fcoeff, symm, stiff, force )
119  integer, intent(in) :: flag
120  type(tcontactstate), intent(in) :: cstate
121  integer, intent(in) :: etype
122  integer, intent(in) :: nnode
123  real(kind=kreal), intent(in) :: ele(3,nnode)
124  real(kind=kreal), intent(in) :: mu
125  real(kind=kreal), intent(in) :: mut
126  real(kind=kreal), intent(in) :: fcoeff
127  logical, intent(in) :: symm
128  real(kind=kreal), intent(out) :: stiff(:,:)
129  real(kind=kreal), intent(out) :: force(:)
130 
131  integer :: i,j
132  real(kind=kreal) :: shapefunc(nnode)
133  real(kind=kreal) :: n(nnode*3+3), dispmat(2,nnode*3+3)
134  real(kind=kreal) :: metric(2,2)
135  real(kind=kreal) :: det, inverse(2,2), ff(2), cff(2)
136  real(kind=kreal) :: dum11(nnode*3+3,nnode*3+3), dum12(nnode*3+3,nnode*3+3)
137  real(kind=kreal) :: dum21(nnode*3+3,nnode*3+3), dum22(nnode*3+3,nnode*3+3)
138  real(kind=kreal) :: tangent(3,2)
139 
140  call getshapefunc( etype, cstate%lpos(1:2), shapefunc )
141  n(1:3) = cstate%direction(1:3)
142  do i=1,nnode
143  n( i*3+1:i*3+3 ) = -shapefunc(i)*cstate%direction(1:3)
144  enddo
145  do j=1,nnode*3+3
146  do i=1,nnode*3+3
147  stiff(i,j) = mu* n(i)*n(j)
148  enddo
149  enddo
150  force(1:nnode*3+3) = n(:)
151 
152  if( fcoeff/=0.d0 .or. flag==contactfslid ) &
153  call dispincrematrix( cstate%lpos(1:2), etype, nnode, ele, tangent, metric, dispmat )
154 
155  ! frictional component
156  if( fcoeff/=0.d0 ) then
157  forall(i=1:nnode*3+3, j=1:nnode*3+3)
158  dum11(i,j) = mut*dispmat(1,i)*dispmat(1,j)
159  dum12(i,j) = mut*dispmat(1,i)*dispmat(2,j)
160  dum21(i,j) = mut*dispmat(2,i)*dispmat(1,j)
161  dum22(i,j) = mut*dispmat(2,i)*dispmat(2,j)
162  end forall
163  stiff(1:nnode*3+3,1:nnode*3+3) &
164  = stiff(1:nnode*3+3,1:nnode*3+3) &
165  + metric(1,1)*dum11 + metric(1,2)*dum12 &
166  + metric(2,1)*dum21 + metric(2,2)*dum22
167 
168  if( cstate%state == contactslip ) then
169  det = metric(1,1)*metric(2,2)-metric(1,2)*metric(2,1)
170  if( det==0.d0 ) stop "Math error in contact stiff calculation"
171  inverse(1,1) = metric(2,2)/det
172  inverse(2,1) = -metric(2,1)/det
173  inverse(1,2) = -metric(1,2)/det
174  inverse(2,2) = metric(1,1)/det
175  ff(:) = cstate%multiplier(2:3)
176  cff(:) = matmul( inverse, ff )
177  ff(:) = ff(:)/dsqrt( ff(1)*ff(1)+ff(2)*ff(2) )
178  cff(:) = cff(:)/dsqrt( cff(1)*cff(1)+cff(2)*cff(2) )
179  stiff(1:nnode*3+3,1:nnode*3+3) = stiff(1:nnode*3+3,1:nnode*3+3) - &
180  ( cff(1)*ff(1)*metric(1,1)+ cff(2)*ff(1)*metric(1,2) )*dum11 - &
181  ( cff(2)*ff(2)*metric(1,2)+ cff(1)*ff(2)*metric(1,1) )*dum21 - &
182  ( cff(1)*ff(1)*metric(1,2)+ cff(2)*ff(1)*metric(2,2) )*dum12 - &
183  ( cff(2)*ff(2)*metric(2,2)+ cff(1)*ff(2)*metric(1,2) )*dum22
184  endif
185  endif
186 
187  end subroutine
188 
190  subroutine tied2stiff( flag, cstate, etype, nnode, mu, mut, stiff, force )
191  integer, intent(in) :: flag
192  type(tcontactstate), intent(in) :: cstate
193  integer, intent(in) :: etype
194  integer, intent(in) :: nnode
195  real(kind=kreal), intent(in) :: mu
196  real(kind=kreal), intent(in) :: mut
197  real(kind=kreal), intent(out) :: stiff(:,:)
198  real(kind=kreal), intent(out) :: force(:)
199 
200  integer :: i,j,k
201  real(kind=kreal) :: shapefunc(nnode)
202  real(kind=kreal) :: n(nnode*3+3)
203 
204  stiff = 0.d0
205 
206  call getshapefunc( etype, cstate%lpos(1:2), shapefunc )
207  n(1) = 1.d0
208  n(2:nnode+1) = -shapefunc(1:nnode)
209 
210  do j=1,nnode+1
211  do k=1,nnode+1
212  do i=1,3
213  stiff(3*k-3+i,3*j-3+i) = mu*n(k)*n(j)
214  enddo
215  enddo
216  enddo
217  force(1:nnode*3+3) = n(:)
218 
219 end subroutine
220 
222  subroutine getmetrictensor( pos, etype, ele, tensor )
223  real(kind=kreal), intent(in) :: pos(2)
224  integer, intent(in) :: etype
225  real(kind=kreal), intent(in) :: ele(:,:)
226  real(kind=kreal), intent(out) :: tensor(2,2)
227 
228  integer :: nn
229  real(kind=kreal) :: tangent(3,2)
230  nn= getnumberofnodes(etype)
231  call tangentbase( etype, nn, pos, ele, tangent )
232  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
233  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
234  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
235  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
236  end subroutine
237 
240  subroutine dispincrematrix( pos, etype, nnode, ele, tangent, tensor, matrix )
241  real(kind=kreal), intent(in) :: pos(2)
242  integer, intent(in) :: etype
243  integer, intent(in) :: nnode
244  real(kind=kreal), intent(in) :: ele(3,nnode)
245  real(kind=kreal), intent(out) :: tangent(3,2)
246  real(kind=kreal), intent(out) :: tensor(2,2)
247  real(kind=kreal), intent(out) :: matrix(2,nnode*3+3)
248 
249  integer :: i,j
250  real(kind=kreal) :: det
251  real(kind=kreal) :: shapefunc(nnode), t1(nnode*3+3), t2(nnode*3+3)
252  call tangentbase( etype, nnode, pos, ele, tangent )
253  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
254  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
255  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
256  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
257  det = tensor(1,1)*tensor(2,2)-tensor(1,2)*tensor(2,1)
258  if( det==0.d0 ) stop "Error in calculate DispIncreMatrix"
259  ! inverse(1,1) = tensor(2,2)/det
260  ! inverse(1,2) = -tensor(1,2)/det
261  ! inverse(2,1) = -tensor(2,1)/det
262  ! inverse(2,2) = tensor(1,1)/det
263 
264  call getshapefunc( etype, pos(:), shapefunc )
265  forall( j=1:3 )
266  t1( j ) = tangent(j,1)
267  t2( j ) = tangent(j,2)
268  end forall
269  forall( i=1:nnode, j=1:3 )
270  t1( i*3+j ) = -tangent(j,1)*shapefunc(i)
271  t2( i*3+j ) = -tangent(j,2)*shapefunc(i)
272  end forall
273  !matrix( 1:2,: ) = matmul( inverse(:,:), matrix )
274  matrix(1,:) = (tensor(2,2)*t1(:)-tensor(1,2)*t2(:))/det
275  matrix(2,:) = (tensor(1,1)*t2(:)-tensor(2,1)*t1(:))/det
276  tangent(:,1) = tangent(:,1)/dsqrt(dot_product(tangent(:,1),tangent(:,1)))
277  tangent(:,2) = tangent(:,2)/dsqrt(dot_product(tangent(:,2),tangent(:,2)))
278  end subroutine
279 
281  subroutine project_point2element(xyz,etype,nn,elemt,reflen,cstate,isin,distclr,ctpos,localclr)
282  real(kind=kreal),intent(in) :: xyz(3)
283  integer, intent(in) :: etype
284  integer, intent(in) :: nn
285  real(kind=kreal),intent(in) :: elemt(:,:)
286  real(kind=kreal),intent(in) :: reflen
287  type(tcontactstate),intent(inout) :: cstate
288  logical, intent(out) :: isin
289  real(kind=kreal), intent(in) :: distclr
290  real(kind=kreal), optional :: ctpos(2)
291  real(kind=kreal), optional :: localclr
292 
293  integer :: count,order, initstate
294  real(kind=kreal) :: determ, inverse(2,2)
295  real(kind=kreal) :: sfunc(nn), curv(3,2,2)
296  real(kind=kreal) :: r(2), dr(2), r_tmp(2) ! natural coordinate
297  real(kind=kreal) :: xyz_out(3) ! curr. projection position
298  real(kind=kreal) :: dist_last,dist_now, dxyz(3) ! dist between the point and its projection
299  real(kind=kreal) :: tangent(3,2) ! base vectors in tangent space
300  real(kind=kreal) :: df(2),d2f(2,2),normal(3)
301  real(kind=kreal),parameter :: eps = 1.0d-8
302  real(kind=kreal) :: clr, tol, factor
303 
304  initstate = cstate%state
305  clr = 1.d-4
306  if( present( localclr ) ) clr=localclr
307  if( present( ctpos ) ) then
308  r(:)= ctpos
309  else
310  call getelementcenter( etype, r(:) )
311  endif
312 
313  tol = 1.0d0
314  do count=1,100
315  call getshapefunc( etype, r, sfunc )
316  xyz_out = matmul( elemt(1:3,1:nn), sfunc )
317  dxyz(1:3) = xyz_out(1:3) - xyz(1:3)
318  dist_last = dot_product( dxyz, dxyz(:) )
319 
320  call tangentbase( etype, nn, r, elemt(1:3,1:nn), tangent )
321  call curvature( etype, nn, r, elemt(1:3,1:nn), curv )
322 
323  ! dF(1:2)
324  df(1:2) = -matmul( dxyz(:), tangent(:,:) )
325  ! d2F(1:2,1:2)
326  d2f(1,1)= dot_product( tangent(:,1), tangent(:,1) ) - dot_product( dxyz, curv(:,1,1) )
327  d2f(1,2)= dot_product( tangent(:,1), tangent(:,2) ) - dot_product( dxyz, curv(:,1,2) )
328  d2f(2,1)= dot_product( tangent(:,2), tangent(:,1) ) - dot_product( dxyz, curv(:,2,1) )
329  d2f(2,2)= dot_product( tangent(:,2), tangent(:,2) ) - dot_product( dxyz, curv(:,2,2) )
330 
331  ! inverse of d2F
332  determ = d2f(1,1)*d2f(2,2) - d2f(1,2)*d2f(2,1)
333  if( determ==0.d0 ) stop "Math error in contact searching"
334  inverse(1,1) = d2f(2,2) / determ
335  inverse(2,2) = d2f(1,1) / determ
336  inverse(1,2) = -d2f(1,2) / determ
337  inverse(2,1) = -d2f(2,1) / determ
338  dr=matmul(inverse,df)
339 
340  tol=dot_product(dr,dr)
341  if( dsqrt(tol)> 3.d0 ) then ! too far away
342  r= -100.d0; exit
343  endif
344 
345  factor = 1.d0
346  do order=1,10
347  r_tmp(1:2) = r(1:2) + factor*dr(1:2)
348  call getshapefunc( etype, r_tmp, sfunc )
349  xyz_out(1:3) = matmul( elemt(1:3,1:nn), sfunc(:) )
350  dxyz(1:3) = xyz(1:3)-xyz_out(:)
351  dist_now = dot_product( dxyz, dxyz )
352  if(dist_now <= dist_last) exit
353  factor = factor*0.7d0
354  enddo
355  r(1:2) = r_tmp(1:2)
356 
357  if( tol<eps ) exit
358  enddo
359 
360  isin = .false.
361  cstate%state = contactfree
362  if( isinsideelement( etype, r, clr )>=0 ) then
363  dxyz(:)=xyz_out(:)-xyz(:)
364  normal(:) = surfacenormal( etype, nn, r, elemt(1:3,1:nn) )
365  normal(:) = normal(:)/dsqrt( dot_product(normal, normal) )
366  do count = 1,3
367  if( dabs(normal(count))<1.d-10 ) normal(count) =0.d0
368  if( dabs(1.d0-dabs(normal(count)))<1.d-10 ) normal(count) =sign(1.d0, normal(count))
369  enddo
370  cstate%distance = dot_product( dxyz, normal )
371 
372  if( cstate%interference_flag == c_if_slave) then
373  if( cstate%init_pos == 0.d0 .and. cstate%distance < cstate%end_pos) then
374  cstate%init_pos = cstate%distance
375  cstate%time_factor = (cstate%end_pos - cstate%distance) / cstate%time_factor
376  cstate%shrink_factor = cstate%distance
377  end if
378  end if
379 
380  if( cstate%interference_flag == 0)then ! not shrink-node
381  if( cstate%distance < distclr*reflen .and. cstate%distance > -5.0d-01*reflen ) isin = .true.
382  else
383  if( cstate%distance < cstate%shrink_factor + distclr*reflen ) isin = .true.
384  end if
385 
386  if( isin ) then
387  if( initstate== contactfree ) then
388  cstate%state = contactstick
389  else
390  cstate%state = initstate
391  endif
392  cstate%gpos(:)=xyz_out(:)
393  cstate%lpos(1:2)=r(:)
394  cstate%direction(:) = normal(:)
395  cstate%wkdist = cstate%distance
396  endif
397  endif
398  end subroutine project_point2element
399 
401  subroutine project_point2solidelement(xyz,etype,nn,elemt,reflen,cstate,isin,distclr,ctpos,localclr)
402  use m_utilities
403  real(kind=kreal),intent(in) :: xyz(3)
404  integer, intent(in) :: etype
405  integer, intent(in) :: nn
406  real(kind=kreal),intent(in) :: elemt(:,:)
407  real(kind=kreal),intent(in) :: reflen
408  type(tcontactstate),intent(inout) :: cstate
409  logical, intent(out) :: isin
410  real(kind=kreal), intent(in) :: distclr
411  real(kind=kreal), optional :: ctpos(3)
412  real(kind=kreal), optional :: localclr
413 
414  integer :: count,order, initstate
415  real(kind=kreal) :: determ, inverse(3,3)
416  real(kind=kreal) :: sfunc(nn), deriv(nn,3)
417  real(kind=kreal) :: r(3), dr(3), r_tmp(3) ! natural coordinate
418  real(kind=kreal) :: xyz_out(3) ! curr. projection position
419  real(kind=kreal) :: dist_last,dist_now, dxyz(3) ! dist between the point and its projection
420  real(kind=kreal) :: tangent(3,2) ! base vectors in tangent space
421  real(kind=kreal) :: df(2),d2f(2,2),normal(3)
422  real(kind=kreal),parameter :: eps = 1.0d-8
423  real(kind=kreal) :: clr, tol, factor
424 
425  initstate = cstate%state
426  clr = 1.d-4
427  if( present( localclr ) ) clr=localclr
428  if( present( ctpos ) ) then
429  r(:)= ctpos
430  else
431  call getelementcenter( etype, r(:) )
432  endif
433 
434  tol = 1.0d0
435  do count=1,100
436  call getshapefunc( etype, r, sfunc )
437  xyz_out = matmul( elemt(1:3,1:nn), sfunc(1:nn) )
438  dxyz(1:3) = xyz_out(1:3) - xyz(1:3)
439  dist_last = dot_product( dxyz, dxyz(:) )
440 
441  call getshapederiv( etype, r, deriv )
442  inverse(1:3,1:3) = matmul(elemt(1:3,1:nn),deriv(1:nn,1:3))
443  call calinverse(3, inverse(1:3,1:3))
444  dr(1:3) = -matmul(inverse(1:3,1:3),dxyz(1:3))
445 
446  tol=dot_product(dr,dr)
447  if( count > 1 .and. dsqrt(tol)> 3.d0 ) then ! too far away
448  r= -100.d0; exit
449  endif
450 
451  factor = 1.d0
452  do order=1,10
453  r_tmp(1:3) = r(1:3) + factor*dr(1:3)
454  call getshapefunc( etype, r_tmp, sfunc )
455  xyz_out(1:3) = matmul( elemt(1:3,1:nn), sfunc(1:nn) )
456  dxyz(1:3) = xyz(1:3)-xyz_out(1:3)
457  dist_now = dot_product( dxyz, dxyz )
458  if(dist_now <= dist_last) exit
459  factor = factor*0.7d0
460  enddo
461  r(1:3) = r_tmp(1:3)
462 
463  if( tol<eps ) exit
464  enddo
465 
466  isin = .false.
467  cstate%state = contactfree
468  if( isinside3delement( etype, r, clr )>=0 ) then
469  dxyz(:)=xyz_out(:)-xyz(:)
470  cstate%distance = dsqrt(dot_product( dxyz, dxyz ))
471 
472  if( cstate%distance < distclr ) isin = .true.
473 
474  if( isin ) then
475  if( initstate== contactfree ) then
476  cstate%state = contactstick
477  else
478  cstate%state = initstate
479  endif
480  cstate%gpos(:)=xyz_out(:)
481  cstate%direction(:) = (/1.d0,0.d0,0.d0/)
482  cstate%lpos(1:3)=r(:)
483  endif
484  endif
485  end subroutine
486 
488  subroutine update_tangentforce(etype,nn,elemt0,elemt,cstate)
489  integer, intent(in) :: etype
490  integer, intent(in) :: nn
491  real(kind=kreal),intent(in) :: elemt0(3,nn)
492  real(kind=kreal),intent(in) :: elemt(3,nn)
493  type(tcontactstate), intent(inout) :: cstate
494 
495  integer :: i
496  real(kind=kreal) :: tangent0(3,2), tangent(3,2) ! base vectors in tangent space
497  real(kind=kreal) :: coeff(2), norm, norm_tmp
498  real(kind=kreal) :: tangentforce_tmp(3)
499 
500  call tangentbase( etype, nn, cstate%lpos(1:2), elemt0, tangent0 )
501  call tangentbase( etype, nn, cstate%lpos(1:2), elemt, tangent )
502 
503  !project tangentforce to base vector tangent0
504  do i=1,2
505  coeff(i) = dot_product(cstate%tangentForce(1:3),tangent0(1:3,i))
506  coeff(i) = coeff(i)/dot_product(tangent0(1:3,i),tangent0(1:3,i))
507  enddo
508  tangentforce_tmp(1:3) = coeff(1)*tangent0(1:3,1) + coeff(2)*tangent0(1:3,2)
509  norm_tmp = dsqrt(dot_product(tangentforce_tmp,tangentforce_tmp))
510  !adjust tangent force of slave point which moved over element boundary
511  if( norm_tmp > 1.d-6 ) then
512  norm = dsqrt(dot_product(cstate%tangentForce,cstate%tangentForce))
513  coeff(1:2) = (norm/norm_tmp)*coeff(1:2)
514  end if
515 
516  !set rotated tangentforce to tangentforce1
517  cstate%tangentForce1(1:3) = coeff(1)*tangent(1:3,1) + coeff(2)*tangent(1:3,2)
518 
519  end subroutine update_tangentforce
520 
521  subroutine set_shrink_factor(ctime, cstate, etime, if_type)
522  real(kind=kreal),intent(in) :: ctime, etime
523  type(tcontactstate), intent(inout) :: cstate
524  integer, intent(in) :: if_type
525 
526  if (if_type == c_if_slave .and. cstate%init_pos == 0.d0) then
527  cstate%shrink_factor = 0.0d0; return
528  end if
529  cstate%shrink_factor = cstate%time_factor*ctime + cstate%init_pos
530  if(ctime >= etime) cstate%shrink_factor = cstate%end_pos
531 
532  end subroutine set_shrink_factor
533 
534 end module m_contact_lib
535 
536 
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:647
elementinfo::getnumberofnodes
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:131
elementinfo::isinside3delement
integer function isinside3delement(fetype, localcoord, clearance)
if a point is inside a surface element -1: No; 0: Yes; >0: Node's (vertex) number
Definition: element.f90:1101
elementinfo::getelementcenter
subroutine getelementcenter(fetype, localcoord)
Return natural coordinate of the center of element.
Definition: element.f90:1012
elementinfo::getshapederiv
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
Definition: element.f90:578
elementinfo::tangentbase
subroutine tangentbase(fetype, nn, localcoord, elecoord, tangent)
Calculate base vector of tangent space of 3d surface.
Definition: element.f90:933
elementinfo::isinsideelement
integer function isinsideelement(fetype, localcoord, clearance)
if a point is inside a surface element -1: No; 0: Yes; >0: Node's (vertex) number
Definition: element.f90:1036
elementinfo::surfacenormal
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:870
elementinfo::curvature
subroutine curvature(fetype, nn, localcoord, elecoord, l2ndderiv, normal, curv)
Calculate curvature tensor at a point along 3d surface.
Definition: element.f90:967
m_contact_lib
This module provide functions of contact stiffness calculation.
Definition: contact_lib.f90:6
m_contact_lib::contacttied
integer, parameter contacttied
contact type or algorithm definition
Definition: contact_lib.f90:21
m_contact_lib::tied2stiff
subroutine tied2stiff(flag, cstate, etype, nnode, mu, mut, stiff, force)
This subroutine calculate contact stiff matrix and contact force.
Definition: contact_lib.f90:191
m_contact_lib::contactfslid
integer, parameter contactfslid
Definition: contact_lib.f90:24
m_contact_lib::print_contact_state
subroutine print_contact_state(fnum, cstate)
Print out contact state.
Definition: contact_lib.f90:87
m_contact_lib::getmetrictensor
subroutine getmetrictensor(pos, etype, ele, tensor)
This subroutine calculate the metric tensor of a elemental surface.
Definition: contact_lib.f90:223
m_contact_lib::contactglued
integer, parameter contactglued
Definition: contact_lib.f90:22
m_contact_lib::contactunknown
integer, parameter contactunknown
Definition: contact_lib.f90:14
m_contact_lib::contact2mpcval
subroutine contact2mpcval(cstate, etype, nnode, mpcval)
Transfer contact condition int mpc bundary conditions.
Definition: contact_lib.f90:98
m_contact_lib::project_point2solidelement
subroutine project_point2solidelement(xyz, etype, nn, elemt, reflen, cstate, isin, distclr, ctpos, localclr)
This subroutine find the projection of a slave point onto master surface.
Definition: contact_lib.f90:402
m_contact_lib::project_point2element
subroutine project_point2element(xyz, etype, nn, elemt, reflen, cstate, isin, distclr, ctpos, localclr)
This subroutine find the projection of a slave point onto master surface.
Definition: contact_lib.f90:282
m_contact_lib::update_tangentforce
subroutine update_tangentforce(etype, nn, elemt0, elemt, cstate)
This subroutine find the projection of a slave point onto master surface.
Definition: contact_lib.f90:489
m_contact_lib::contact_state_init
subroutine contact_state_init(cstate)
Initializer.
Definition: contact_lib.f90:57
m_contact_lib::contact2stiff
subroutine contact2stiff(flag, cstate, etype, nnode, ele, mu, mut, fcoeff, symm, stiff, force)
This subroutine calculate contact stiff matrix and contact force.
Definition: contact_lib.f90:119
m_contact_lib::set_shrink_factor
subroutine set_shrink_factor(ctime, cstate, etime, if_type)
Definition: contact_lib.f90:522
m_contact_lib::contactslip
integer, parameter contactslip
Definition: contact_lib.f90:18
m_contact_lib::contactsslid
integer, parameter contactsslid
Definition: contact_lib.f90:23
m_contact_lib::contactfree
integer, parameter contactfree
contact state definition
Definition: contact_lib.f90:16
m_contact_lib::c_if_slave
integer, parameter c_if_slave
contact interference type
Definition: contact_lib.f90:27
m_contact_lib::dispincrematrix
subroutine dispincrematrix(pos, etype, nnode, ele, tangent, tensor, matrix)
This subroutine calculate the relation between global disp and displacement along natural coordinate ...
Definition: contact_lib.f90:241
m_contact_lib::contact_state_copy
subroutine contact_state_copy(cstate1, cstate2)
Copy.
Definition: contact_lib.f90:80
m_contact_lib::contactstick
integer, parameter contactstick
Definition: contact_lib.f90:17
m_contact_lib::c_if_master
integer, parameter c_if_master
Definition: contact_lib.f90:28
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_utilities::calinverse
subroutine calinverse(NN, A)
calculate inverse of matrix a
Definition: utilities.f90:295
m_contact_lib::tcontactstate
This structure records contact status.
Definition: contact_lib.f90:30