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, initstate
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  initstate = cstate%state
43  clr = 1.d-4
44  if( present( localclr ) ) clr=localclr
45  if( present( ctpos ) ) then
46  r(:)= ctpos
47  else
48  call getelementcenter( etype, r(:) )
49  endif
50 
51  tol = 1.0d0
52  do count=1,100
53  call getshapefunc( etype, r, sfunc )
54  xyz_out = matmul( elemt(1:3,1:nn), sfunc )
55  dxyz(1:3) = xyz_out(1:3) - xyz(1:3)
56  dist_last = dot_product( dxyz, dxyz(:) )
57 
58  call tangentbase( etype, nn, r, elemt(1:3,1:nn), tangent )
59  call curvature( etype, nn, r, elemt(1:3,1:nn), curv )
60 
61  ! dF(1:2)
62  df(1:2) = -matmul( dxyz(:), tangent(:,:) )
63  ! d2F(1:2,1:2)
64  d2f(1,1)= dot_product( tangent(:,1), tangent(:,1) ) - dot_product( dxyz, curv(:,1,1) )
65  d2f(1,2)= dot_product( tangent(:,1), tangent(:,2) ) - dot_product( dxyz, curv(:,1,2) )
66  d2f(2,1)= dot_product( tangent(:,2), tangent(:,1) ) - dot_product( dxyz, curv(:,2,1) )
67  d2f(2,2)= dot_product( tangent(:,2), tangent(:,2) ) - dot_product( dxyz, curv(:,2,2) )
68 
69  ! inverse of d2F
70  determ = d2f(1,1)*d2f(2,2) - d2f(1,2)*d2f(2,1)
71  if( determ==0.d0 ) stop "Math error in contact searching"
72  inverse(1,1) = d2f(2,2) / determ
73  inverse(2,2) = d2f(1,1) / determ
74  inverse(1,2) = -d2f(1,2) / determ
75  inverse(2,1) = -d2f(2,1) / determ
76  dr=matmul(inverse,df)
77 
78  tol=dot_product(dr,dr)
79  if( dsqrt(tol)> 3.d0 ) then ! too far away
80  r= -100.d0; exit
81  endif
82 
83  factor = 1.d0
84  do order=1,10
85  r_tmp(1:2) = r(1:2) + factor*dr(1:2)
86  call getshapefunc( etype, r_tmp, sfunc )
87  xyz_out(1:3) = matmul( elemt(1:3,1:nn), sfunc(:) )
88  dxyz(1:3) = xyz(1:3)-xyz_out(:)
89  dist_now = dot_product( dxyz, dxyz )
90  if(dist_now <= dist_last) exit
91  factor = factor*0.7d0
92  enddo
93  r(1:2) = r_tmp(1:2)
94 
95  if( tol<eps ) exit
96  enddo
97 
98  isin = .false.
99  cstate%state = contactfree
100  if( isinsideelement( etype, r, clr )>=0 ) then
101  dxyz(:)=xyz_out(:)-xyz(:)
102  normal(:) = surfacenormal( etype, nn, r, elemt(1:3,1:nn) )
103  normal(:) = normal(:)/dsqrt( dot_product(normal, normal) )
104  do count = 1,3
105  if( dabs(normal(count))<1.d-10 ) normal(count) =0.d0
106  if( dabs(1.d0-dabs(normal(count)))<1.d-10 ) normal(count) =sign(1.d0, normal(count))
107  enddo
108  cstate%distance = dot_product( dxyz, normal )
109 
110  if( cstate%interference_flag == c_if_slave) then
111  if( cstate%init_pos == 0.d0 .and. cstate%distance < cstate%end_pos) then
112  cstate%init_pos = cstate%distance
113  cstate%time_factor = (cstate%end_pos - cstate%distance) / cstate%time_factor
114  cstate%shrink_factor = cstate%distance
115  end if
116  end if
117 
118  if( cstate%interference_flag == 0)then ! not shrink-node
119  if( cstate%distance < distclr*reflen .and. cstate%distance > -5.0d-01*reflen ) isin = .true.
120  else
121  if( cstate%distance < cstate%shrink_factor + distclr*reflen ) isin = .true.
122  end if
123 
124  if( isin ) then
125  if( initstate== contactfree ) then
126  cstate%state = contactstick
127  else
128  cstate%state = initstate
129  endif
130  cstate%gpos(:)=xyz_out(:)
131  cstate%lpos(1:2)=r(:)
132  cstate%direction(:) = normal(:)
133  cstate%wkdist = cstate%distance
134  endif
135  endif
136  end subroutine project_point2element
137 
139  subroutine project_point2solidelement(xyz,etype,nn,elemt,reflen,cstate,isin,distclr,ctpos,localclr)
140  use m_utilities
141  real(kind=kreal),intent(in) :: xyz(3)
142  integer, intent(in) :: etype
143  integer, intent(in) :: nn
144  real(kind=kreal),intent(in) :: elemt(:,:)
145  real(kind=kreal),intent(in) :: reflen
146  type(tcontactstate),intent(inout) :: cstate
147  logical, intent(out) :: isin
148  real(kind=kreal), intent(in) :: distclr
149  real(kind=kreal), optional :: ctpos(3)
150  real(kind=kreal), optional :: localclr
151 
152  integer :: count,order, initstate
153  real(kind=kreal) :: inverse(3,3)
154  real(kind=kreal) :: sfunc(nn), deriv(nn,3)
155  real(kind=kreal) :: r(3), dr(3), r_tmp(3) ! natural coordinate
156  real(kind=kreal) :: xyz_out(3) ! curr. projection position
157  real(kind=kreal) :: dist_last,dist_now, dxyz(3) ! dist between the point and its projection
158  real(kind=kreal),parameter :: eps = 1.0d-8
159  real(kind=kreal) :: clr, tol, factor
160 
161  initstate = cstate%state
162  clr = 1.d-4
163  if( present( localclr ) ) clr=localclr
164  if( present( ctpos ) ) then
165  r(:)= ctpos
166  else
167  call getelementcenter( etype, r(:) )
168  endif
169 
170  tol = 1.0d0
171  do count=1,100
172  call getshapefunc( etype, r, sfunc )
173  xyz_out = matmul( elemt(1:3,1:nn), sfunc(1:nn) )
174  dxyz(1:3) = xyz_out(1:3) - xyz(1:3)
175  dist_last = dot_product( dxyz, dxyz(:) )
176 
177  call getshapederiv( etype, r, deriv )
178  inverse(1:3,1:3) = matmul(elemt(1:3,1:nn),deriv(1:nn,1:3))
179  call calinverse(3, inverse(1:3,1:3))
180  dr(1:3) = -matmul(inverse(1:3,1:3),dxyz(1:3))
181 
182  tol=dot_product(dr,dr)
183  if( count > 1 .and. dsqrt(tol)> 3.d0 ) then ! too far away
184  r= -100.d0; exit
185  endif
186 
187  factor = 1.d0
188  do order=1,10
189  r_tmp(1:3) = r(1:3) + factor*dr(1:3)
190  call getshapefunc( etype, r_tmp, sfunc )
191  xyz_out(1:3) = matmul( elemt(1:3,1:nn), sfunc(1:nn) )
192  dxyz(1:3) = xyz(1:3)-xyz_out(1:3)
193  dist_now = dot_product( dxyz, dxyz )
194  if(dist_now <= dist_last) exit
195  factor = factor*0.7d0
196  enddo
197  r(1:3) = r_tmp(1:3)
198 
199  if( tol<eps ) exit
200  enddo
201 
202  isin = .false.
203  cstate%state = contactfree
204  if( isinside3delement( etype, r, clr )>=0 ) then
205  dxyz(:)=xyz_out(:)-xyz(:)
206  cstate%distance = dsqrt(dot_product( dxyz, dxyz ))
207 
208  if( cstate%distance < distclr ) isin = .true.
209 
210  if( isin ) then
211  if( initstate== contactfree ) then
212  cstate%state = contactstick
213  else
214  cstate%state = initstate
215  endif
216  cstate%gpos(:)=xyz_out(:)
217  cstate%direction(:) = (/1.d0,0.d0,0.d0/)
218  cstate%lpos(1:3)=r(:)
219  endif
220  endif
221  end subroutine
222 
225  subroutine project_point2surfelement(xyz, surf, currpos, cstate, isin, distclr, ctpos, localclr, smoothing)
226  real(kind=kreal), intent(in) :: xyz(3)
227  type(tsurfelement), intent(in) :: surf
228  real(kind=kreal), intent(in) :: currpos(:)
229  type(tcontactstate), intent(inout) :: cstate
230  logical, intent(out) :: isin
231  real(kind=kreal), intent(in) :: distclr
232  real(kind=kreal), optional, intent(in) :: ctpos(2)
233  real(kind=kreal), optional, intent(in) :: localclr
234  integer(kind=kint), intent(in) :: smoothing
235 
236  integer(kind=kint) :: nn, j, iss, etype_use, nn_use
237  real(kind=kreal) :: elem(3, l_max_elem_node)
238  real(kind=kreal) :: elem_tri3n(3, 3)
239 
240  ! Extract element coordinates from surf structure
241  nn = size(surf%nodes)
242  do j = 1, nn
243  iss = surf%nodes(j)
244  elem(1:3, j) = currpos(3*iss-2:3*iss)
245  enddo
246 
247  if (smoothing == kcsnagata) then
248  if (surf%etype == fe_tri3n) then
249  elem_tri3n = elem(1:3, 1:3)
250  call reorder_tri3n_to_tri6n(elem_tri3n, surf%intermediate_points, elem(1:3,1:6))
251  etype_use = fe_tri6n
252  nn_use = 6
253  else if (surf%etype == fe_quad4n) then
254  do j = 1, nn
255  elem(1:3, nn+j) = surf%intermediate_points(1:3, j)
256  enddo
257  etype_use = fe_quad8n
258  nn_use = 8
259  else
260  etype_use = surf%etype
261  nn_use = nn
262  endif
263  else
264  etype_use = surf%etype
265  nn_use = nn
266  endif
267 
268  call project_point2element(xyz, etype_use, nn_use, elem, surf%reflen, cstate, &
269  isin, distclr, ctpos, localclr)
270 
271  end subroutine project_point2surfelement
272 
274  subroutine getmetrictensor( pos, etype, ele, tensor )
275  real(kind=kreal), intent(in) :: pos(2)
276  integer, intent(in) :: etype
277  real(kind=kreal), intent(in) :: ele(:,:)
278  real(kind=kreal), intent(out) :: tensor(2,2)
279 
280  integer :: nn
281  real(kind=kreal) :: tangent(3,2)
282  nn= getnumberofnodes(etype)
283  call tangentbase( etype, nn, pos, ele, tangent )
284  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
285  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
286  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
287  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
288  end subroutine
289 
292  subroutine dispincrematrix( pos, etype, nnode, ele, tangent, tensor, matrix )
293  real(kind=kreal), intent(in) :: pos(2)
294  integer, intent(in) :: etype
295  integer, intent(in) :: nnode
296  real(kind=kreal), intent(in) :: ele(3,nnode)
297  real(kind=kreal), intent(out) :: tangent(3,2)
298  real(kind=kreal), intent(out) :: tensor(2,2)
299  real(kind=kreal), intent(out) :: matrix(2,nnode*3+3)
300 
301  integer :: i,j
302  real(kind=kreal) :: det
303  real(kind=kreal) :: shapefunc(nnode), t1(nnode*3+3), t2(nnode*3+3)
304  call tangentbase( etype, nnode, pos, ele, tangent )
305  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
306  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
307  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
308  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
309  det = tensor(1,1)*tensor(2,2)-tensor(1,2)*tensor(2,1)
310  if( det==0.d0 ) stop "Error in calculate DispIncreMatrix"
311  ! inverse(1,1) = tensor(2,2)/det
312  ! inverse(1,2) = -tensor(1,2)/det
313  ! inverse(2,1) = -tensor(2,1)/det
314  ! inverse(2,2) = tensor(1,1)/det
315 
316  call getshapefunc( etype, pos(:), shapefunc )
317  forall( j=1:3 )
318  t1( j ) = tangent(j,1)
319  t2( j ) = tangent(j,2)
320  end forall
321  forall( i=1:nnode, j=1:3 )
322  t1( i*3+j ) = -tangent(j,1)*shapefunc(i)
323  t2( i*3+j ) = -tangent(j,2)*shapefunc(i)
324  end forall
325  !matrix( 1:2,: ) = matmul( inverse(:,:), matrix )
326  matrix(1,:) = (tensor(2,2)*t1(:)-tensor(1,2)*t2(:))/det
327  matrix(2,:) = (tensor(1,1)*t2(:)-tensor(2,1)*t1(:))/det
328  tangent(:,1) = tangent(:,1)/dsqrt(dot_product(tangent(:,1),tangent(:,1)))
329  tangent(:,2) = tangent(:,2)/dsqrt(dot_product(tangent(:,2),tangent(:,2)))
330  end subroutine
331 
333  subroutine update_direction( nslave, contact, currpos )
334  integer, intent(in) :: nslave
335  type( tcontact ), intent(inout) :: contact
336  real(kind=kreal), intent(in) :: currpos(:)
337 
338  integer(kind=kint) :: slave, sid0, etype, iSS
339  real(kind=kreal) :: coord(3)
340  logical :: isin
341  type(tcontactstate) :: cstate_tmp
342 
343  slave = contact%slave(nslave)
344  coord(:) = currpos(3*slave-2:3*slave)
346  sid0 = contact%states(nslave)%surface
347 
348  cstate_tmp = contact%states(nslave)
349  call project_point2surfelement( coord, contact%master(sid0), currpos, &
350  cstate_tmp, isin, contact%cparam%DISTCLR_NOCHECK, &
351  contact%states(nslave)%lpos, contact%cparam%CLR_SAME_ELEM, &
352  smoothing=contact%smoothing )
353 
354  if( isin ) then
355  etype = contact%master(sid0)%etype
356  iss = isinsideelement( etype, cstate_tmp%lpos, contact%cparam%CLR_CAL_NORM )
357  if( iss>0 ) &
358  call cal_node_normal( cstate_tmp%surface, iss, contact%master, currpos, &
359  cstate_tmp%lpos, cstate_tmp%direction(:) )
360  endif
361 
362  contact%states(nslave)%direction = cstate_tmp%direction
363 
364  end subroutine
365 
367  subroutine cal_node_normal( csurf, isin, surf, currpos, lpos, normal )
368  use elementinfo, only:getvertexcoord, surfacenormal
369  integer, intent(in) :: csurf
370  integer, intent(in) :: isin
371  type(tsurfelement), intent(in) :: surf(:)
372  real(kind=kreal), intent(in) :: currpos(:)
373  real(kind=kreal), intent(in) :: lpos(:)
374  real(kind=kreal), intent(out) :: normal(3)
375  integer(kind=kint) :: cnode, i, j, cnt, nd1, gn, etype, iss, nn,cgn
376  real(kind=kreal) :: cnpos(2), elem(3, l_max_elem_node)
377  integer(kind=kint) :: cnode1, cnode2, gn1, gn2, nsurf, cgn1, cgn2, isin_n
378  real(kind=kreal) :: x, normal_n(3), lpos_n(2)
379 
380  if( 1 <= isin .and. isin <= 4 ) then ! corner
381  cnode = isin
382  gn = surf(csurf)%nodes(cnode)
383  etype = surf(csurf)%etype
384  call getvertexcoord( etype, cnode, cnpos )
385  nn = size( surf(csurf)%nodes )
386  do j=1,nn
387  iss = surf(csurf)%nodes(j)
388  elem(1:3,j)=currpos(3*iss-2:3*iss)
389  enddo
390  normal = surfacenormal( etype, nn, cnpos, elem )
391  cnt = 1
392  do i=1,surf(csurf)%n_neighbor
393  nd1 = surf(csurf)%neighbor(i)
394  nn = size( surf(nd1)%nodes )
395  etype = surf(nd1)%etype
396  cgn = 0
397  do j=1,nn
398  iss = surf(nd1)%nodes(j)
399  elem(1:3,j)=currpos(3*iss-2:3*iss)
400  if( iss==gn ) cgn=j
401  enddo
402  if( cgn>0 ) then
403  call getvertexcoord( etype, cgn, cnpos )
404  !normal = normal+SurfaceNormal( etype, nn, cnpos, elem )
405  normal_n = surfacenormal( etype, nn, cnpos, elem )
406  normal = normal+normal_n
407  cnt = cnt+1
408  endif
409  enddo
410  !normal = normal/cnt !!-???
411  elseif( 12 <= isin .and. isin <= 41 ) then ! edge
412  cnode1 = isin / 10
413  cnode2 = mod(isin, 10)
414  gn1 = surf(csurf)%nodes(cnode1)
415  gn2 = surf(csurf)%nodes(cnode2)
416  etype = surf(csurf)%etype
417  nn = size( surf(csurf)%nodes )
418  do j=1,nn
419  iss = surf(csurf)%nodes(j)
420  elem(1:3,j)=currpos(3*iss-2:3*iss)
421  enddo
422  normal = surfacenormal( etype, nn, lpos, elem )
423  select case (etype)
424  case (fe_tri3n, fe_tri6n, fe_tri6nc)
425  if ( isin==12 ) then; x=lpos(2)-lpos(1)
426  elseif( isin==23 ) then; x=1.d0-2.d0*lpos(2)
427  elseif( isin==31 ) then; x=2.d0*lpos(1)-1.d0
428  else; stop "Error: cal_node_normal: invalid isin"
429  endif
430  case (fe_quad4n, fe_quad8n)
431  if ( isin==12 ) then; x=lpos(1)
432  elseif( isin==23 ) then; x=lpos(2)
433  elseif( isin==34 ) then; x=-lpos(1)
434  elseif( isin==41 ) then; x=-lpos(2)
435  else; stop "Error: cal_node_normal: invalid isin"
436  endif
437  end select
438  ! find neighbor surf that includes cnode1 and cnode2
439  nsurf = 0
440  neib_loop: do i=1, surf(csurf)%n_neighbor
441  nd1 = surf(csurf)%neighbor(i)
442  nn = size( surf(nd1)%nodes )
443  etype = surf(nd1)%etype
444  cgn1 = 0
445  cgn2 = 0
446  do j=1,nn
447  iss = surf(nd1)%nodes(j)
448  elem(1:3,j)=currpos(3*iss-2:3*iss)
449  if( iss==gn1 ) cgn1=j
450  if( iss==gn2 ) cgn2=j
451  enddo
452  if( cgn1>0 .and. cgn2>0 ) then
453  nsurf = nd1
454  isin_n = 10*cgn2 + cgn1
455  x = -x
456  select case (etype)
457  case (fe_tri3n, fe_tri6n, fe_tri6nc)
458  if ( isin_n==12 ) then; lpos_n(1)=0.5d0*(1.d0-x); lpos_n(2)=0.5d0*(1.d0+x)
459  elseif( isin_n==23 ) then; lpos_n(1)=0.d0; lpos_n(2)=0.5d0*(1.d0-x)
460  elseif( isin_n==31 ) then; lpos_n(1)=0.5d0*(1.d0+x); lpos_n(2)=0.d0
461  else; stop "Error: cal_node_normal: invalid isin_n"
462  endif
463  case (fe_quad4n, fe_quad8n)
464  if ( isin_n==12 ) then; lpos_n(1)= x; lpos_n(2)=-1.d0
465  elseif( isin_n==23 ) then; lpos_n(1)= 1.d0; lpos_n(2)= x
466  elseif( isin_n==34 ) then; lpos_n(1)=-x; lpos_n(2)= 1.d0
467  elseif( isin_n==41 ) then; lpos_n(1)=-1.d0; lpos_n(2)=-x
468  else; stop "Error: cal_node_normal: invalid isin_n"
469  endif
470  end select
471  !normal = normal + SurfaceNormal( etype, nn, lpos_n, elem )
472  normal_n = surfacenormal( etype, nn, lpos_n, elem )
473  normal = normal+normal_n
474  exit neib_loop
475  endif
476  enddo neib_loop
477  !if( nsurf==0 ) write(0,*) "Warning: cal_node_normal: neighbor surf not found"
478  !normal = normal/2
479  endif
480  normal = normal/ dsqrt( dot_product( normal, normal ) )
481  end subroutine cal_node_normal
482 
483 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: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::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:578
elementinfo::getvertexcoord
subroutine getvertexcoord(fetype, cnode, localcoord)
Get the natural coord of a vertex node.
Definition: element.f90:1143
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:933
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: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
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:140
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:226
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:334
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:368
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:293
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:275
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