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