FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
fstr_contact_def.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 !-------------------------------------------------------------------------------
13 
14  use hecmw
15  use msurfelement
16  use m_contact_lib
18  use bucket_search
19 
20  implicit none
21 
22  include 'fstr_ctrl_util_f.inc'
23 
25  type tcontact
26  ! following contact definition
27  character(len=HECMW_NAME_LEN) :: name
28  integer :: ctype
29  integer :: group
30  character(len=HECMW_NAME_LEN) :: pair_name
31  integer :: surf_id1, surf_id2
32  integer :: surf_id1_sgrp
33  type(tsurfelement), pointer :: master(:)=>null()
34  integer, pointer :: slave(:)=>null()
35  real(kind=kreal) :: fcoeff
36  real(kind=kreal) :: tpenalty
37 
38  ! following algorithm
39  ! -1: not initialized
40  ! 1: TIED-Just rigidly fixed the two surfaces
41  ! 2: GLUED-Distance between the two surfaces to zero and glue them
42  ! 3: SSLID-Small sliding contact( no position but with contact state change)
43  ! 4: FSLID-Finite sliding contact (both changes in contact state and poisition possible)
44  integer :: algtype
45 
46  logical :: mpced
47  logical :: symmetric
48 
49  ! following contact state
50  type(tcontactstate), pointer :: states(:)=>null()
51 
52  type(fstrst_contact_comm) :: comm
53  type(bucketdb) :: master_bktdb
54  end type tcontact
55 
57  logical :: active
58  integer(kind=kint) :: contact2free
59  integer(kind=kint) :: contact2neighbor
60  integer(kind=kint) :: contact2difflpos
61  integer(kind=kint) :: free2contact
62  integer(kind=kint) :: contactnode_previous
63  integer(kind=kint) :: contactnode_current
65 
66  real(kind=kreal), parameter :: clearance = 1.d-4
67  real(kind=kreal), parameter :: clr_same_elem = 5.d-3
68  real(kind=kreal), parameter :: clr_difflpos = 1.d-2
69  real(kind=kreal), parameter :: clr_cal_norm = 1.d-4
70  real(kind=kreal), parameter :: distclr_init = 1.d-6
71  real(kind=kreal), parameter :: distclr_free =-1.d-6
72  real(kind=kreal), parameter :: distclr_nocheck = 1.d0
74  real(kind=kreal), parameter :: box_exp_rate = 1.05d0
75  real(kind=kreal), parameter :: tensile_force =-1.d-8
76 
77  private :: is_mpc_available
78  private :: is_active_contact
79  private :: cal_node_normal
80  private :: track_contact_position
81  private :: reset_contact_force
82 
83  private :: clearance
84  private :: clr_same_elem
85  private :: clr_difflpos
86  private :: clr_cal_norm
87  private :: distclr_free
88  private :: distclr_nocheck
89  private :: box_exp_rate
90  private :: tensile_force
91 
92 contains
93 
95  logical function is_mpc_available( contact )
96  type(tcontact), intent(in) :: contact
97  is_mpc_available = .true.
98  if( contact%fcoeff/=0.d0 ) is_mpc_available = .false.
99  end function
100 
102  subroutine fstr_write_contact( file, contact )
103  integer(kind=kint), intent(in) :: file
104  type(tcontact), intent(in) :: contact
105  integer :: i
106  write(file,*) "CONTACT:", contact%ctype,contact%group,trim(contact%pair_name),contact%fcoeff
107  write(file,*) "---Slave----"
108  write(file,*) "num.slave",size(contact%slave)
109  if( associated(contact%slave) ) then
110  do i=1,size(contact%slave)
111  write(file, *) contact%slave(i)
112  enddo
113  endif
114  write(file,*) "----master---"
115  write(file,*) "num.master",size(contact%master)
116  if( associated(contact%master) ) then
117  do i=1,size(contact%master)
118  call write_surf( file, contact%master(i) )
119  enddo
120  endif
121  end subroutine
122 
124  subroutine fstr_contact_finalize( contact )
125  type(tcontact), intent(inout) :: contact
126  integer :: i
127  if( associated( contact%slave ) ) deallocate(contact%slave)
128  if( associated( contact%master ) ) then
129  do i=1,size( contact%master )
130  call finalize_surf( contact%master(i) )
131  enddo
132  deallocate(contact%master)
133  endif
134  if( associated(contact%states) ) deallocate(contact%states)
135  call fstr_contact_comm_finalize(contact%comm)
136  call bucketdb_finalize( contact%master_bktDB )
137  end subroutine
138 
140  logical function fstr_contact_check( contact, hecMESH )
141  type(tcontact), intent(inout) :: contact
142  type(hecmwst_local_mesh), pointer :: hecmesh
143 
144  integer :: i
145  logical :: isfind
146 
147  fstr_contact_check = .false.
148 
149  ! if contact pair exist?
150  isfind = .false.
151  do i=1,hecmesh%contact_pair%n_pair
152  if( hecmesh%contact_pair%name(i) == contact%pair_name ) then
153  contact%ctype = hecmesh%contact_pair%type(i)
154  contact%surf_id1 = hecmesh%contact_pair%slave_grp_id(i)
155  contact%surf_id2 = hecmesh%contact_pair%master_grp_id(i)
156  contact%surf_id1_sgrp = hecmesh%contact_pair%slave_orisgrp_id(i)
157  isfind = .true.
158  endif
159  enddo
160  if( .not. isfind ) return;
161  if( contact%fcoeff<=0.d0 ) contact%fcoeff=0.d0
162  if( contact%ctype/=1 .and. contact%ctype/=2 ) return
163  if( contact%group<=0 ) return
164 
165  fstr_contact_check = .true.
166  end function
167 
169  logical function fstr_contact_init( contact, hecMESH,myrank )
170  type(tcontact), intent(inout) :: contact
171  type(hecmwst_local_mesh), pointer :: hecmesh
172  integer(kint),intent(in),optional :: myrank
173 
174  integer :: i, j, is, ie, cgrp, nsurf, nslave, ic, ic_type, iss, nn, ii
175  integer :: count,nodeid
176 
177  fstr_contact_init = .false.
178  ! master surface
179  cgrp = contact%surf_id2
180  if( cgrp<=0 ) return
181  is= hecmesh%surf_group%grp_index(cgrp-1) + 1
182  ie= hecmesh%surf_group%grp_index(cgrp )
183 
184  if(present(myrank)) then
185  ! PARA_CONTACT
186  count = 0
187  do i=is,ie
188  ic = hecmesh%surf_group%grp_item(2*i-1)
189  if(hecmesh%elem_ID(ic*2) /= myrank) cycle
190  count = count + 1
191  enddo
192  allocate( contact%master(count) )
193  count = 0
194  do i=is,ie
195  ic = hecmesh%surf_group%grp_item(2*i-1)
196  if(hecmesh%elem_ID(ic*2) /= myrank) cycle
197  count = count + 1
198  nsurf = hecmesh%surf_group%grp_item(2*i)
199  ic_type = hecmesh%elem_type(ic)
200  call initialize_surf( ic, ic_type, nsurf, contact%master(count) )
201  iss = hecmesh%elem_node_index(ic-1)
202  do j=1, size( contact%master(count)%nodes )
203  nn = contact%master(count)%nodes(j)
204  contact%master(count)%nodes(j) = hecmesh%elem_node_item( iss+nn )
205  enddo
206  enddo
207 
208  else
209  ! not PARA_CONTACT
210  allocate( contact%master(ie-is+1) )
211  do i=is,ie
212  ic = hecmesh%surf_group%grp_item(2*i-1)
213  nsurf = hecmesh%surf_group%grp_item(2*i)
214  ic_type = hecmesh%elem_type(ic)
215  call initialize_surf( ic, ic_type, nsurf, contact%master(i-is+1) )
216  iss = hecmesh%elem_node_index(ic-1)
217  do j=1, size( contact%master(i-is+1)%nodes )
218  nn = contact%master(i-is+1)%nodes(j)
219  contact%master(i-is+1)%nodes(j) = hecmesh%elem_node_item( iss+nn )
220  enddo
221  enddo
222 
223  endif
224 
225  call update_surface_reflen( contact%master, hecmesh%node )
226 
227  ! slave surface
228  ! if( contact%ctype==1 ) then
229  cgrp = contact%surf_id1
230  if( cgrp<=0 ) return
231  is= hecmesh%node_group%grp_index(cgrp-1) + 1
232  ie= hecmesh%node_group%grp_index(cgrp )
233  nslave = 0
234  do i=is,ie
235  nodeid = hecmesh%global_node_ID(hecmesh%node_group%grp_item(i))
236  if(present(myrank)) then
237  ! PARA_CONTACT
238  nslave = nslave + 1
239  else
240  ! not PARA_CONTACT
241  if( hecmesh%node_group%grp_item(i) <= hecmesh%nn_internal) then
242  nslave = nslave + 1
243  endif
244  endif
245  enddo
246  allocate( contact%slave(nslave) )
247  allocate( contact%states(nslave) )
248  ii = 0
249  do i=is,ie
250  if(.not.present(myrank)) then
251  ! not PARA_CONTACT
252  if( hecmesh%node_group%grp_item(i) > hecmesh%nn_internal) cycle
253  endif
254  ii = ii + 1
255  contact%slave(ii) = hecmesh%node_group%grp_item(i)
256  contact%states(ii)%state = -1
257  contact%states(ii)%multiplier(:) = 0.d0
258  contact%states(ii)%tangentForce(:) = 0.d0
259  contact%states(ii)%tangentForce1(:) = 0.d0
260  contact%states(ii)%tangentForce_trial(:) = 0.d0
261  contact%states(ii)%tangentForce_final(:) = 0.d0
262  contact%states(ii)%reldisp(:) = 0.d0
263  enddo
264  ! endif
265 
266  ! contact state
267  ! allocate( contact%states(nslave) )
268  do i=1,nslave
269  call contact_state_init( contact%states(i) )
270  enddo
271 
272  ! neighborhood of surface group
273  call update_surface_box_info( contact%master, hecmesh%node )
274  call bucketdb_init( contact%master_bktDB )
275  call update_surface_bucket_info( contact%master, contact%master_bktDB )
276  call find_surface_neighbor( contact%master, contact%master_bktDB )
277 
278  ! initialize contact communication table
279  call fstr_contact_comm_init( contact%comm, hecmesh, 1, nslave, contact%slave )
280 
281  contact%symmetric = .true.
282  fstr_contact_init = .true.
283  end function
284 
286  subroutine clear_contact_state( contact )
287  type(tcontact), intent(inout) :: contact
288  integer :: i
289  if( .not. associated(contact%states) ) return
290  do i=1,size( contact%states )
291  contact%states(i)%state = -1
292  enddo
293  end subroutine
294 
296  logical function is_active_contact( acgrp, contact )
297  integer, intent(in) :: acgrp(:)
298  type(tcontact), intent(in) :: contact
299  if( any( acgrp==contact%group ) ) then
300  is_active_contact = .true.
301  else
302  is_active_contact = .false.
303  endif
304  end function
305 
309  subroutine scan_contact_state( flag_ctAlgo, contact, currpos, currdisp, ndforce, infoCTChange, &
310  nodeID, elemID, is_init, active, mu, B )
311  character(len=9), intent(in) :: flag_ctAlgo
312  type( tcontact ), intent(inout) :: contact
313  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
314  real(kind=kreal), intent(in) :: currpos(:)
315  real(kind=kreal), intent(in) :: currdisp(:)
316  real(kind=kreal), intent(in) :: ndforce(:)
317  integer(kind=kint), intent(in) :: nodeID(:)
318  integer(kind=kint), intent(in) :: elemID(:)
319  logical, intent(in) :: is_init
320  logical, intent(out) :: active
321  real(kind=kreal), intent(in) :: mu
322  real(kind=kreal), optional, target :: b(:)
323 
324  real(kind=kreal) :: distclr
325  integer(kind=kint) :: slave, id, etype
326  integer(kind=kint) :: nn, i, j, iSS, nactive
327  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node )
328  real(kind=kreal) :: nlforce, slforce(3)
329  logical :: isin
330  integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
331  !
332  integer, pointer :: indexMaster(:),indexCand(:)
333  integer :: nMaster,idm,nMasterMax,bktID,nCand
334  logical :: is_cand, is_present_B
335  real(kind=kreal), pointer :: bp(:)
336 
337  if( contact%algtype<=2 ) return
338 
339  if( is_init ) then
340  distclr = distclr_init
341  else
342  distclr = distclr_free
343  endif
344 
345  allocate(contact_surf(size(nodeid)))
346  allocate(states_prev(size(contact%slave)))
347  contact_surf(:) = size(elemid)+1
348  do i = 1, size(contact%slave)
349  states_prev(i) = contact%states(i)%state
350  enddo
351 
352  call update_surface_box_info( contact%master, currpos )
353  call update_surface_bucket_info( contact%master, contact%master_bktDB )
354 
355  ! for gfortran-10: optional parameter seems not allowed within omp parallel
356  is_present_b = present(b)
357  if( is_present_b ) bp => b
358 
359  !$omp parallel do &
360  !$omp& default(none) &
361  !$omp& private(i,slave,slforce,id,nlforce,coord,indexMaster,nMaster,nn,j,iSS,elem,is_cand,idm,etype,isin, &
362  !$omp& bktID,nCand,indexCand) &
363  !$omp& firstprivate(nMasterMax,is_present_B) &
364  !$omp& shared(contact,ndforce,flag_ctAlgo,infoCTChange,currpos,currdisp,mu,nodeID,elemID,Bp,distclr,contact_surf) &
365  !$omp& reduction(.or.:active) &
366  !$omp& schedule(dynamic,1)
367  do i= 1, size(contact%slave)
368  slave = contact%slave(i)
369  if( contact%states(i)%state==contactstick .or. contact%states(i)%state==contactslip ) then
370  slforce(1:3)=ndforce(3*slave-2:3*slave)
371  id = contact%states(i)%surface
372  nlforce = contact%states(i)%multiplier(1)
373 
374  if( nlforce < tensile_force ) then
375  contact%states(i)%state = contactfree
376  contact%states(i)%multiplier(:) = 0.d0
377  write(*,'(A,i10,A,i10,A,e12.3)') "Node",nodeid(slave)," free from contact with element", &
378  elemid(contact%master(id)%eid), " with tensile force ", nlforce
379  cycle
380  endif
381  if( contact%algtype /= contactfslid .or. (.not. is_present_b) ) then ! small slide problem
382  contact_surf(contact%slave(i)) = -id
383  else
384  call track_contact_position( flag_ctalgo, i, contact, currpos, currdisp, mu, infoctchange, nodeid, elemid, bp )
385  if( contact%states(i)%state /= contactfree ) then
386  contact_surf(contact%slave(i)) = -contact%states(i)%surface
387  endif
388  endif
389 
390  else if( contact%states(i)%state==contactfree ) then
391  coord(:) = currpos(3*slave-2:3*slave)
392 
393  ! get master candidates from bucketDB
394  bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
395  ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
396  if (ncand == 0) cycle
397  allocate(indexcand(ncand))
398  call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
399 
400  nmastermax = ncand
401  allocate(indexmaster(nmastermax))
402  nmaster = 0
403 
404  ! narrow down candidates
405  do idm= 1, ncand
406  id = indexcand(idm)
407  if (.not. is_in_surface_box( contact%master(id), coord(1:3), box_exp_rate )) cycle
408  nmaster = nmaster + 1
409  indexmaster(nmaster) = id
410  enddo
411  deallocate(indexcand)
412 
413  if(nmaster == 0) then
414  deallocate(indexmaster)
415  cycle
416  endif
417 
418  do idm = 1,nmaster
419  id = indexmaster(idm)
420  etype = contact%master(id)%etype
421  nn = size( contact%master(id)%nodes )
422  do j=1,nn
423  iss = contact%master(id)%nodes(j)
424  elem(1:3,j)=currpos(3*iss-2:3*iss)
425  enddo
426  call project_point2element( coord,etype,nn,elem,contact%master(id)%reflen,contact%states(i), &
427  isin,distclr,localclr=clearance )
428  if( .not. isin ) cycle
429  contact%states(i)%surface = id
430  contact%states(i)%multiplier(:) = 0.d0
431  iss = isinsideelement( etype, contact%states(i)%lpos, clr_cal_norm )
432  if( iss>0 ) &
433  call cal_node_normal( id, iss, contact%master, currpos, contact%states(i)%lpos, &
434  contact%states(i)%direction(:) )
435  contact_surf(contact%slave(i)) = id
436  write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3,A,3f7.3)') "Node",nodeid(slave)," contact with element", &
437  elemid(contact%master(id)%eid), &
438  " with distance ", contact%states(i)%distance," at ",contact%states(i)%lpos(:), &
439  " along direction ", contact%states(i)%direction
440  exit
441  enddo
442  deallocate(indexmaster)
443  endif
444  enddo
445  !$omp end parallel do
446 
447  call fstr_contact_comm_allreduce_i(contact%comm, contact_surf, hecmw_min)
448  nactive = 0
449  do i = 1, size(contact%slave)
450  if (contact%states(i)%state /= contactfree) then ! any slave in contact
451  if (abs(contact_surf(contact%slave(i))) /= contact%states(i)%surface) then ! that is in contact with other surface
452  contact%states(i)%state = contactfree ! should be freed
453  write(*,'(A,i10,A,i6,A,i6,A)') "Node",nodeid(contact%slave(i)), &
454  " in rank",hecmw_comm_get_rank()," freed due to duplication"
455  else
456  nactive = nactive + 1
457  endif
458  endif
459  if (states_prev(i) == contactfree .and. contact%states(i)%state /= contactfree) then
460  infoctchange%free2contact = infoctchange%free2contact + 1
461  elseif (states_prev(i) /= contactfree .and. contact%states(i)%state == contactfree) then
462  infoctchange%contact2free = infoctchange%contact2free + 1
463  endif
464  enddo
465  active = (nactive > 0)
466  deallocate(contact_surf)
467  deallocate(states_prev)
468  end subroutine scan_contact_state
469 
471  subroutine cal_node_normal( csurf, isin, surf, currpos, lpos, normal )
473  integer, intent(in) :: csurf
474  integer, intent(in) :: isin
475  type(tsurfelement), intent(in) :: surf(:)
476  real(kind=kreal), intent(in) :: currpos(:)
477  real(kind=kreal), intent(in) :: lpos(:)
478  real(kind=kreal), intent(out) :: normal(3)
479  integer(kind=kint) :: cnode, i, j, cnt, nd1, gn, etype, iSS, nn,cgn
480  real(kind=kreal) :: cnpos(2), elem(3, l_max_elem_node )
481  integer(kind=kint) :: cnode1, cnode2, gn1, gn2, nsurf, cgn1, cgn2, isin_n
482  real(kind=kreal) :: x=0, normal_n(3), lpos_n(2)
483 
484  if( 1 <= isin .and. isin <= 4 ) then ! corner
485  cnode = isin
486  gn = surf(csurf)%nodes(cnode)
487  etype = surf(csurf)%etype
488  call getvertexcoord( etype, cnode, cnpos )
489  nn = size( surf(csurf)%nodes )
490  do j=1,nn
491  iss = surf(csurf)%nodes(j)
492  elem(1:3,j)=currpos(3*iss-2:3*iss)
493  enddo
494  normal = surfacenormal( etype, nn, cnpos, elem )
495  cnt = 1
496  do i=1,surf(csurf)%n_neighbor
497  nd1 = surf(csurf)%neighbor(i)
498  nn = size( surf(nd1)%nodes )
499  etype = surf(nd1)%etype
500  cgn = 0
501  do j=1,nn
502  iss = surf(nd1)%nodes(j)
503  elem(1:3,j)=currpos(3*iss-2:3*iss)
504  if( iss==gn ) cgn=j
505  enddo
506  if( cgn>0 ) then
507  call getvertexcoord( etype, cgn, cnpos )
508  !normal = normal+SurfaceNormal( etype, nn, cnpos, elem )
509  normal_n = surfacenormal( etype, nn, cnpos, elem )
510  normal = normal+normal_n
511  cnt = cnt+1
512  endif
513  enddo
514  !normal = normal/cnt !!-???
515  elseif( 12 <= isin .and. isin <= 41 ) then ! edge
516  cnode1 = isin / 10
517  cnode2 = mod(isin, 10)
518  gn1 = surf(csurf)%nodes(cnode1)
519  gn2 = surf(csurf)%nodes(cnode2)
520  etype = surf(csurf)%etype
521  nn = size( surf(csurf)%nodes )
522  do j=1,nn
523  iss = surf(csurf)%nodes(j)
524  elem(1:3,j)=currpos(3*iss-2:3*iss)
525  enddo
526  normal = surfacenormal( etype, nn, lpos, elem )
527  select case (etype)
528  case (fe_tri3n, fe_tri6n, fe_tri6nc)
529  if ( isin==12 ) then; x=lpos(2)-lpos(1)
530  elseif( isin==23 ) then; x=1.d0-2.d0*lpos(2)
531  elseif( isin==31 ) then; x=2.d0*lpos(1)-1.d0
532  else; stop "Error: cal_node_normal: invalid isin"
533  endif
534  case (fe_quad4n, fe_quad8n)
535  if ( isin==12 ) then; x=lpos(1)
536  elseif( isin==23 ) then; x=lpos(2)
537  elseif( isin==34 ) then; x=-lpos(1)
538  elseif( isin==41 ) then; x=-lpos(2)
539  else; stop "Error: cal_node_normal: invalid isin"
540  endif
541  end select
542  ! find neighbor surf that includes cnode1 and cnode2
543  nsurf = 0
544  neib_loop: do i=1, surf(csurf)%n_neighbor
545  nd1 = surf(csurf)%neighbor(i)
546  nn = size( surf(nd1)%nodes )
547  etype = surf(nd1)%etype
548  cgn1 = 0
549  cgn2 = 0
550  do j=1,nn
551  iss = surf(nd1)%nodes(j)
552  elem(1:3,j)=currpos(3*iss-2:3*iss)
553  if( iss==gn1 ) cgn1=j
554  if( iss==gn2 ) cgn2=j
555  enddo
556  if( cgn1>0 .and. cgn2>0 ) then
557  nsurf = nd1
558  isin_n = 10*cgn2 + cgn1
559  x = -x
560  select case (etype)
561  case (fe_tri3n, fe_tri6n, fe_tri6nc)
562  if ( isin_n==12 ) then; lpos_n(1)=0.5d0*(1.d0-x); lpos_n(2)=0.5d0*(1.d0+x)
563  elseif( isin_n==23 ) then; lpos_n(1)=0.d0; lpos_n(2)=0.5d0*(1.d0-x)
564  elseif( isin_n==31 ) then; lpos_n(1)=0.5d0*(1.d0+x); lpos_n(2)=0.d0
565  else; stop "Error: cal_node_normal: invalid isin_n"
566  endif
567  case (fe_quad4n, fe_quad8n)
568  if ( isin_n==12 ) then; lpos_n(1)= x; lpos_n(2)=-1.d0
569  elseif( isin_n==23 ) then; lpos_n(1)= 1.d0; lpos_n(2)= x
570  elseif( isin_n==34 ) then; lpos_n(1)=-x; lpos_n(2)= 1.d0
571  elseif( isin_n==41 ) then; lpos_n(1)=-1.d0; lpos_n(2)=-x
572  else; stop "Error: cal_node_normal: invalid isin_n"
573  endif
574  end select
575  !normal = normal + SurfaceNormal( etype, nn, lpos_n, elem )
576  normal_n = surfacenormal( etype, nn, lpos_n, elem )
577  normal = normal+normal_n
578  exit neib_loop
579  endif
580  enddo neib_loop
581  !if( nsurf==0 ) write(0,*) "Warning: cal_node_normal: neighbor surf not found"
582  !normal = normal/2
583  endif
584  normal = normal/ dsqrt( dot_product( normal, normal ) )
585  end subroutine cal_node_normal
586 
588  subroutine track_contact_position( flag_ctAlgo, nslave, contact, currpos, currdisp, mu, infoCTChange, nodeID, elemID, B )
589  character(len=9), intent(in) :: flag_ctAlgo
590  integer, intent(in) :: nslave
591  type( tcontact ), intent(inout) :: contact
592  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
593  real(kind=kreal), intent(in) :: currpos(:)
594  real(kind=kreal), intent(in) :: currdisp(:)
595  real(kind=kreal), intent(in) :: mu
596  integer(kind=kint), intent(in) :: nodeID(:)
597  integer(kind=kint), intent(in) :: elemID(:)
598  real(kind=kreal), intent(inout) :: b(:)
599 
600  integer(kind=kint) :: slave, sid0, sid, etype
601  integer(kind=kint) :: nn, i, j, iSS
602  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node ), elem0(3, l_max_elem_node )
603  logical :: isin
604  real(kind=kreal) :: opos(2), odirec(3)
605  integer(kind=kint) :: bktID, nCand, idm
606  integer(kind=kint), allocatable :: indexCand(:)
607 
608  sid = 0
609 
610  slave = contact%slave(nslave)
611  coord(:) = currpos(3*slave-2:3*slave)
613  sid0 = contact%states(nslave)%surface
614  opos = contact%states(nslave)%lpos
615  odirec = contact%states(nslave)%direction
616  etype = contact%master(sid0)%etype
617  nn = getnumberofnodes( etype )
618  do j=1,nn
619  iss = contact%master(sid0)%nodes(j)
620  elem(1:3,j)=currpos(3*iss-2:3*iss)
621  elem0(1:3,j)=currpos(3*iss-2:3*iss)-currdisp(3*iss-2:3*iss)
622  enddo
623  call project_point2element( coord,etype,nn,elem,contact%master(sid0)%reflen,contact%states(nslave), &
624  isin,distclr_nocheck,contact%states(nslave)%lpos,clr_same_elem )
625  if( .not. isin ) then
626  do i=1, contact%master(sid0)%n_neighbor
627  sid = contact%master(sid0)%neighbor(i)
628  etype = contact%master(sid)%etype
629  nn = getnumberofnodes( etype )
630  do j=1,nn
631  iss = contact%master(sid)%nodes(j)
632  elem(1:3,j)=currpos(3*iss-2:3*iss)
633  enddo
634  call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
635  isin,distclr_nocheck,localclr=clearance )
636  if( isin ) then
637  contact%states(nslave)%surface = sid
638  exit
639  endif
640  enddo
641  endif
642 
643  if( .not. isin ) then ! such case is considered to rarely or never occur
644  write(*,*) 'Warning: contact moved beyond neighbor elements'
645  ! get master candidates from bucketDB
646  bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
647  ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
648  if (ncand > 0) then
649  allocate(indexcand(ncand))
650  call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
651  do idm= 1, ncand
652  sid = indexcand(idm)
653  if( sid==sid0 ) cycle
654  if( any(sid==contact%master(sid0)%neighbor(:)) ) cycle
655  if (.not. is_in_surface_box( contact%master(sid), coord(1:3), box_exp_rate )) cycle
656  etype = contact%master(sid)%etype
657  nn = size( contact%master(sid)%nodes )
658  do j=1,nn
659  iss = contact%master(sid)%nodes(j)
660  elem(1:3,j)=currpos(3*iss-2:3*iss)
661  enddo
662  call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
663  isin,distclr_nocheck,localclr=clearance )
664  if( isin ) then
665  contact%states(nslave)%surface = sid
666  exit
667  endif
668  enddo
669  deallocate(indexcand)
670  endif
671  endif
672 
673  if( isin ) then
674  if( contact%states(nslave)%surface==sid0 ) then
675  if(any(dabs(contact%states(nslave)%lpos(:)-opos(:)) >= clr_difflpos)) then
676  !$omp atomic
677  infoctchange%contact2difflpos = infoctchange%contact2difflpos + 1
678  endif
679  else
680  write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3)') "Node",nodeid(slave)," move to contact with", &
681  elemid(contact%master(sid)%eid), " with distance ", &
682  contact%states(nslave)%distance," at ",contact%states(nslave)%lpos(:)
683  !$omp atomic
684  infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
685  if( flag_ctalgo=='ALagrange' ) &
686  call reset_contact_force( contact, currpos, nslave, sid0, opos, odirec, b )
687  endif
688  if( flag_ctalgo=='SLagrange' ) call update_tangentforce(etype,nn,elem0,elem,contact%states(nslave))
689  iss = isinsideelement( etype, contact%states(nslave)%lpos, clr_cal_norm )
690  if( iss>0 ) &
691  call cal_node_normal( contact%states(nslave)%surface, iss, contact%master, currpos, &
692  contact%states(nslave)%lpos, contact%states(nslave)%direction(:) )
693  else if( .not. isin ) then
694  write(*,'(A,i10,A)') "Node",nodeid(slave)," move out of contact"
695  contact%states(nslave)%state = contactfree
696  contact%states(nslave)%multiplier(:) = 0.d0
697  endif
698 
699  end subroutine track_contact_position
700 
703  subroutine reset_contact_force( contact, currpos, lslave, omaster, opos, odirec, B )
704  type( tcontact ), intent(inout) :: contact
705  real(kind=kreal), intent(in) :: currpos(:)
706  integer, intent(in) :: lslave
707  integer, intent(in) :: omaster
708  real(kind=kreal), intent(in) :: opos(2)
709  real(kind=kreal), intent(in) :: odirec(3)
710  real(kind=kreal), intent(inout) :: b(:)
711 
712  integer(kind=kint) :: slave, etype, master
713  integer(kind=kint) :: nn, j, iSS
714  real(kind=kreal) :: nrlforce, fcoeff, tangent(3,2)
715  real(kind=kreal) :: elemcrd(3, l_max_elem_node )
716  real(kind=kreal) :: shapefunc(l_max_surface_node)
717  real(kind=kreal) :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
718  real(kind=kreal) :: fric(2), f3(l_max_elem_node*3+3)
719  integer(kind=kint) :: i, idx0
720 
721  slave = contact%slave(lslave)
722  fcoeff = contact%fcoeff
723  ! clear contact force in former contacted element
724  nrlforce = contact%states(lslave)%multiplier(1)
725  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)+nrlforce*odirec
726  nn = size( contact%master(omaster)%nodes )
727  etype = contact%master(omaster)%etype
728  call getshapefunc( etype, opos(:), shapefunc )
729  do j=1,nn
730  iss = contact%master(omaster)%nodes(j)
731  ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)-nrlforce*shapefunc(j)*odirec
732  idx0 = 3*(iss-1)
733  do i=1,3
734  !$omp atomic
735  b(idx0+i) = b(idx0+i)-nrlforce*shapefunc(j)*odirec(i)
736  enddo
737  enddo
738  if( fcoeff/=0.d0 ) then
739  do j=1,nn
740  iss = contact%master(omaster)%nodes(j)
741  elemcrd(:,j) = currpos(3*iss-2:3*iss)
742  enddo
743  call dispincrematrix( opos(:), etype, nn, elemcrd, tangent, &
744  metric, dispmat )
745  fric(1:2) = contact%states(lslave)%multiplier(2:3)
746  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
747  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)+f3(1:3)
748  do j=1,nn
749  iss = contact%master(omaster)%nodes(j)
750  ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)+f3(3*j+1:3*j+3)
751  idx0 = 3*(iss-1)
752  do i=1,3
753  !$omp atomic
754  b(idx0+i) = b(idx0+i)+f3(3*j+i)
755  enddo
756  enddo
757  endif
758 
759  ! reset contact force in new contacting element
760  master = contact%states(lslave)%surface
761  nn = size( contact%master(master)%nodes )
762  etype = contact%master(master)%etype
763  call getshapefunc( etype, contact%states(lslave)%lpos(:), shapefunc )
764  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-nrlforce*contact%states(lslave)%direction
765  do j=1,nn
766  iss = contact%master(master)%nodes(j)
767  ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)+nrlforce* &
768  ! shapefunc(j)*contact%states(lslave)%direction
769  idx0 = 3*(iss-1)
770  do i=1,3
771  !$omp atomic
772  b(idx0+i) = b(idx0+i)+nrlforce* &
773  shapefunc(j)*contact%states(lslave)%direction(i)
774  enddo
775  enddo
776  if( fcoeff/=0.d0 ) then
777  do j=1,nn
778  iss = contact%master(master)%nodes(j)
779  elemcrd(:,j) = currpos(3*iss-2:3*iss)
780  enddo
781  call dispincrematrix( contact%states(lslave)%lpos, etype, nn, elemcrd, tangent, &
782  metric, dispmat )
783  fric(1:2) = contact%states(lslave)%multiplier(2:3)
784  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
785  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-f3(1:3)
786  do j=1,nn
787  iss = contact%master(master)%nodes(j)
788  ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)-f3(3*j+1:3*j+3)
789  idx0 = 3*(iss-1)
790  do i=1,3
791  !$omp atomic
792  b(idx0+i) = b(idx0+i)-f3(3*j+i)
793  enddo
794  enddo
795  endif
796 
797  end subroutine reset_contact_force
798 
802  subroutine calcu_contact_force0( contact, coord, disp, ddisp, fcoeff, mu, &
803  mut, B )
804  type( tcontact ), intent(inout) :: contact
805  real(kind=kreal), intent(in) :: coord(:)
806  real(kind=kreal), intent(in) :: disp(:)
807  real(kind=kreal), intent(in) :: ddisp(:)
808  real(kind=kreal), intent(in) :: fcoeff
809  real(kind=kreal), intent(in) :: mu, mut
810  real(kind=kreal), intent(inout) :: b(:)
811 
812  integer(kind=kint) :: slave, etype, master
813  integer(kind=kint) :: nn, i, j, iSS
814  real(kind=kreal) :: nrlforce, elemdisp(3,l_max_elem_node), tangent(3,2)
815  real(kind=kreal) :: dg(3), elemg(3), elemcrd(3, l_max_elem_node )
816  real(kind=kreal) :: dgn, dxi(2), dxy(2), shapefunc(l_max_surface_node)
817  real(kind=kreal) :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
818  real(kind=kreal) :: fric(2), f3(3*l_max_elem_node+3), edisp(3*l_max_elem_node+3)
819 
820  do i= 1, size(contact%slave)
821  if( contact%states(i)%state==contactfree ) cycle ! not in contact
822  slave = contact%slave(i)
823  edisp(1:3) = ddisp(3*slave-2:3*slave)
824  master = contact%states(i)%surface
825 
826  nn = size( contact%master(master)%nodes )
827  etype = contact%master(master)%etype
828  do j=1,nn
829  iss = contact%master(master)%nodes(j)
830  elemdisp(:,j) = ddisp(3*iss-2:3*iss)
831  edisp(3*j+1:3*j+3) = ddisp(3*iss-2:3*iss)
832  elemcrd(:,j) = coord(3*iss-2:3*iss)+disp(3*iss-2:3*iss)
833  enddo
834  call getshapefunc( etype, contact%states(i)%lpos(:), shapefunc )
835 
836  ! normal component
837  elemg = 0.d0
838  do j=1,nn
839  elemg(:) = elemg(:)+shapefunc(j)*elemdisp(:,j)
840  enddo
841  dg(:) = ddisp(3*slave-2:3*slave) - elemg(:)
842  dgn = dot_product( contact%states(i)%direction, dg )
843  nrlforce = contact%states(i)%multiplier(1)- mu*(contact%states(i)%wkdist-dgn)
844  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-nrlforce*contact%states(i)%direction
845 
846  do j=1,nn
847  iss = contact%master(master)%nodes(j)
848  b(3*iss-2:3*iss) = b(3*iss-2:3*iss)+nrlforce* &
849  shapefunc(j)*contact%states(i)%direction
850  enddo
851 
852  if( fcoeff==0.d0 ) cycle
853  ! tangent component
854  call dispincrematrix( contact%states(i)%lpos, etype, nn, elemcrd, tangent, &
855  metric, dispmat )
856  dxi(1) = dot_product( dispmat(1,1:nn*3+3), edisp(1:nn*3+3) )
857  dxi(2) = dot_product( dispmat(2,1:nn*3+3), edisp(1:nn*3+3) )
858  dxy(:) = matmul( metric, dxi )
859  fric(1:2) = contact%states(i)%multiplier(2:3) + mut*dxy(1:2)
860  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
861 
862  if( contact%states(i)%state==contactslip ) then
863  dgn = dsqrt( f3(1)*f3(1)+f3(2)*f3(2)+f3(3)*f3(3) )
864  f3(:) = f3(:)*fcoeff*contact%states(i)%multiplier(1)/dgn
865  endif
866  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-f3(1:3)
867  do j=1,nn
868  iss = contact%master(master)%nodes(j)
869  b(3*iss-2:3*iss) = b(3*iss-2:3*iss)-f3(3*j+1:3*j+3)
870  enddo
871  enddo
872  end subroutine calcu_contact_force0
873 
874 
877  subroutine update_contact_multiplier( contact, coord, disp, ddisp, fcoeff, mu, mut, &
878  gnt, ctchanged )
879  type( tcontact ), intent(inout) :: contact
880  real(kind=kreal), intent(in) :: coord(:)
881  real(kind=kreal), intent(in) :: disp(:)
882  real(kind=kreal), intent(in) :: ddisp(:)
883  real(kind=kreal), intent(in) :: fcoeff
884  real(kind=kreal), intent(in) :: mu, mut
885  real(kind=kreal), intent(out) :: gnt(2)
886  logical, intent(inout) :: ctchanged
887 
888  integer(kind=kint) :: slave, etype, master
889  integer(kind=kint) :: nn, i, j, iss, cnt
890  real(kind=kreal) :: elemdisp(3,l_max_elem_node), tangent(3,2)
891  real(kind=kreal) :: dg(3), elemg(3), elemcrd(3, l_max_elem_node )
892  real(kind=kreal) :: dgn, dxi(2), dxy(2), shapefunc(l_max_surface_node)
893  real(kind=kreal) :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
894  real(kind=kreal) :: lgnt(2), fric(2), f3(3*l_max_elem_node+3), edisp(3*l_max_elem_node+3)
895 
896  cnt =0; lgnt(:)=0.d0
897  do i= 1, size(contact%slave)
898  if( contact%states(i)%state==contactfree ) cycle ! not in contact
899  slave = contact%slave(i)
900  edisp(1:3) = ddisp(3*slave-2:3*slave)
901  master = contact%states(i)%surface
902 
903  nn = size( contact%master(master)%nodes )
904  etype = contact%master(master)%etype
905  do j=1,nn
906  iss = contact%master(master)%nodes(j)
907  elemdisp(:,j) = ddisp(3*iss-2:3*iss)
908  edisp(3*j+1:3*j+3) = ddisp(3*iss-2:3*iss)
909  elemcrd(:,j) = coord(3*iss-2:3*iss)+disp(3*iss-2:3*iss)
910  enddo
911  call getshapefunc( etype, contact%states(i)%lpos(:), shapefunc )
912 
913  ! normal component
914  elemg = 0.d0
915  do j=1,nn
916  elemg(:) = elemg(:)+shapefunc(j)*elemdisp(:,j)
917  enddo
918  dg(:) = ddisp(3*slave-2:3*slave) - elemg(:)
919  dgn = dot_product( contact%states(i)%direction, dg )
920  contact%states(i)%wkdist = contact%states(i)%wkdist-dgn
921  contact%states(i)%multiplier(1) = contact%states(i)%multiplier(1)- mu*contact%states(i)%wkdist
922  contact%states(i)%distance = contact%states(i)%distance - dgn
923  cnt = cnt+1
924  lgnt(1) = lgnt(1)- contact%states(i)%wkdist
925 
926  if( fcoeff==0.d0 ) cycle
927  ! tangent component
928  call dispincrematrix( contact%states(i)%lpos, etype, nn, elemcrd, tangent, &
929  metric, dispmat )
930  dxi(1) = dot_product( dispmat(1,1:nn*3+3), edisp(1:nn*3+3) )
931  dxi(2) = dot_product( dispmat(2,1:nn*3+3), edisp(1:nn*3+3) )
932  dxy(:) = matmul( metric, dxi )
933  fric(1:2) = contact%states(i)%multiplier(2:3) + mut*dxy(1:2)
934  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
935  dgn = dsqrt( f3(1)*f3(1)+f3(2)*f3(2)+f3(3)*f3(3) )
936  if( contact%states(i)%multiplier(1)>0.d0 ) then
937  if( dgn > fcoeff*contact%states(i)%multiplier(1) ) then
938  if( contact%states(i)%state==contactstick ) then
939  ctchanged= .true.
940  print *, "Node", slave, "to slip state",dgn, fcoeff*contact%states(i)%multiplier(1)
941  endif
942  contact%states(i)%state = contactslip
943  fric(:) = fric(:)*fcoeff*contact%states(i)%multiplier(1)/dgn
944  else
945  if( contact%states(i)%state==contactslip ) then
946  ctchanged= .true.
947  print *, "Node", slave, "to stick state",dgn, fcoeff*contact%states(i)%multiplier(1)
948  endif
949  contact%states(i)%state = contactstick
950  endif
951  endif
952  contact%states(i)%multiplier(2:3) = fric(:)
953  dxy(:) = matmul( dg, tangent )
954  lgnt(2) = lgnt(2)+dsqrt( dxy(1)*dxy(1)+dxy(2)*dxy(2) )
955  enddo
956  if(cnt>0) lgnt(:) = lgnt(:)/cnt
957  gnt = gnt + lgnt
958  end subroutine update_contact_multiplier
959 
961  subroutine ass_contact_force( contact, coord, disp, B )
962  type( tcontact ), intent(in) :: contact
963  real(kind=kreal), intent(in) :: coord(:)
964  real(kind=kreal), intent(in) :: disp(:)
965  real(kind=kreal), intent(inout) :: b(:)
966 
967  integer(kind=kint) :: slave, etype, master
968  integer(kind=kint) :: nn, i, j, iss
969  real(kind=kreal) :: fcoeff, nrlforce, tangent(3,2)
970  real(kind=kreal) :: elemcrd(3, l_max_elem_node )
971  real(kind=kreal) :: shapefunc(l_max_surface_node)
972  real(kind=kreal) :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
973  real(kind=kreal) :: fric(2), f3(l_max_elem_node*3+3)
974  fcoeff = contact%fcoeff
975  do i= 1, size(contact%slave)
976  if( contact%states(i)%state==contactfree ) cycle ! not in contact
977  slave = contact%slave(i)
978  master = contact%states(i)%surface
979 
980  nn = size( contact%master(master)%nodes )
981  etype = contact%master(master)%etype
982  call getshapefunc( etype, contact%states(i)%lpos(:), shapefunc )
983 
984  ! normal component
985  nrlforce = contact%states(i)%multiplier(1)
986  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-nrlforce*contact%states(i)%direction
987 
988  do j=1,nn
989  iss = contact%master(master)%nodes(j)
990  b(3*iss-2:3*iss) = b(3*iss-2:3*iss)+nrlforce* &
991  shapefunc(j)*contact%states(i)%direction
992  enddo
993 
994  if( fcoeff==0.d0 ) cycle
995  ! tangent component
996  do j=1,nn
997  iss = contact%master(master)%nodes(j)
998  elemcrd(:,j) = coord(3*iss-2:3*iss)+disp(3*iss-2:3*iss)
999  enddo
1000  call dispincrematrix( contact%states(i)%lpos, etype, nn, elemcrd, tangent, &
1001  metric, dispmat )
1002 
1003  fric(1:2) = contact%states(i)%multiplier(2:3)
1004  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
1005  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-f3(1:3)
1006  do j=1,nn
1007  iss = contact%master(master)%nodes(j)
1008  b(3*iss-2:3*iss) = b(3*iss-2:3*iss)-f3(3*j+1:3*j+3)
1009  enddo
1010  enddo
1011 
1012  end subroutine ass_contact_force
1013 
1015  subroutine set_contact_state_vector( contact, dt, relvel_vec, state_vec )
1016  type( tcontact ), intent(in) :: contact
1017  real(kind=kreal), intent(in) :: dt
1018  real(kind=kreal), intent(inout) :: relvel_vec(:)
1019  real(kind=kreal), intent(inout) :: state_vec(:)
1020 
1021  integer(kind=kint) :: i, slave
1022 
1023  do i= 1, size(contact%slave)
1024  slave = contact%slave(i)
1025  if( state_vec(slave) < 0.1d0 .or. contact%states(i)%state > 0 ) &
1026  & state_vec(slave) = dble(contact%states(i)%state)
1027 
1028  if( contact%states(i)%state==contactfree ) cycle ! not in contact
1029  if( dt < 1.d-16 ) cycle ! too small delta t
1030  relvel_vec(3*slave-2:3*slave) = contact%states(i)%reldisp(1:3)/dt
1031  enddo
1032 
1033  end subroutine set_contact_state_vector
1034 
1035  subroutine update_contact_tangentforce( contact )
1036  type( tcontact ), intent(inout) :: contact
1037 
1038  integer(kind=kint) :: i
1039 
1040  do i= 1, size(contact%slave)
1041  if( contact%states(i)%state==contactfree ) then
1042  contact%states(i)%tangentForce(1:3) = 0.d0
1043  contact%states(i)%tangentForce_trial(1:3) = 0.d0
1044  contact%states(i)%tangentForce_final(1:3) = 0.d0
1045  else
1046  contact%states(i)%tangentForce(1:3) = contact%states(i)%tangentForce_final(1:3)
1047  end if
1048  contact%states(i)%tangentForce1(1:3) = contact%states(i)%tangentForce(1:3)
1049  enddo
1050  end subroutine update_contact_tangentforce
1051 
1053  subroutine track_contact_position_exp( nslave, contact, currpos, currdisp, infoCTChange, nodeID, elemID )
1054  integer, intent(in) :: nslave
1055  type( tcontact ), intent(inout) :: contact
1056  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
1057  real(kind=kreal), intent(in) :: currpos(:)
1058  real(kind=kreal), intent(in) :: currdisp(:)
1059  integer(kind=kint), intent(in) :: nodeID(:)
1060  integer(kind=kint), intent(in) :: elemID(:)
1061 
1062  integer(kind=kint) :: slave, sid0, sid, etype
1063  integer(kind=kint) :: nn, i, j, iSS
1064  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node ), elem0(3, l_max_elem_node )
1065  logical :: isin
1066  real(kind=kreal) :: opos(2), odirec(3)
1067  integer(kind=kint) :: bktID, nCand, idm
1068  integer(kind=kint), allocatable :: indexCand(:)
1069 
1070  sid = 0
1071 
1072  slave = contact%slave(nslave)
1073  coord(:) = currpos(3*slave-2:3*slave)
1075  sid0 = contact%states(nslave)%surface
1076  opos = contact%states(nslave)%lpos
1077  odirec = contact%states(nslave)%direction
1078  etype = contact%master(sid0)%etype
1079  nn = getnumberofnodes( etype )
1080  do j=1,nn
1081  iss = contact%master(sid0)%nodes(j)
1082  elem(1:3,j)=currpos(3*iss-2:3*iss)
1083  elem0(1:3,j)=currpos(3*iss-2:3*iss)-currdisp(3*iss-2:3*iss)
1084  enddo
1085  call project_point2element( coord,etype,nn,elem,contact%master(sid0)%reflen,contact%states(nslave), &
1086  isin,distclr_nocheck,contact%states(nslave)%lpos,clr_same_elem )
1087  if( .not. isin ) then
1088  do i=1, contact%master(sid0)%n_neighbor
1089  sid = contact%master(sid0)%neighbor(i)
1090  etype = contact%master(sid)%etype
1091  nn = getnumberofnodes( etype )
1092  do j=1,nn
1093  iss = contact%master(sid)%nodes(j)
1094  elem(1:3,j)=currpos(3*iss-2:3*iss)
1095  enddo
1096  call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
1097  isin,distclr_nocheck,localclr=clearance )
1098  if( isin ) then
1099  contact%states(nslave)%surface = sid
1100  exit
1101  endif
1102  enddo
1103  endif
1104 
1105  if( .not. isin ) then ! such case is considered to rarely or never occur
1106  write(*,*) 'Warning: contact moved beyond neighbor elements'
1107  ! get master candidates from bucketDB
1108  bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
1109  ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
1110  if (ncand > 0) then
1111  allocate(indexcand(ncand))
1112  call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
1113  do idm= 1, ncand
1114  sid = indexcand(idm)
1115  if( sid==sid0 ) cycle
1116  if( any(sid==contact%master(sid0)%neighbor(:)) ) cycle
1117  if (.not. is_in_surface_box( contact%master(sid), coord(1:3), box_exp_rate )) cycle
1118  etype = contact%master(sid)%etype
1119  nn = size( contact%master(sid)%nodes )
1120  do j=1,nn
1121  iss = contact%master(sid)%nodes(j)
1122  elem(1:3,j)=currpos(3*iss-2:3*iss)
1123  enddo
1124  call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
1125  isin,distclr_nocheck,localclr=clearance )
1126  if( isin ) then
1127  contact%states(nslave)%surface = sid
1128  exit
1129  endif
1130  enddo
1131  deallocate(indexcand)
1132  endif
1133  endif
1134 
1135  if( isin ) then
1136  if( contact%states(nslave)%surface==sid0 ) then
1137  if(any(dabs(contact%states(nslave)%lpos(:)-opos(:)) >= clr_difflpos)) then
1138  !$omp atomic
1139  infoctchange%contact2difflpos = infoctchange%contact2difflpos + 1
1140  endif
1141  else
1142  write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3)') "Node",nodeid(slave)," move to contact with", &
1143  elemid(contact%master(sid)%eid), " with distance ", &
1144  contact%states(nslave)%distance," at ",contact%states(nslave)%lpos(:)
1145  !$omp atomic
1146  infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
1147  endif
1148  iss = isinsideelement( etype, contact%states(nslave)%lpos, clr_cal_norm )
1149  if( iss>0 ) then
1150  call cal_node_normal( contact%states(nslave)%surface, iss, contact%master, currpos, &
1151  contact%states(nslave)%lpos, contact%states(nslave)%direction(:) )
1152  endif
1153  else if( .not. isin ) then
1154  write(*,'(A,i10,A)') "Node",nodeid(slave)," move out of contact"
1155  contact%states(nslave)%state = contactfree
1156  contact%states(nslave)%multiplier(:) = 0.d0
1157  endif
1158 
1159  end subroutine track_contact_position_exp
1160 
1164  subroutine scan_contact_state_exp( contact, currpos, currdisp, infoCTChange, &
1165  nodeID, elemID, is_init, active )
1166  type( tcontact ), intent(inout) :: contact
1167  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
1168  real(kind=kreal), intent(in) :: currpos(:)
1169  real(kind=kreal), intent(in) :: currdisp(:)
1170  integer(kind=kint), intent(in) :: nodeID(:)
1171  integer(kind=kint), intent(in) :: elemID(:)
1172  logical, intent(in) :: is_init
1173  logical, intent(out) :: active
1174 
1175  real(kind=kreal) :: distclr
1176  integer(kind=kint) :: slave, id, etype
1177  integer(kind=kint) :: nn, i, j, iSS, nactive
1178  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node )
1179  real(kind=kreal) :: nlforce
1180  logical :: isin
1181  integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
1182  !
1183  integer, pointer :: indexMaster(:),indexCand(:)
1184  integer :: nMaster,idm,nMasterMax,bktID,nCand
1185  logical :: is_cand
1186 
1187  if( is_init ) then
1188  distclr = distclr_init
1189  else
1190  distclr = distclr_free
1191  endif
1192 
1193  allocate(contact_surf(size(nodeid)))
1194  allocate(states_prev(size(contact%slave)))
1195  contact_surf(:) = size(elemid)+1
1196  do i = 1, size(contact%slave)
1197  states_prev(i) = contact%states(i)%state
1198  enddo
1199 
1200  call update_surface_box_info( contact%master, currpos )
1201  call update_surface_bucket_info( contact%master, contact%master_bktDB )
1202 
1203  !$omp parallel do &
1204  !$omp& default(none) &
1205  !$omp& private(i,slave,id,nlforce,coord,indexMaster,nMaster,nn,j,iSS,elem,is_cand,idm,etype,isin, &
1206  !$omp& bktID,nCand,indexCand) &
1207  !$omp& firstprivate(nMasterMax) &
1208  !$omp& shared(contact,infoCTChange,currpos,currdisp,nodeID,elemID,distclr,contact_surf) &
1209  !$omp& reduction(.or.:active) &
1210  !$omp& schedule(dynamic,1)
1211  do i= 1, size(contact%slave)
1212  slave = contact%slave(i)
1213  if( contact%states(i)%state==contactstick .or. contact%states(i)%state==contactslip ) then
1214  call track_contact_position_exp( i, contact, currpos, currdisp, infoctchange, nodeid, elemid )
1215  if( contact%states(i)%state /= contactfree ) then
1216  contact_surf(contact%slave(i)) = -contact%states(i)%surface
1217  endif
1218  else if( contact%states(i)%state==contactfree ) then
1219  coord(:) = currpos(3*slave-2:3*slave)
1220 
1221  ! get master candidates from bucketDB
1222  bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
1223  ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
1224  if (ncand == 0) cycle
1225  allocate(indexcand(ncand))
1226  call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
1227 
1228  nmastermax = ncand
1229  allocate(indexmaster(nmastermax))
1230  nmaster = 0
1231 
1232  ! narrow down candidates
1233  do idm= 1, ncand
1234  id = indexcand(idm)
1235  if (.not. is_in_surface_box( contact%master(id), coord(1:3), box_exp_rate )) cycle
1236  nmaster = nmaster + 1
1237  indexmaster(nmaster) = id
1238  enddo
1239  deallocate(indexcand)
1240 
1241  if(nmaster == 0) then
1242  deallocate(indexmaster)
1243  cycle
1244  endif
1245 
1246  do idm = 1,nmaster
1247  id = indexmaster(idm)
1248  etype = contact%master(id)%etype
1249  nn = size( contact%master(id)%nodes )
1250  do j=1,nn
1251  iss = contact%master(id)%nodes(j)
1252  elem(1:3,j)=currpos(3*iss-2:3*iss)
1253  enddo
1254  call project_point2element( coord,etype,nn,elem,contact%master(id)%reflen,contact%states(i), &
1255  isin,distclr,localclr=clearance )
1256  if( .not. isin ) cycle
1257  contact%states(i)%surface = id
1258  contact%states(i)%multiplier(:) = 0.d0
1259  iss = isinsideelement( etype, contact%states(i)%lpos, clr_cal_norm )
1260  if( iss>0 ) &
1261  call cal_node_normal( id, iss, contact%master, currpos, contact%states(i)%lpos, &
1262  contact%states(i)%direction(:) )
1263  contact_surf(contact%slave(i)) = id
1264  write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3,A,3f7.3)') "Node",nodeid(slave)," contact with element", &
1265  elemid(contact%master(id)%eid), &
1266  " with distance ", contact%states(i)%distance," at ",contact%states(i)%lpos(:), &
1267  " along direction ", contact%states(i)%direction
1268  exit
1269  enddo
1270  deallocate(indexmaster)
1271  endif
1272  enddo
1273  !$omp end parallel do
1274 
1275  call fstr_contact_comm_allreduce_i(contact%comm, contact_surf, hecmw_min)
1276  nactive = 0
1277  do i = 1, size(contact%slave)
1278  if (contact%states(i)%state /= contactfree) then ! any slave in contact
1279  if (abs(contact_surf(contact%slave(i))) /= contact%states(i)%surface) then ! that is in contact with other surface
1280  contact%states(i)%state = contactfree ! should be freed
1281  write(*,'(A,i10,A,i6,A,i6,A)') "Node",nodeid(contact%slave(i)), &
1282  " in rank",hecmw_comm_get_rank()," freed due to duplication"
1283  else
1284  nactive = nactive + 1
1285  endif
1286  endif
1287  if (states_prev(i) == contactfree .and. contact%states(i)%state /= contactfree) then
1288  infoctchange%free2contact = infoctchange%free2contact + 1
1289  elseif (states_prev(i) /= contactfree .and. contact%states(i)%state == contactfree) then
1290  infoctchange%contact2free = infoctchange%contact2free + 1
1291  endif
1292  enddo
1293  active = (nactive > 0)
1294  deallocate(contact_surf)
1295  deallocate(states_prev)
1296  end subroutine scan_contact_state_exp
1297 
1298 end module mcontactdef
This module provides bucket-search functionality It provides definition of bucket info and its access...
integer(kind=kint) function, public bucketdb_getbucketid(bktdb, x)
Get bucket ID that includes given point.
subroutine, public bucketdb_finalize(bktdb)
Finalizer.
integer(kind=kint) function, public bucketdb_getnumcand(bktdb, bid)
Get number of candidates within neighboring buckets of a given bucket Bucket ID has to be obtained wi...
subroutine, public bucketdb_init(bktdb)
Initializer.
subroutine, public bucketdb_getcand(bktdb, bid, ncand, cand)
Get candidates within neighboring buckets of a given bucket Number of candidates has to be obtained w...
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getvertexcoord(fetype, cnode, localcoord)
Get the natural coord of a vertex node.
Definition: element.f90:1086
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:863
Definition: hecmw.f90:6
This module provide functions of contact stiffness calculation.
Definition: contact_lib.f90:6
integer, parameter contactfslid
Definition: contact_lib.f90:24
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.
subroutine update_tangentforce(etype, nn, elemt0, elemt, cstate)
This subroutine find the projection of a slave point onto master surface.
subroutine contact_state_init(cstate)
Initializer.
Definition: contact_lib.f90:48
integer, parameter contactslip
Definition: contact_lib.f90:18
integer, parameter contactfree
contact state definition
Definition: contact_lib.f90:16
subroutine dispincrematrix(pos, etype, nnode, ele, tangent, tensor, matrix)
This subroutine calculate the relation between global disp and displacement along natural coordinate ...
integer, parameter contactstick
Definition: contact_lib.f90:17
subroutine, public fstr_contact_comm_init(conComm, hecMESH, ndof, n_contact_dof, contact_dofs)
subroutine, public fstr_contact_comm_finalize(conComm)
subroutine, public fstr_contact_comm_allreduce_i(conComm, vec, op)
This module manage the data structure for contact calculation.
subroutine set_contact_state_vector(contact, dt, relvel_vec, state_vec)
This subroutine setup contact output nodal vectors.
real(kind=kreal), parameter distclr_init
dist clearance for initial scan
subroutine update_contact_multiplier(contact, coord, disp, ddisp, fcoeff, mu, mut, gnt, ctchanged)
This subroutine update lagrangian multiplier and the distance between contacting nodes.
logical function fstr_contact_check(contact, hecMESH)
Check the consistency with given mesh of contact defintiion.
subroutine clear_contact_state(contact)
Reset contact state all to free.
subroutine update_contact_tangentforce(contact)
subroutine track_contact_position_exp(nslave, contact, currpos, currdisp, infoCTChange, nodeID, elemID)
This subroutine tracks down next contact position after a finite slide.
subroutine scan_contact_state_exp(contact, currpos, currdisp, infoCTChange, nodeID, elemID, is_init, active)
This subroutine update contact states, which include.
subroutine fstr_write_contact(file, contact)
Write out contact definition.
subroutine fstr_contact_finalize(contact)
Finalizer.
logical function fstr_contact_init(contact, hecMESH, myrank)
Initializer of tContactState.
subroutine ass_contact_force(contact, coord, disp, B)
This subroutine assemble contact force into contacing nodes.
subroutine calcu_contact_force0(contact, coord, disp, ddisp, fcoeff, mu, mut, B)
This subroutine update contact condition as follows:
subroutine scan_contact_state(flag_ctAlgo, contact, currpos, currdisp, ndforce, infoCTChange, nodeID, elemID, is_init, active, mu, B)
This subroutine update contact states, which include.
This module manage surface elements in 3D It provides basic definition of surface elements (trianglar...
Definition: surf_ele.f90:8
subroutine initialize_surf(eid, etype, nsurf, surf)
Initializer.
Definition: surf_ele.f90:38
subroutine update_surface_reflen(surf, coord)
Compute reference length of surface elements.
Definition: surf_ele.f90:212
subroutine write_surf(file, surf)
Write out elemental surface.
Definition: surf_ele.f90:68
subroutine update_surface_box_info(surf, currpos)
Update info of cubic box including surface elements.
Definition: surf_ele.f90:232
subroutine find_surface_neighbor(surf, bktDB)
Find neighboring surface elements.
Definition: surf_ele.f90:80
subroutine update_surface_bucket_info(surf, bktDB)
Update bucket info for searching surface elements.
Definition: surf_ele.f90:272
subroutine finalize_surf(surf)
Memeory management subroutine.
Definition: surf_ele.f90:61
logical function is_in_surface_box(surf, x0, exp_rate)
Check if given point is within cubic box including surface element.
Definition: surf_ele.f90:257
This structure records contact status.
Definition: contact_lib.f90:27
Structure to includes all info needed by contact calculation.
Structure to define surface group.
Definition: surf_ele.f90:19