FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
fstr_contact_geom.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_fstr_contact_geom
7  use hecmw
8  use elementinfo
9  use mcontactdef
10  use m_fstr_contact_smoothing
11  implicit none
12 
13  public :: cal_node_normal
14 
15 contains
16 
17  ! TODO: Move from fstr_contact_lib.f90
19  subroutine project_point2element(xyz,etype,nn,elemt,reflen,cstate,isin,distclr,ctpos,localclr)
20  real(kind=kreal),intent(in) :: xyz(3)
21  integer, intent(in) :: etype
22  integer, intent(in) :: nn
23  real(kind=kreal),intent(in) :: elemt(:,:)
24  real(kind=kreal),intent(in) :: reflen
25  type(tcontactstate),intent(inout) :: cstate
26  logical, intent(out) :: isin
27  real(kind=kreal), intent(in) :: distclr
28  real(kind=kreal), optional :: ctpos(2)
29  real(kind=kreal), optional :: localclr
30 
31  integer :: count,order
32  real(kind=kreal) :: determ, inverse(2,2)
33  real(kind=kreal) :: sfunc(nn), curv(3,2,2)
34  real(kind=kreal) :: r(2), dr(2), r_tmp(2) ! natural coordinate
35  real(kind=kreal) :: xyz_out(3) ! curr. projection position
36  real(kind=kreal) :: dist_last,dist_now, dxyz(3) ! dist between the point and its projection
37  real(kind=kreal) :: tangent(3,2) ! base vectors in tangent space
38  real(kind=kreal) :: df(2),d2f(2,2),normal(3)
39  real(kind=kreal),parameter :: eps = 1.0d-8
40  real(kind=kreal) :: clr, tol, factor
41 
42  clr = 1.d-4
43  if( present( localclr ) ) clr=localclr
44  if( present( ctpos ) ) then
45  r(:)= ctpos
46  else
47  call getelementcenter( etype, r(:) )
48  endif
49 
50  tol = 1.0d0
51  do count=1,100
52  call getshapefunc( etype, r, sfunc )
53  xyz_out = matmul( elemt(1:3,1:nn), sfunc )
54  dxyz(1:3) = xyz_out(1:3) - xyz(1:3)
55  dist_last = dot_product( dxyz, dxyz(:) )
56 
57  call tangentbase( etype, nn, r, elemt(1:3,1:nn), tangent )
58  call curvature( etype, nn, r, elemt(1:3,1:nn), curv )
59 
60  ! dF(1:2)
61  df(1:2) = -matmul( dxyz(:), tangent(:,:) )
62  ! d2F(1:2,1:2)
63  d2f(1,1)= dot_product( tangent(:,1), tangent(:,1) ) - dot_product( dxyz, curv(:,1,1) )
64  d2f(1,2)= dot_product( tangent(:,1), tangent(:,2) ) - dot_product( dxyz, curv(:,1,2) )
65  d2f(2,1)= dot_product( tangent(:,2), tangent(:,1) ) - dot_product( dxyz, curv(:,2,1) )
66  d2f(2,2)= dot_product( tangent(:,2), tangent(:,2) ) - dot_product( dxyz, curv(:,2,2) )
67 
68  ! inverse of d2F
69  determ = d2f(1,1)*d2f(2,2) - d2f(1,2)*d2f(2,1)
70  if( determ==0.d0 ) stop "Math error in contact searching"
71  inverse(1,1) = d2f(2,2) / determ
72  inverse(2,2) = d2f(1,1) / determ
73  inverse(1,2) = -d2f(1,2) / determ
74  inverse(2,1) = -d2f(2,1) / determ
75  dr=matmul(inverse,df)
76 
77  tol=dot_product(dr,dr)
78  if( dsqrt(tol)> 3.d0 ) then ! too far away
79  r= -100.d0; exit
80  endif
81 
82  factor = 1.d0
83  do order=1,10
84  r_tmp(1:2) = r(1:2) + factor*dr(1:2)
85  call getshapefunc( etype, r_tmp, sfunc )
86  xyz_out(1:3) = matmul( elemt(1:3,1:nn), sfunc(:) )
87  dxyz(1:3) = xyz(1:3)-xyz_out(:)
88  dist_now = dot_product( dxyz, dxyz )
89  if(dist_now <= dist_last) exit
90  factor = factor*0.7d0
91  enddo
92  r(1:2) = r_tmp(1:2)
93 
94  if( tol<eps ) exit
95  enddo
96 
97  ! The contact state itself is left untouched until the candidate is
98  ! accepted, so that a failed attempt does not destroy the stick or slip
99  ! state that a subsequent candidate has to inherit.
100  isin = .false.
101  if( isinsideelement( etype, r, clr )>=0 ) then
102  dxyz(:)=xyz_out(:)-xyz(:)
103  normal(:) = surfacenormal( etype, nn, r, elemt(1:3,1:nn) )
104  normal(:) = normal(:)/dsqrt( dot_product(normal, normal) )
105  do count = 1,3
106  if( dabs(normal(count))<1.d-10 ) normal(count) =0.d0
107  if( dabs(1.d0-dabs(normal(count)))<1.d-10 ) normal(count) =sign(1.d0, normal(count))
108  enddo
109  cstate%distance = dot_product( dxyz, normal )
110 
111  if( cstate%interference_flag == c_if_slave) then
112  if( cstate%init_pos == 0.d0 .and. cstate%distance < cstate%end_pos) then
113  cstate%init_pos = cstate%distance
114  cstate%time_factor = (cstate%end_pos - cstate%distance) / cstate%time_factor
115  cstate%shrink_factor = cstate%distance
116  end if
117  end if
118 
119  if( cstate%interference_flag == 0)then ! not shrink-node
120  if( cstate%distance < distclr*reflen .and. cstate%distance > -5.0d-01*reflen ) isin = .true.
121  else
122  if( cstate%distance < cstate%shrink_factor + distclr*reflen ) isin = .true.
123  end if
124 
125  if( isin ) then
126  if( cstate%state == contactfree ) cstate%state = contactstick
127  cstate%gpos(:)=xyz_out(:)
128  cstate%lpos(1:2)=r(:)
129  cstate%direction(:) = normal(:)
130  cstate%wkdist = cstate%distance
131  endif
132  endif
133  end subroutine project_point2element
134 
136  subroutine project_point2solidelement(xyz,etype,nn,elemt,reflen,cstate,isin,distclr,ctpos,localclr)
137  use m_utilities
138  real(kind=kreal),intent(in) :: xyz(3)
139  integer, intent(in) :: etype
140  integer, intent(in) :: nn
141  real(kind=kreal),intent(in) :: elemt(:,:)
142  real(kind=kreal),intent(in) :: reflen
143  type(tcontactstate),intent(inout) :: cstate
144  logical, intent(out) :: isin
145  real(kind=kreal), intent(in) :: distclr
146  real(kind=kreal), optional :: ctpos(3)
147  real(kind=kreal), optional :: localclr
148 
149  integer :: count,order
150  real(kind=kreal) :: inverse(3,3)
151  real(kind=kreal) :: sfunc(nn), deriv(nn,3)
152  real(kind=kreal) :: r(3), dr(3), r_tmp(3) ! natural coordinate
153  real(kind=kreal) :: xyz_out(3) ! curr. projection position
154  real(kind=kreal) :: dist_last,dist_now, dxyz(3) ! dist between the point and its projection
155  real(kind=kreal),parameter :: eps = 1.0d-8
156  real(kind=kreal) :: clr, tol, factor
157 
158  clr = 1.d-4
159  if( present( localclr ) ) clr=localclr
160  if( present( ctpos ) ) then
161  r(:)= ctpos
162  else
163  call getelementcenter( etype, r(:) )
164  endif
165 
166  tol = 1.0d0
167  do count=1,100
168  call getshapefunc( etype, r, sfunc )
169  xyz_out = matmul( elemt(1:3,1:nn), sfunc(1:nn) )
170  dxyz(1:3) = xyz_out(1:3) - xyz(1:3)
171  dist_last = dot_product( dxyz, dxyz(:) )
172 
173  call getshapederiv( etype, r, deriv )
174  inverse(1:3,1:3) = matmul(elemt(1:3,1:nn),deriv(1:nn,1:3))
175  call calinverse(3, inverse(1:3,1:3))
176  dr(1:3) = -matmul(inverse(1:3,1:3),dxyz(1:3))
177 
178  tol=dot_product(dr,dr)
179  if( count > 1 .and. dsqrt(tol)> 3.d0 ) then ! too far away
180  r= -100.d0; exit
181  endif
182 
183  factor = 1.d0
184  do order=1,10
185  r_tmp(1:3) = r(1:3) + factor*dr(1:3)
186  call getshapefunc( etype, r_tmp, sfunc )
187  xyz_out(1:3) = matmul( elemt(1:3,1:nn), sfunc(1:nn) )
188  dxyz(1:3) = xyz(1:3)-xyz_out(1:3)
189  dist_now = dot_product( dxyz, dxyz )
190  if(dist_now <= dist_last) exit
191  factor = factor*0.7d0
192  enddo
193  r(1:3) = r_tmp(1:3)
194 
195  if( tol<eps ) exit
196  enddo
197 
198  isin = .false.
199  if( isinside3delement( etype, r, clr )>=0 ) then
200  dxyz(:)=xyz_out(:)-xyz(:)
201  cstate%distance = dsqrt(dot_product( dxyz, dxyz ))
202 
203  if( cstate%distance < distclr ) isin = .true.
204 
205  if( isin ) then
206  if( cstate%state == contactfree ) cstate%state = contactstick
207  cstate%gpos(:)=xyz_out(:)
208  cstate%direction(:) = (/1.d0,0.d0,0.d0/)
209  cstate%lpos(1:3)=r(:)
210  endif
211  endif
212  end subroutine
213 
216  subroutine project_point2surfelement(xyz, surf, currpos, cstate, isin, distclr, ctpos, localclr, smoothing)
217  real(kind=kreal), intent(in) :: xyz(3)
218  type(tsurfelement), intent(in) :: surf
219  real(kind=kreal), intent(in) :: currpos(:)
220  type(tcontactstate), intent(inout) :: cstate
221  logical, intent(out) :: isin
222  real(kind=kreal), intent(in) :: distclr
223  real(kind=kreal), optional, intent(in) :: ctpos(2)
224  real(kind=kreal), optional, intent(in) :: localclr
225  integer(kind=kint), intent(in) :: smoothing
226 
227  integer(kind=kint) :: nn, j, iss, etype_use, nn_use
228  real(kind=kreal) :: elem(3, l_max_elem_node)
229  real(kind=kreal) :: elem_tri3n(3, 3)
230 
231  ! Extract element coordinates from surf structure
232  nn = size(surf%nodes)
233  do j = 1, nn
234  iss = surf%nodes(j)
235  elem(1:3, j) = currpos(3*iss-2:3*iss)
236  enddo
237 
238  if (smoothing == kcsnagata) then
239  if (surf%etype == fe_tri3n) then
240  elem_tri3n = elem(1:3, 1:3)
241  call reorder_tri3n_to_tri6n(elem_tri3n, surf%intermediate_points, elem(1:3,1:6))
242  etype_use = fe_tri6n
243  nn_use = 6
244  else if (surf%etype == fe_quad4n) then
245  do j = 1, nn
246  elem(1:3, nn+j) = surf%intermediate_points(1:3, j)
247  enddo
248  etype_use = fe_quad8n
249  nn_use = 8
250  else
251  etype_use = surf%etype
252  nn_use = nn
253  endif
254  else
255  etype_use = surf%etype
256  nn_use = nn
257  endif
258 
259  call project_point2element(xyz, etype_use, nn_use, elem, surf%reflen, cstate, &
260  isin, distclr, ctpos, localclr)
261 
262  end subroutine project_point2surfelement
263 
265  subroutine getmetrictensor( pos, etype, ele, tensor )
266  real(kind=kreal), intent(in) :: pos(2)
267  integer, intent(in) :: etype
268  real(kind=kreal), intent(in) :: ele(:,:)
269  real(kind=kreal), intent(out) :: tensor(2,2)
270 
271  integer :: nn
272  real(kind=kreal) :: tangent(3,2)
273  nn= getnumberofnodes(etype)
274  call tangentbase( etype, nn, pos, ele, tangent )
275  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
276  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
277  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
278  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
279  end subroutine
280 
283  subroutine dispincrematrix( pos, etype, nnode, ele, tangent, tensor, matrix )
284  real(kind=kreal), intent(in) :: pos(2)
285  integer, intent(in) :: etype
286  integer, intent(in) :: nnode
287  real(kind=kreal), intent(in) :: ele(3,nnode)
288  real(kind=kreal), intent(out) :: tangent(3,2)
289  real(kind=kreal), intent(out) :: tensor(2,2)
290  real(kind=kreal), intent(out) :: matrix(2,nnode*3+3)
291 
292  integer :: i,j
293  real(kind=kreal) :: det
294  real(kind=kreal) :: shapefunc(nnode), t1(nnode*3+3), t2(nnode*3+3)
295  call tangentbase( etype, nnode, pos, ele, tangent )
296  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
297  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
298  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
299  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
300  det = tensor(1,1)*tensor(2,2)-tensor(1,2)*tensor(2,1)
301  if( det==0.d0 ) stop "Error in calculate DispIncreMatrix"
302  ! inverse(1,1) = tensor(2,2)/det
303  ! inverse(1,2) = -tensor(1,2)/det
304  ! inverse(2,1) = -tensor(2,1)/det
305  ! inverse(2,2) = tensor(1,1)/det
306 
307  call getshapefunc( etype, pos(:), shapefunc )
308  forall( j=1:3 )
309  t1( j ) = tangent(j,1)
310  t2( j ) = tangent(j,2)
311  end forall
312  forall( i=1:nnode, j=1:3 )
313  t1( i*3+j ) = -tangent(j,1)*shapefunc(i)
314  t2( i*3+j ) = -tangent(j,2)*shapefunc(i)
315  end forall
316  !matrix( 1:2,: ) = matmul( inverse(:,:), matrix )
317  matrix(1,:) = (tensor(2,2)*t1(:)-tensor(1,2)*t2(:))/det
318  matrix(2,:) = (tensor(1,1)*t2(:)-tensor(2,1)*t1(:))/det
319  tangent(:,1) = tangent(:,1)/dsqrt(dot_product(tangent(:,1),tangent(:,1)))
320  tangent(:,2) = tangent(:,2)/dsqrt(dot_product(tangent(:,2),tangent(:,2)))
321  end subroutine
322 
324  subroutine update_direction( nslave, contact, currpos )
325  integer, intent(in) :: nslave
326  type( tcontact ), intent(inout) :: contact
327  real(kind=kreal), intent(in) :: currpos(:)
328 
329  integer(kind=kint) :: slave, sid0, etype, iSS
330  real(kind=kreal) :: coord(3)
331  logical :: isin
332  type(tcontactstate) :: cstate_tmp
333 
334  slave = contact%slave(nslave)
335  coord(:) = currpos(3*slave-2:3*slave)
337  sid0 = contact%states(nslave)%surface
338 
339  cstate_tmp = contact%states(nslave)
340  call project_point2surfelement( coord, contact%master(sid0), currpos, &
341  cstate_tmp, isin, contact%cparam%DISTCLR_NOCHECK, &
342  contact%states(nslave)%lpos, contact%cparam%CLR_SAME_ELEM, &
343  smoothing=contact%smoothing )
344 
345  if( isin ) then
346  etype = contact%master(sid0)%etype
347  iss = isinsideelement( etype, cstate_tmp%lpos, contact%cparam%CLR_CAL_NORM )
348  if( iss>0 .and. contact%smoothing /= kcsnagata ) &
349  call cal_node_normal( cstate_tmp%surface, iss, contact%master, currpos, &
350  cstate_tmp%lpos, cstate_tmp%direction(:) )
351  endif
352 
353  contact%states(nslave)%direction = cstate_tmp%direction
354 
355  end subroutine
356 
358  subroutine cal_node_normal( csurf, isin, surf, currpos, lpos, normal )
359  use elementinfo, only:getvertexcoord, surfacenormal
360  integer, intent(in) :: csurf
361  integer, intent(in) :: isin
362  type(tsurfelement), intent(in) :: surf(:)
363  real(kind=kreal), intent(in) :: currpos(:)
364  real(kind=kreal), intent(in) :: lpos(:)
365  real(kind=kreal), intent(out) :: normal(3)
366  integer(kind=kint) :: cnode, i, j, cnt, nd1, gn, etype, iss, nn,cgn
367  real(kind=kreal) :: cnpos(2), elem(3, l_max_elem_node)
368  integer(kind=kint) :: cnode1, cnode2, gn1, gn2, nsurf, cgn1, cgn2, isin_n
369  real(kind=kreal) :: x, normal_n(3), lpos_n(2)
370 
371  if( 1 <= isin .and. isin <= 4 ) then ! corner
372  cnode = isin
373  gn = surf(csurf)%nodes(cnode)
374  etype = surf(csurf)%etype
375  call getvertexcoord( etype, cnode, cnpos )
376  nn = size( surf(csurf)%nodes )
377  do j=1,nn
378  iss = surf(csurf)%nodes(j)
379  elem(1:3,j)=currpos(3*iss-2:3*iss)
380  enddo
381  normal = surfacenormal( etype, nn, cnpos, elem )
382  cnt = 1
383  do i=1,surf(csurf)%n_neighbor
384  nd1 = surf(csurf)%neighbor(i)
385  nn = size( surf(nd1)%nodes )
386  etype = surf(nd1)%etype
387  cgn = 0
388  do j=1,nn
389  iss = surf(nd1)%nodes(j)
390  elem(1:3,j)=currpos(3*iss-2:3*iss)
391  if( iss==gn ) cgn=j
392  enddo
393  if( cgn>0 ) then
394  call getvertexcoord( etype, cgn, cnpos )
395  !normal = normal+SurfaceNormal( etype, nn, cnpos, elem )
396  normal_n = surfacenormal( etype, nn, cnpos, elem )
397  normal = normal+normal_n
398  cnt = cnt+1
399  endif
400  enddo
401  !normal = normal/cnt !!-???
402  elseif( 12 <= isin .and. isin <= 41 ) then ! edge
403  cnode1 = isin / 10
404  cnode2 = mod(isin, 10)
405  gn1 = surf(csurf)%nodes(cnode1)
406  gn2 = surf(csurf)%nodes(cnode2)
407  etype = surf(csurf)%etype
408  nn = size( surf(csurf)%nodes )
409  do j=1,nn
410  iss = surf(csurf)%nodes(j)
411  elem(1:3,j)=currpos(3*iss-2:3*iss)
412  enddo
413  normal = surfacenormal( etype, nn, lpos, elem )
414  select case (etype)
415  case (fe_tri3n, fe_tri6n, fe_tri6nc)
416  if ( isin==12 ) then; x=lpos(2)-lpos(1)
417  elseif( isin==23 ) then; x=1.d0-2.d0*lpos(2)
418  elseif( isin==31 ) then; x=2.d0*lpos(1)-1.d0
419  else; stop "Error: cal_node_normal: invalid isin"
420  endif
421  case (fe_quad4n, fe_quad8n)
422  if ( isin==12 ) then; x=lpos(1)
423  elseif( isin==23 ) then; x=lpos(2)
424  elseif( isin==34 ) then; x=-lpos(1)
425  elseif( isin==41 ) then; x=-lpos(2)
426  else; stop "Error: cal_node_normal: invalid isin"
427  endif
428  end select
429  ! find neighbor surf that includes cnode1 and cnode2
430  nsurf = 0
431  neib_loop: do i=1, surf(csurf)%n_neighbor
432  nd1 = surf(csurf)%neighbor(i)
433  nn = size( surf(nd1)%nodes )
434  etype = surf(nd1)%etype
435  cgn1 = 0
436  cgn2 = 0
437  do j=1,nn
438  iss = surf(nd1)%nodes(j)
439  elem(1:3,j)=currpos(3*iss-2:3*iss)
440  if( iss==gn1 ) cgn1=j
441  if( iss==gn2 ) cgn2=j
442  enddo
443  if( cgn1>0 .and. cgn2>0 ) then
444  nsurf = nd1
445  isin_n = 10*cgn2 + cgn1
446  x = -x
447  select case (etype)
448  case (fe_tri3n, fe_tri6n, fe_tri6nc)
449  if ( isin_n==12 ) then; lpos_n(1)=0.5d0*(1.d0-x); lpos_n(2)=0.5d0*(1.d0+x)
450  elseif( isin_n==23 ) then; lpos_n(1)=0.d0; lpos_n(2)=0.5d0*(1.d0-x)
451  elseif( isin_n==31 ) then; lpos_n(1)=0.5d0*(1.d0+x); lpos_n(2)=0.d0
452  else; stop "Error: cal_node_normal: invalid isin_n"
453  endif
454  case (fe_quad4n, fe_quad8n)
455  if ( isin_n==12 ) then; lpos_n(1)= x; lpos_n(2)=-1.d0
456  elseif( isin_n==23 ) then; lpos_n(1)= 1.d0; lpos_n(2)= x
457  elseif( isin_n==34 ) then; lpos_n(1)=-x; lpos_n(2)= 1.d0
458  elseif( isin_n==41 ) then; lpos_n(1)=-1.d0; lpos_n(2)=-x
459  else; stop "Error: cal_node_normal: invalid isin_n"
460  endif
461  end select
462  !normal = normal + SurfaceNormal( etype, nn, lpos_n, elem )
463  normal_n = surfacenormal( etype, nn, lpos_n, elem )
464  normal = normal+normal_n
465  exit neib_loop
466  endif
467  enddo neib_loop
468  !if( nsurf==0 ) write(0,*) "Warning: cal_node_normal: neighbor surf not found"
469  !normal = normal/2
470  endif
471  normal = normal/ dsqrt( dot_product( normal, normal ) )
472  end subroutine cal_node_normal
473 
474 end module m_fstr_contact_geom
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:696
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:1150
elementinfo::getelementcenter
subroutine getelementcenter(fetype, localcoord)
Return natural coordinate of the center of element.
Definition: element.f90:1061
elementinfo::fe_tri6n
integer, parameter fe_tri6n
Definition: element.f90:70
elementinfo::getshapederiv
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
Definition: element.f90:627
elementinfo::getvertexcoord
subroutine getvertexcoord(fetype, cnode, localcoord)
Get the natural coord of a vertex node.
Definition: element.f90:1192
elementinfo::fe_tri3n
integer, parameter fe_tri3n
Definition: element.f90:69
elementinfo::tangentbase
subroutine tangentbase(fetype, nn, localcoord, elecoord, tangent)
Calculate base vector of tangent space of 3d surface.
Definition: element.f90:982
elementinfo::fe_quad4n
integer, parameter fe_quad4n
Definition: element.f90:72
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:1085
elementinfo::surfacenormal
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:919
elementinfo::curvature
subroutine curvature(fetype, nn, localcoord, elecoord, l2ndderiv, normal, curv)
Calculate curvature tensor at a point along 3d surface.
Definition: element.f90:1016
elementinfo::fe_quad8n
integer, parameter fe_quad8n
Definition: element.f90:73
hecmw
Definition: hecmw.f90:6
m_fstr_contact_geom
This module provides geometric calculations for contact.
Definition: fstr_contact_geom.f90:6
m_fstr_contact_geom::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: fstr_contact_geom.f90:137
m_fstr_contact_geom::project_point2surfelement
subroutine project_point2surfelement(xyz, surf, currpos, cstate, isin, distclr, ctpos, localclr, smoothing)
Wrapper for project_Point2Element that takes tSurfElement structure This subroutine handles element c...
Definition: fstr_contact_geom.f90:217
m_fstr_contact_geom::update_direction
subroutine update_direction(nslave, contact, currpos)
This subroutine tracks down next contact position after a finite slide.
Definition: fstr_contact_geom.f90:325
m_fstr_contact_geom::cal_node_normal
subroutine, public cal_node_normal(csurf, isin, surf, currpos, lpos, normal)
Calculate averaged nodal normal.
Definition: fstr_contact_geom.f90:359
m_fstr_contact_geom::dispincrematrix
subroutine dispincrematrix(pos, etype, nnode, ele, tangent, tensor, matrix)
This subroutine calculate the relation between global disp and displacement along natural coordinate ...
Definition: fstr_contact_geom.f90:284
m_fstr_contact_geom::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: fstr_contact_geom.f90:20
m_fstr_contact_geom::getmetrictensor
subroutine getmetrictensor(pos, etype, ele, tensor)
This subroutine calculate the metric tensor of a elemental surface.
Definition: fstr_contact_geom.f90:266
m_fstr_contact_smoothing
Contact surface smoothing using Nagata patch interpolation.
Definition: fstr_contact_smoothing.f90:6
m_fstr_contact_smoothing::reorder_tri3n_to_tri6n
subroutine, public reorder_tri3n_to_tri6n(vertices_in, midpoints_in, nodes_out)
Reorder triangle vertices and midpoints from fe_tri3n to fe_tri6n order fe_tri3n: v1(xi),...
Definition: fstr_contact_smoothing.f90:32
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
mcontactdef
This module manages the data structure for contact calculation.
Definition: fstr_contact_def.F90:12
mcontactdef::c_if_slave
integer, parameter c_if_slave
contact interference type
Definition: fstr_contact_def.F90:50
mcontactdef::cgn
real(kind=kreal), save cgn
convergent condition of penetration
Definition: fstr_contact_def.F90:24
mcontactdef::contactfree
integer, parameter contactfree
contact state definition
Definition: fstr_contact_def.F90:34
mcontactdef::contactstick
integer, parameter contactstick
Definition: fstr_contact_def.F90:36
mcontactdef::kcsnagata
integer, parameter kcsnagata
Definition: fstr_contact_def.F90:47
mcontactdef::tcontact
Structure to includes all info needed by contact calculation.
Definition: fstr_contact_def.F90:82
mcontactdef::tcontactstate
This structure records contact status.
Definition: fstr_contact_def.F90:58