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