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,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,reflen,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  real(kind=kreal),intent(in) :: reflen
144  type(tcontactstate),intent(inout) :: cstate
145  logical, intent(out) :: isin
146  real(kind=kreal), intent(in) :: distclr
147  real(kind=kreal), optional :: ctpos(3)
148  real(kind=kreal), optional :: localclr
149 
150  integer :: count,order
151  real(kind=kreal) :: inverse(3,3)
152  real(kind=kreal) :: sfunc(nn), deriv(nn,3)
153  real(kind=kreal) :: r(3), dr(3), r_tmp(3) ! natural coordinate
154  real(kind=kreal) :: xyz_out(3) ! curr. projection position
155  real(kind=kreal) :: dist_last,dist_now, dxyz(3) ! dist between the point and its projection
156  real(kind=kreal),parameter :: eps = 1.0d-8
157  real(kind=kreal) :: clr, tol, factor
158 
159  clr = 1.d-4
160  if( present( localclr ) ) clr=localclr
161  if( present( ctpos ) ) then
162  r(:)= ctpos
163  else
164  call getelementcenter( etype, r(:) )
165  endif
166 
167  tol = 1.0d0
168  do count=1,100
169  call getshapefunc( etype, r, sfunc )
170  xyz_out = matmul( elemt(1:3,1:nn), sfunc(1:nn) )
171  dxyz(1:3) = xyz_out(1:3) - xyz(1:3)
172  dist_last = dot_product( dxyz, dxyz(:) )
173 
174  call getshapederiv( etype, r, deriv )
175  inverse(1:3,1:3) = matmul(elemt(1:3,1:nn),deriv(1:nn,1:3))
176  call calinverse(3, inverse(1:3,1:3))
177  dr(1:3) = -matmul(inverse(1:3,1:3),dxyz(1:3))
178 
179  tol=dot_product(dr,dr)
180  if( count > 1 .and. dsqrt(tol)> 3.d0 ) then ! too far away
181  r= -100.d0; exit
182  endif
183 
184  factor = 1.d0
185  do order=1,10
186  r_tmp(1:3) = r(1:3) + factor*dr(1:3)
187  call getshapefunc( etype, r_tmp, sfunc )
188  xyz_out(1:3) = matmul( elemt(1:3,1:nn), sfunc(1:nn) )
189  dxyz(1:3) = xyz(1:3)-xyz_out(1:3)
190  dist_now = dot_product( dxyz, dxyz )
191  if(dist_now <= dist_last) exit
192  factor = factor*0.7d0
193  enddo
194  r(1:3) = r_tmp(1:3)
195 
196  if( tol<eps ) exit
197  enddo
198 
199  isin = .false.
200  if( isinside3delement( etype, r, clr )>=0 ) then
201  dxyz(:)=xyz_out(:)-xyz(:)
202  cstate%distance = dsqrt(dot_product( dxyz, dxyz ))
203 
204  if( cstate%distance < distclr ) isin = .true.
205 
206  if( isin ) then
207  if( cstate%state == contactfree ) cstate%state = contactstick
208  cstate%gpos(:)=xyz_out(:)
209  cstate%direction(:) = (/1.d0,0.d0,0.d0/)
210  cstate%lpos(1:3)=r(:)
211  endif
212  endif
213  end subroutine
214 
217  subroutine project_point2surfelement(xyz, surf, currpos, cstate, isin, distclr, penclr, ctpos, localclr, smoothing)
218  real(kind=kreal), intent(in) :: xyz(3)
219  type(tsurfelement), intent(in) :: surf
220  real(kind=kreal), intent(in) :: currpos(:)
221  type(tcontactstate), intent(inout) :: cstate
222  logical, intent(out) :: isin
223  real(kind=kreal), intent(in) :: distclr
224  real(kind=kreal), intent(in) :: penclr
225  real(kind=kreal), optional, intent(in) :: ctpos(2)
226  real(kind=kreal), optional, intent(in) :: localclr
227  integer(kind=kint), intent(in) :: smoothing
228 
229  integer(kind=kint) :: nn, j, iss, etype_use, nn_use
230  real(kind=kreal) :: elem(3, l_max_elem_node)
231  real(kind=kreal) :: elem_tri3n(3, 3)
232 
233  ! Extract element coordinates from surf structure
234  nn = size(surf%nodes)
235  do j = 1, nn
236  iss = surf%nodes(j)
237  elem(1:3, j) = currpos(3*iss-2:3*iss)
238  enddo
239 
240  if (smoothing == kcsnagata) then
241  if (surf%etype == fe_tri3n) then
242  elem_tri3n = elem(1:3, 1:3)
243  call reorder_tri3n_to_tri6n(elem_tri3n, surf%intermediate_points, elem(1:3,1:6))
244  etype_use = fe_tri6n
245  nn_use = 6
246  else if (surf%etype == fe_quad4n) then
247  do j = 1, nn
248  elem(1:3, nn+j) = surf%intermediate_points(1:3, j)
249  enddo
250  etype_use = fe_quad8n
251  nn_use = 8
252  else
253  etype_use = surf%etype
254  nn_use = nn
255  endif
256  else
257  etype_use = surf%etype
258  nn_use = nn
259  endif
260 
261  call project_point2element(xyz, etype_use, nn_use, elem, surf%reflen, cstate, &
262  isin, distclr, penclr, ctpos, localclr)
263 
264  end subroutine project_point2surfelement
265 
273  subroutine project_point2surfelement_ss( coord, master_surf, sSurf, ncoord, currpos, &
274  cstate, isin, distclr, penclr, ctpos, localclr )
275  real(kind=kreal), intent(in) :: coord(3)
276  type(tsurfelement), intent(in) :: master_surf
277  type(tcontactsurf), intent(in) :: sSurf
278  real(kind=kreal), intent(in) :: ncoord(2)
279  real(kind=kreal), intent(in) :: currpos(:)
280  type(tcontactstate), intent(inout) :: cstate
281  logical, intent(out) :: isin
282  real(kind=kreal), intent(in) :: distclr
283  real(kind=kreal), intent(in) :: penclr
284  real(kind=kreal), optional, intent(in) :: ctpos(2)
285  real(kind=kreal), optional, intent(in) :: localclr
286 
287  integer(kind=kint) :: nnode_s, j, islave
288  real(kind=kreal) :: slave_pos(3, l_max_elem_node), normal(3)
289 
290  ! kcsNONE: the mortar integral is built on the flat master facets. SMOOTHING= is a
291  ! NODE-SURF feature and is not combined with MORTAR=YES.
292  call project_point2surfelement( coord, master_surf, currpos, cstate, isin, distclr, penclr, ctpos, localclr, kcsnone )
293  if( .not. isin ) return
294 
295  nnode_s = size(ssurf%nodes)
296  do j = 1, nnode_s
297  islave = ssurf%nodes(j)
298  slave_pos(:,j) = currpos(3*islave-2:3*islave)
299  enddo
300  normal(:) = - surfacenormal( ssurf%etype, nnode_s, ncoord, slave_pos )
301  cstate%direction(:) = normal(:) / dsqrt( dot_product(normal, normal) )
302  end subroutine project_point2surfelement_ss
303 
305  subroutine getmetrictensor( pos, etype, ele, tensor )
306  real(kind=kreal), intent(in) :: pos(2)
307  integer, intent(in) :: etype
308  real(kind=kreal), intent(in) :: ele(:,:)
309  real(kind=kreal), intent(out) :: tensor(2,2)
310 
311  integer :: nn
312  real(kind=kreal) :: tangent(3,2)
313  nn= getnumberofnodes(etype)
314  call tangentbase( etype, nn, pos, ele, tangent )
315  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
316  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
317  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
318  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
319  end subroutine
320 
323  subroutine dispincrematrix( pos, etype, nnode, ele, tangent, tensor, matrix )
324  real(kind=kreal), intent(in) :: pos(2)
325  integer, intent(in) :: etype
326  integer, intent(in) :: nnode
327  real(kind=kreal), intent(in) :: ele(3,nnode)
328  real(kind=kreal), intent(out) :: tangent(3,2)
329  real(kind=kreal), intent(out) :: tensor(2,2)
330  real(kind=kreal), intent(out) :: matrix(2,nnode*3+3)
331 
332  integer :: i,j
333  real(kind=kreal) :: det
334  real(kind=kreal) :: shapefunc(nnode), t1(nnode*3+3), t2(nnode*3+3)
335  call tangentbase( etype, nnode, pos, ele, tangent )
336  tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
337  tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
338  tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
339  tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
340  det = tensor(1,1)*tensor(2,2)-tensor(1,2)*tensor(2,1)
341  if( det==0.d0 ) stop "Error in calculate DispIncreMatrix"
342  ! inverse(1,1) = tensor(2,2)/det
343  ! inverse(1,2) = -tensor(1,2)/det
344  ! inverse(2,1) = -tensor(2,1)/det
345  ! inverse(2,2) = tensor(1,1)/det
346 
347  call getshapefunc( etype, pos(:), shapefunc )
348  forall( j=1:3 )
349  t1( j ) = tangent(j,1)
350  t2( j ) = tangent(j,2)
351  end forall
352  forall( i=1:nnode, j=1:3 )
353  t1( i*3+j ) = -tangent(j,1)*shapefunc(i)
354  t2( i*3+j ) = -tangent(j,2)*shapefunc(i)
355  end forall
356  !matrix( 1:2,: ) = matmul( inverse(:,:), matrix )
357  matrix(1,:) = (tensor(2,2)*t1(:)-tensor(1,2)*t2(:))/det
358  matrix(2,:) = (tensor(1,1)*t2(:)-tensor(2,1)*t1(:))/det
359  tangent(:,1) = tangent(:,1)/dsqrt(dot_product(tangent(:,1),tangent(:,1)))
360  tangent(:,2) = tangent(:,2)/dsqrt(dot_product(tangent(:,2),tangent(:,2)))
361  end subroutine
362 
364  subroutine update_direction( nslave, contact, currpos )
365  integer, intent(in) :: nslave
366  type( tcontact ), intent(inout) :: contact
367  real(kind=kreal), intent(in) :: currpos(:)
368 
369  integer(kind=kint) :: slave, sid0, etype, iSS
370  real(kind=kreal) :: coord(3)
371  logical :: isin
372  type(tcontactstate) :: cstate_tmp
373 
374  slave = contact%slave(nslave)
375  coord(:) = currpos(3*slave-2:3*slave)
377  sid0 = contact%states(nslave)%surface
378 
379  cstate_tmp = contact%states(nslave)
380  call project_point2surfelement( coord, contact%master(sid0), currpos, &
381  cstate_tmp, isin, contact%cparam%DISTCLR_NOCHECK, contact%cparam%PENCLR_NOCHECK, &
382  contact%states(nslave)%lpos, contact%cparam%CLR_SAME_ELEM, &
383  smoothing=contact%smoothing )
384 
385  if( isin ) then
386  etype = contact%master(sid0)%etype
387  iss = isinsideelement( etype, cstate_tmp%lpos, contact%cparam%CLR_CAL_NORM )
388  if( iss>0 .and. contact%smoothing /= kcsnagata ) &
389  call cal_node_normal( cstate_tmp%surface, iss, contact%master, currpos, &
390  cstate_tmp%lpos, cstate_tmp%direction(:) )
391  endif
392 
393  contact%states(nslave)%direction = cstate_tmp%direction
394 
395  end subroutine
396 
398  subroutine cal_node_normal( csurf, isin, surf, currpos, lpos, normal )
399  use elementinfo, only:getvertexcoord, surfacenormal
400  integer, intent(in) :: csurf
401  integer, intent(in) :: isin
402  type(tsurfelement), intent(in) :: surf(:)
403  real(kind=kreal), intent(in) :: currpos(:)
404  real(kind=kreal), intent(in) :: lpos(:)
405  real(kind=kreal), intent(out) :: normal(3)
406  integer(kind=kint) :: cnode, i, j, cnt, nd1, gn, etype, iss, nn,cgn
407  real(kind=kreal) :: cnpos(2), elem(3, l_max_elem_node)
408  integer(kind=kint) :: cnode1, cnode2, gn1, gn2, nsurf, cgn1, cgn2, isin_n
409  real(kind=kreal) :: x, normal_n(3), lpos_n(2)
410 
411  if( 1 <= isin .and. isin <= 4 ) then ! corner
412  cnode = isin
413  gn = surf(csurf)%nodes(cnode)
414  etype = surf(csurf)%etype
415  call getvertexcoord( etype, cnode, cnpos )
416  nn = size( surf(csurf)%nodes )
417  do j=1,nn
418  iss = surf(csurf)%nodes(j)
419  elem(1:3,j)=currpos(3*iss-2:3*iss)
420  enddo
421  normal = surfacenormal( etype, nn, cnpos, elem )
422  cnt = 1
423  do i=1,surf(csurf)%n_neighbor
424  nd1 = surf(csurf)%neighbor(i)
425  nn = size( surf(nd1)%nodes )
426  etype = surf(nd1)%etype
427  cgn = 0
428  do j=1,nn
429  iss = surf(nd1)%nodes(j)
430  elem(1:3,j)=currpos(3*iss-2:3*iss)
431  if( iss==gn ) cgn=j
432  enddo
433  if( cgn>0 ) then
434  call getvertexcoord( etype, cgn, cnpos )
435  !normal = normal+SurfaceNormal( etype, nn, cnpos, elem )
436  normal_n = surfacenormal( etype, nn, cnpos, elem )
437  normal = normal+normal_n
438  cnt = cnt+1
439  endif
440  enddo
441  !normal = normal/cnt !!-???
442  elseif( 12 <= isin .and. isin <= 41 ) then ! edge
443  cnode1 = isin / 10
444  cnode2 = mod(isin, 10)
445  gn1 = surf(csurf)%nodes(cnode1)
446  gn2 = surf(csurf)%nodes(cnode2)
447  etype = surf(csurf)%etype
448  nn = size( surf(csurf)%nodes )
449  do j=1,nn
450  iss = surf(csurf)%nodes(j)
451  elem(1:3,j)=currpos(3*iss-2:3*iss)
452  enddo
453  normal = surfacenormal( etype, nn, lpos, elem )
454  select case (etype)
455  case (fe_tri3n, fe_tri6n, fe_tri6nc)
456  if ( isin==12 ) then; x=lpos(2)-lpos(1)
457  elseif( isin==23 ) then; x=1.d0-2.d0*lpos(2)
458  elseif( isin==31 ) then; x=2.d0*lpos(1)-1.d0
459  else; stop "Error: cal_node_normal: invalid isin"
460  endif
461  case (fe_quad4n, fe_quad8n)
462  if ( isin==12 ) then; x=lpos(1)
463  elseif( isin==23 ) then; x=lpos(2)
464  elseif( isin==34 ) then; x=-lpos(1)
465  elseif( isin==41 ) then; x=-lpos(2)
466  else; stop "Error: cal_node_normal: invalid isin"
467  endif
468  end select
469  ! find neighbor surf that includes cnode1 and cnode2
470  nsurf = 0
471  neib_loop: do i=1, surf(csurf)%n_neighbor
472  nd1 = surf(csurf)%neighbor(i)
473  nn = size( surf(nd1)%nodes )
474  etype = surf(nd1)%etype
475  cgn1 = 0
476  cgn2 = 0
477  do j=1,nn
478  iss = surf(nd1)%nodes(j)
479  elem(1:3,j)=currpos(3*iss-2:3*iss)
480  if( iss==gn1 ) cgn1=j
481  if( iss==gn2 ) cgn2=j
482  enddo
483  if( cgn1>0 .and. cgn2>0 ) then
484  nsurf = nd1
485  isin_n = 10*cgn2 + cgn1
486  x = -x
487  select case (etype)
488  case (fe_tri3n, fe_tri6n, fe_tri6nc)
489  if ( isin_n==12 ) then; lpos_n(1)=0.5d0*(1.d0-x); lpos_n(2)=0.5d0*(1.d0+x)
490  elseif( isin_n==23 ) then; lpos_n(1)=0.d0; lpos_n(2)=0.5d0*(1.d0-x)
491  elseif( isin_n==31 ) then; lpos_n(1)=0.5d0*(1.d0+x); lpos_n(2)=0.d0
492  else; stop "Error: cal_node_normal: invalid isin_n"
493  endif
494  case (fe_quad4n, fe_quad8n)
495  if ( isin_n==12 ) then; lpos_n(1)= x; lpos_n(2)=-1.d0
496  elseif( isin_n==23 ) then; lpos_n(1)= 1.d0; lpos_n(2)= x
497  elseif( isin_n==34 ) then; lpos_n(1)=-x; lpos_n(2)= 1.d0
498  elseif( isin_n==41 ) then; lpos_n(1)=-1.d0; lpos_n(2)=-x
499  else; stop "Error: cal_node_normal: invalid isin_n"
500  endif
501  end select
502  !normal = normal + SurfaceNormal( etype, nn, lpos_n, elem )
503  normal_n = surfacenormal( etype, nn, lpos_n, elem )
504  normal = normal+normal_n
505  exit neib_loop
506  endif
507  enddo neib_loop
508  !if( nsurf==0 ) write(0,*) "Warning: cal_node_normal: neighbor surf not found"
509  !normal = normal/2
510  endif
511  normal = normal/ dsqrt( dot_product( normal, normal ) )
512  end subroutine cal_node_normal
513 
514 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_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:138
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::update_direction
subroutine update_direction(nslave, contact, currpos)
This subroutine tracks down next contact position after a finite slide.
Definition: fstr_contact_geom.f90:365
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:399
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:324
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:218
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:275
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:306
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::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:130
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:98