FrontISTR  6.0.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,penclr,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), intent(in) :: penclr
29  real(kind=kreal), optional :: ctpos(2)
30  real(kind=kreal), optional :: localclr
31 
32  integer :: count,order
33  real(kind=kreal) :: determ, inverse(2,2)
34  real(kind=kreal) :: sfunc(nn), curv(3,2,2)
35  real(kind=kreal) :: r(2), dr(2), r_tmp(2) ! natural coordinate
36  real(kind=kreal) :: xyz_out(3) ! curr. projection position
37  real(kind=kreal) :: dist_last,dist_now, dxyz(3) ! dist between the point and its projection
38  real(kind=kreal) :: tangent(3,2) ! base vectors in tangent space
39  real(kind=kreal) :: df(2),d2f(2,2),normal(3)
40  real(kind=kreal),parameter :: eps = 1.0d-8
41  real(kind=kreal) :: clr, tol, factor
42 
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  ! The contact state itself is left untouched until the candidate is
99  ! accepted, so that a failed attempt does not destroy the stick or slip
100  ! state that a subsequent candidate has to inherit.
101  isin = .false.
102  if( isinsideelement( etype, r, clr )>=0 ) then
103  dxyz(:)=xyz_out(:)-xyz(:)
104  normal(:) = surfacenormal( etype, nn, r, elemt(1:3,1:nn) )
105  normal(:) = normal(:)/dsqrt( dot_product(normal, normal) )
106  do count = 1,3
107  if( dabs(normal(count))<1.d-10 ) normal(count) =0.d0
108  if( dabs(1.d0-dabs(normal(count)))<1.d-10 ) normal(count) =sign(1.d0, normal(count))
109  enddo
110  cstate%distance = dot_product( dxyz, normal )
111 
112  if( cstate%interference_flag == c_if_slave) then
113  if( cstate%init_pos == 0.d0 .and. cstate%distance < cstate%end_pos) then
114  cstate%init_pos = cstate%distance
115  cstate%time_factor = (cstate%end_pos - cstate%distance) / cstate%time_factor
116  cstate%shrink_factor = cstate%distance
117  end if
118  end if
119 
120  if( cstate%interference_flag == 0)then ! not shrink-node
121  if( cstate%distance < distclr*reflen .and. cstate%distance > -penclr*reflen ) isin = .true.
122  else
123  if( cstate%distance < cstate%shrink_factor + distclr*reflen ) isin = .true.
124  end if
125 
126  if( isin ) then
127  if( cstate%state == contactfree ) cstate%state = contactstick
128  cstate%gpos(:)=xyz_out(:)
129  cstate%lpos(1:2)=r(:)
130  cstate%direction(:) = normal(:)
131  cstate%wkdist = cstate%distance
132  endif
133  endif
134  end subroutine project_point2element
135 
137  subroutine project_point2solidelement(xyz,etype,nn,elemt,cstate,isin,distclr,ctpos,localclr)
138  use m_utilities
139  real(kind=kreal),intent(in) :: xyz(3)
140  integer, intent(in) :: etype
141  integer, intent(in) :: nn
142  real(kind=kreal),intent(in) :: elemt(:,:)
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, penclr, 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), intent(in) :: penclr
224  real(kind=kreal), optional, intent(in) :: ctpos(2)
225  real(kind=kreal), optional, intent(in) :: localclr
226  integer(kind=kint), intent(in) :: smoothing
227 
228  integer(kind=kint) :: nn, j, iss, etype_use, nn_use
229  real(kind=kreal) :: elem(3, l_max_elem_node)
230  real(kind=kreal) :: elem_tri3n(3, 3)
231 
232  ! Extract element coordinates from surf structure
233  nn = size(surf%nodes)
234  do j = 1, nn
235  iss = surf%nodes(j)
236  elem(1:3, j) = currpos(3*iss-2:3*iss)
237  enddo
238 
239  if (smoothing == kcsnagata) then
240  if (surf%etype == fe_tri3n) then
241  elem_tri3n = elem(1:3, 1:3)
242  call reorder_tri3n_to_tri6n(elem_tri3n, surf%intermediate_points, elem(1:3,1:6))
243  etype_use = fe_tri6n
244  nn_use = 6
245  else if (surf%etype == fe_quad4n) then
246  do j = 1, nn
247  elem(1:3, nn+j) = surf%intermediate_points(1:3, j)
248  enddo
249  etype_use = fe_quad8n
250  nn_use = 8
251  else
252  etype_use = surf%etype
253  nn_use = nn
254  endif
255  else
256  etype_use = surf%etype
257  nn_use = nn
258  endif
259 
260  call project_point2element(xyz, etype_use, nn_use, elem, surf%reflen, cstate, &
261  isin, distclr, penclr, ctpos, localclr)
262 
263  end subroutine project_point2surfelement
264 
272  subroutine project_point2surfelement_ss( coord, master_surf, sSurf, ncoord, currpos, &
273  cstate, isin, distclr, penclr, ctpos, localclr )
274  real(kind=kreal), intent(in) :: coord(3)
275  type(tsurfelement), intent(in) :: master_surf
276  type(tcontactsurf), intent(in) :: sSurf
277  real(kind=kreal), intent(in) :: ncoord(2)
278  real(kind=kreal), intent(in) :: currpos(:)
279  type(tcontactstate), intent(inout) :: cstate
280  logical, intent(out) :: isin
281  real(kind=kreal), intent(in) :: distclr
282  real(kind=kreal), intent(in) :: penclr
283  real(kind=kreal), optional, intent(in) :: ctpos(2)
284  real(kind=kreal), optional, intent(in) :: localclr
285 
286  integer(kind=kint) :: nnode_s, j, islave
287  real(kind=kreal) :: slave_pos(3, l_max_elem_node), normal(3)
288 
289  ! kcsNONE: the mortar integral is built on the flat master facets. SMOOTHING= is a
290  ! NODE-SURF feature and is not combined with MORTAR=YES.
291  call project_point2surfelement( coord, master_surf, currpos, cstate, isin, distclr, penclr, ctpos, localclr, kcsnone )
292  if( .not. isin ) return
293 
294  nnode_s = size(ssurf%nodes)
295  do j = 1, nnode_s
296  islave = ssurf%nodes(j)
297  slave_pos(:,j) = currpos(3*islave-2:3*islave)
298  enddo
299  normal(:) = - surfacenormal( ssurf%etype, nnode_s, ncoord, slave_pos )
300  cstate%direction(:) = normal(:) / dsqrt( dot_product(normal, normal) )
301  end subroutine project_point2surfelement_ss
302 
304  subroutine getmetrictensor( pos, etype, ele, tensor )
305  real(kind=kreal), intent(in) :: pos(2)
306  integer, intent(in) :: etype
307  real(kind=kreal), intent(in) :: ele(:,:)
308  real(kind=kreal), intent(out) :: tensor(2,2)
309 
310  integer :: nn
311  real(kind=kreal) :: tangent(3,2)
312  nn= getnumberofnodes(etype)
313  call tangentbase( etype, nn, pos, ele, tangent )
314  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
315  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
316  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
317  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
318  end subroutine
319 
322  subroutine dispincrematrix( pos, etype, nnode, ele, tangent, tensor, matrix )
323  real(kind=kreal), intent(in) :: pos(2)
324  integer, intent(in) :: etype
325  integer, intent(in) :: nnode
326  real(kind=kreal), intent(in) :: ele(3,nnode)
327  real(kind=kreal), intent(out) :: tangent(3,2)
328  real(kind=kreal), intent(out) :: tensor(2,2)
329  real(kind=kreal), intent(out) :: matrix(2,nnode*3+3)
330 
331  integer :: i,j
332  real(kind=kreal) :: det
333  real(kind=kreal) :: shapefunc(nnode), t1(nnode*3+3), t2(nnode*3+3)
334  call tangentbase( etype, nnode, pos, ele, tangent )
335  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
336  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
337  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
338  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
339  det = tensor(1,1)*tensor(2,2)-tensor(1,2)*tensor(2,1)
340  if( det==0.d0 ) stop "Error in calculate DispIncreMatrix"
341  ! inverse(1,1) = tensor(2,2)/det
342  ! inverse(1,2) = -tensor(1,2)/det
343  ! inverse(2,1) = -tensor(2,1)/det
344  ! inverse(2,2) = tensor(1,1)/det
345 
346  call getshapefunc( etype, pos(:), shapefunc )
347  forall( j=1:3 )
348  t1( j ) = tangent(j,1)
349  t2( j ) = tangent(j,2)
350  end forall
351  forall( i=1:nnode, j=1:3 )
352  t1( i*3+j ) = -tangent(j,1)*shapefunc(i)
353  t2( i*3+j ) = -tangent(j,2)*shapefunc(i)
354  end forall
355  !matrix( 1:2,: ) = matmul( inverse(:,:), matrix )
356  matrix(1,:) = (tensor(2,2)*t1(:)-tensor(1,2)*t2(:))/det
357  matrix(2,:) = (tensor(1,1)*t2(:)-tensor(2,1)*t1(:))/det
358  tangent(:,1) = tangent(:,1)/dsqrt(dot_product(tangent(:,1),tangent(:,1)))
359  tangent(:,2) = tangent(:,2)/dsqrt(dot_product(tangent(:,2),tangent(:,2)))
360  end subroutine
361 
363  subroutine update_direction( nslave, contact, currpos )
364  integer, intent(in) :: nslave
365  type( tcontact ), intent(inout) :: contact
366  real(kind=kreal), intent(in) :: currpos(:)
367 
368  integer(kind=kint) :: slave, sid0, etype, iSS
369  real(kind=kreal) :: coord(3)
370  logical :: isin
371  type(tcontactstate) :: cstate_tmp
372 
373  slave = contact%slave(nslave)
374  coord(:) = currpos(3*slave-2:3*slave)
376  sid0 = contact%states(nslave)%surface
377 
378  cstate_tmp = contact%states(nslave)
379  call project_point2surfelement( coord, contact%master(sid0), currpos, &
380  cstate_tmp, isin, contact%cparam%DISTCLR_NOCHECK, contact%cparam%PENCLR_NOCHECK, &
381  contact%states(nslave)%lpos, contact%cparam%CLR_SAME_ELEM, &
382  smoothing=contact%smoothing )
383 
384  if( isin ) then
385  etype = contact%master(sid0)%etype
386  iss = isinsideelement( etype, cstate_tmp%lpos, contact%cparam%CLR_CAL_NORM )
387  if( iss>0 .and. contact%smoothing /= kcsnagata ) &
388  call cal_node_normal( cstate_tmp%surface, iss, contact%master, currpos, &
389  cstate_tmp%lpos, cstate_tmp%direction(:) )
390  endif
391 
392  contact%states(nslave)%direction = cstate_tmp%direction
393 
394  end subroutine
395 
397  subroutine cal_node_normal( csurf, isin, surf, currpos, lpos, normal )
398  use elementinfo, only:getvertexcoord, surfacenormal
399  integer, intent(in) :: csurf
400  integer, intent(in) :: isin
401  type(tsurfelement), intent(in) :: surf(:)
402  real(kind=kreal), intent(in) :: currpos(:)
403  real(kind=kreal), intent(in) :: lpos(:)
404  real(kind=kreal), intent(out) :: normal(3)
405  integer(kind=kint) :: cnode, i, j, cnt, nd1, gn, etype, iss, nn,cgn
406  real(kind=kreal) :: cnpos(2), elem(3, l_max_elem_node)
407  integer(kind=kint) :: cnode1, cnode2, gn1, gn2, nsurf, cgn1, cgn2, isin_n
408  real(kind=kreal) :: x, normal_n(3), lpos_n(2)
409 
410  if( 1 <= isin .and. isin <= 4 ) then ! corner
411  cnode = isin
412  gn = surf(csurf)%nodes(cnode)
413  etype = surf(csurf)%etype
414  call getvertexcoord( etype, cnode, cnpos )
415  nn = size( surf(csurf)%nodes )
416  do j=1,nn
417  iss = surf(csurf)%nodes(j)
418  elem(1:3,j)=currpos(3*iss-2:3*iss)
419  enddo
420  normal = surfacenormal( etype, nn, cnpos, elem )
421  cnt = 1
422  do i=1,surf(csurf)%n_neighbor
423  nd1 = surf(csurf)%neighbor(i)
424  nn = size( surf(nd1)%nodes )
425  etype = surf(nd1)%etype
426  cgn = 0
427  do j=1,nn
428  iss = surf(nd1)%nodes(j)
429  elem(1:3,j)=currpos(3*iss-2:3*iss)
430  if( iss==gn ) cgn=j
431  enddo
432  if( cgn>0 ) then
433  call getvertexcoord( etype, cgn, cnpos )
434  !normal = normal+SurfaceNormal( etype, nn, cnpos, elem )
435  normal_n = surfacenormal( etype, nn, cnpos, elem )
436  normal = normal+normal_n
437  cnt = cnt+1
438  endif
439  enddo
440  !normal = normal/cnt !!-???
441  elseif( 12 <= isin .and. isin <= 41 ) then ! edge
442  cnode1 = isin / 10
443  cnode2 = mod(isin, 10)
444  gn1 = surf(csurf)%nodes(cnode1)
445  gn2 = surf(csurf)%nodes(cnode2)
446  etype = surf(csurf)%etype
447  nn = size( surf(csurf)%nodes )
448  do j=1,nn
449  iss = surf(csurf)%nodes(j)
450  elem(1:3,j)=currpos(3*iss-2:3*iss)
451  enddo
452  normal = surfacenormal( etype, nn, lpos, elem )
453  select case (etype)
454  case (fe_tri3n, fe_tri6n, fe_tri6nc)
455  if ( isin==12 ) then; x=lpos(2)-lpos(1)
456  elseif( isin==23 ) then; x=1.d0-2.d0*lpos(2)
457  elseif( isin==31 ) then; x=2.d0*lpos(1)-1.d0
458  else; stop "Error: cal_node_normal: invalid isin"
459  endif
460  case (fe_quad4n, fe_quad8n)
461  if ( isin==12 ) then; x=lpos(1)
462  elseif( isin==23 ) then; x=lpos(2)
463  elseif( isin==34 ) then; x=-lpos(1)
464  elseif( isin==41 ) then; x=-lpos(2)
465  else; stop "Error: cal_node_normal: invalid isin"
466  endif
467  end select
468  ! find neighbor surf that includes cnode1 and cnode2
469  nsurf = 0
470  neib_loop: do i=1, surf(csurf)%n_neighbor
471  nd1 = surf(csurf)%neighbor(i)
472  nn = size( surf(nd1)%nodes )
473  etype = surf(nd1)%etype
474  cgn1 = 0
475  cgn2 = 0
476  do j=1,nn
477  iss = surf(nd1)%nodes(j)
478  elem(1:3,j)=currpos(3*iss-2:3*iss)
479  if( iss==gn1 ) cgn1=j
480  if( iss==gn2 ) cgn2=j
481  enddo
482  if( cgn1>0 .and. cgn2>0 ) then
483  nsurf = nd1
484  isin_n = 10*cgn2 + cgn1
485  x = -x
486  select case (etype)
487  case (fe_tri3n, fe_tri6n, fe_tri6nc)
488  if ( isin_n==12 ) then; lpos_n(1)=0.5d0*(1.d0-x); lpos_n(2)=0.5d0*(1.d0+x)
489  elseif( isin_n==23 ) then; lpos_n(1)=0.d0; lpos_n(2)=0.5d0*(1.d0-x)
490  elseif( isin_n==31 ) then; lpos_n(1)=0.5d0*(1.d0+x); lpos_n(2)=0.d0
491  else; stop "Error: cal_node_normal: invalid isin_n"
492  endif
493  case (fe_quad4n, fe_quad8n)
494  if ( isin_n==12 ) then; lpos_n(1)= x; lpos_n(2)=-1.d0
495  elseif( isin_n==23 ) then; lpos_n(1)= 1.d0; lpos_n(2)= x
496  elseif( isin_n==34 ) then; lpos_n(1)=-x; lpos_n(2)= 1.d0
497  elseif( isin_n==41 ) then; lpos_n(1)=-1.d0; lpos_n(2)=-x
498  else; stop "Error: cal_node_normal: invalid isin_n"
499  endif
500  end select
501  !normal = normal + SurfaceNormal( etype, nn, lpos_n, elem )
502  normal_n = surfacenormal( etype, nn, lpos_n, elem )
503  normal = normal+normal_n
504  exit neib_loop
505  endif
506  enddo neib_loop
507  !if( nsurf==0 ) write(0,*) "Warning: cal_node_normal: neighbor surf not found"
508  !normal = normal/2
509  endif
510  normal = normal/ dsqrt( dot_product( normal, normal ) )
511  end subroutine cal_node_normal
512 
513 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:754
elementinfo::getnumberofnodes
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:132
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:1247
elementinfo::getelementcenter
subroutine getelementcenter(fetype, localcoord)
Return natural coordinate of the center of element.
Definition: element.f90:1158
elementinfo::fe_tri6n
integer, parameter fe_tri6n
Definition: element.f90:71
elementinfo::getshapederiv
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
Definition: element.f90:685
elementinfo::getvertexcoord
subroutine getvertexcoord(fetype, cnode, localcoord)
Get the natural coord of a vertex node.
Definition: element.f90:1289
elementinfo::fe_tri3n
integer, parameter fe_tri3n
Definition: element.f90:70
elementinfo::tangentbase
subroutine tangentbase(fetype, nn, localcoord, elecoord, tangent)
Calculate base vector of tangent space of 3d surface.
Definition: element.f90:1079
elementinfo::fe_quad4n
integer, parameter fe_quad4n
Definition: element.f90:73
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:1182
elementinfo::surfacenormal
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:1016
elementinfo::curvature
subroutine curvature(fetype, nn, localcoord, elecoord, l2ndderiv, normal, curv)
Calculate curvature tensor at a point along 3d surface.
Definition: element.f90:1113
elementinfo::fe_quad8n
integer, parameter fe_quad8n
Definition: element.f90:74
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_point2element
subroutine project_point2element(xyz, etype, nn, elemt, reflen, cstate, isin, distclr, penclr, ctpos, localclr)
This subroutine find the projection of a slave point onto master surface.
Definition: fstr_contact_geom.f90:20
m_fstr_contact_geom::project_point2solidelement
subroutine project_point2solidelement(xyz, etype, nn, elemt, cstate, isin, distclr, ctpos, localclr)
This subroutine find the projection of a slave point onto master surface.
Definition: fstr_contact_geom.f90:138
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:364
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:398
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:323
m_fstr_contact_geom::project_point2surfelement
subroutine project_point2surfelement(xyz, surf, currpos, cstate, isin, distclr, penclr, 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::project_point2surfelement_ss
subroutine project_point2surfelement_ss(coord, master_surf, sSurf, ncoord, currpos, cstate, isin, distclr, penclr, ctpos, localclr)
Mortar projection wrapper. Projects a slave integration point onto a master surface (project_Point2Su...
Definition: fstr_contact_geom.f90:274
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:305
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:270
mcontactdef
This module manages the data structure for contact calculation.
Definition: fstr_contact_def.F90:12
mcontactdef::kcsnone
integer, parameter kcsnone
contact smoothing type
Definition: fstr_contact_def.F90:51
mcontactdef::c_if_slave
integer, parameter c_if_slave
contact interference type
Definition: fstr_contact_def.F90:66
mcontactdef::cgn
real(kind=kreal), save cgn
convergent condition of penetration
Definition: fstr_contact_def.F90:23
mcontactdef::contactfree
integer, parameter contactfree
contact state definition
Definition: fstr_contact_def.F90:33
mcontactdef::contactstick
integer, parameter contactstick
Definition: fstr_contact_def.F90:35
mcontactdef::kcsnagata
integer, parameter kcsnagata
Definition: fstr_contact_def.F90:52
mcontactdef::tcontact
Structure to includes all info needed by contact calculation.
Definition: fstr_contact_def.F90:128
mcontactdef::tcontactstate
This structure records contact status.
Definition: fstr_contact_def.F90:74
mcontactdef::tcontactsurf
Structure to define a slave surface segment of a SURF-SURF (mortar) contact pair.
Definition: fstr_contact_def.F90:99