FrontISTR  5.7.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 !-------------------------------------------------------------------------------
12 module mcontactdef
13 
14  use hecmw
15  use msurfelement
16  use m_contact_lib
17  use m_hecmw_contact_comm
18  use bucket_search
19  use mcontactparam
20 
21  implicit none
22 
23  include 'fstr_ctrl_util_f.inc'
24 
26  type tcontact
27  ! following contact definition
28  character(len=HECMW_NAME_LEN) :: name
29  integer :: ctype
30  integer :: group
31  character(len=HECMW_NAME_LEN) :: pair_name
32  integer :: surf_id1, surf_id2
33  integer :: surf_id1_sgrp
34  type(tsurfelement), pointer :: master(:)=>null()
35  integer, pointer :: slave(:)=>null()
36  real(kind=kreal) :: fcoeff
37  real(kind=kreal) :: tpenalty
38 
39  ! following algorithm
40  ! -1: not initialized
41  ! 1: TIED-Just rigidly fixed the two surfaces
42  ! 2: GLUED-Distance between the two surfaces to zero and glue them
43  ! 3: SSLID-Small sliding contact( no position but with contact state change)
44  ! 4: FSLID-Finite sliding contact (both changes in contact state and position possible)
45  integer :: algtype
46 
47  logical :: mpced
48  logical :: symmetric
49 
50  ! following contact state
51  type(tcontactstate), pointer :: states(:)=>null()
52 
53  type(hecmwst_contact_comm) :: comm
54  type(bucketdb) :: master_bktdb
55 
56  type(tcontactparam), pointer :: cparam=>null()
57  end type tcontact
58 
59  type fstr_info_contactchange
60  logical :: active
61  integer(kind=kint) :: contact2free
62  integer(kind=kint) :: contact2neighbor
63  integer(kind=kint) :: contact2difflpos
64  integer(kind=kint) :: free2contact
65  integer(kind=kint) :: contactnode_previous
66  integer(kind=kint) :: contactnode_current
67  end type fstr_info_contactchange
68 
69  private :: is_mpc_available
70  private :: is_active_contact
71  private :: cal_node_normal
72  private :: track_contact_position
73  private :: reset_contact_force
74 
75 contains
76 
78  logical function is_mpc_available( contact )
79  type(tcontact), intent(in) :: contact
80  is_mpc_available = .true.
81  if( contact%fcoeff/=0.d0 ) is_mpc_available = .false.
82  end function
83 
85  subroutine fstr_write_contact( file, contact )
86  integer(kind=kint), intent(in) :: file
87  type(tcontact), intent(in) :: contact
88  integer :: i
89  write(file,*) "CONTACT:", contact%ctype,contact%group,trim(contact%pair_name),contact%fcoeff
90  write(file,*) "---Slave----"
91  write(file,*) "num.slave",size(contact%slave)
92  if( associated(contact%slave) ) then
93  do i=1,size(contact%slave)
94  write(file, *) contact%slave(i)
95  enddo
96  endif
97  write(file,*) "----master---"
98  write(file,*) "num.master",size(contact%master)
99  if( associated(contact%master) ) then
100  do i=1,size(contact%master)
101  call write_surf( file, contact%master(i) )
102  enddo
103  endif
104  end subroutine
105 
107  subroutine fstr_contact_finalize( contact )
108  type(tcontact), intent(inout) :: contact
109  integer :: i
110  if( associated( contact%slave ) ) deallocate(contact%slave)
111  if( associated( contact%master ) ) then
112  do i=1,size( contact%master )
113  call finalize_surf( contact%master(i) )
114  enddo
115  deallocate(contact%master)
116  endif
117  if( associated(contact%states) ) deallocate(contact%states)
118  call hecmw_contact_comm_finalize(contact%comm)
119  call bucketdb_finalize( contact%master_bktDB )
120  end subroutine
121 
123  logical function fstr_contact_check( contact, hecMESH )
124  type(tcontact), intent(inout) :: contact
125  type(hecmwst_local_mesh), pointer :: hecmesh
126 
127  integer :: i
128  logical :: isfind
129 
130  fstr_contact_check = .false.
131 
132  ! if contact pair exist?
133  isfind = .false.
134  do i=1,hecmesh%contact_pair%n_pair
135  if( hecmesh%contact_pair%name(i) == contact%pair_name ) then
136  contact%ctype = hecmesh%contact_pair%type(i)
137  contact%surf_id1 = hecmesh%contact_pair%slave_grp_id(i)
138  contact%surf_id2 = hecmesh%contact_pair%master_grp_id(i)
139  contact%surf_id1_sgrp = hecmesh%contact_pair%slave_orisgrp_id(i)
140  isfind = .true.
141  endif
142  enddo
143  if( .not. isfind ) return;
144  if( contact%fcoeff<=0.d0 ) contact%fcoeff=0.d0
145  if( contact%ctype < 1 .and. contact%ctype > 3 ) return
146  if( contact%group<=0 ) return
147 
148  fstr_contact_check = .true.
149  end function
150 
152  logical function fstr_contact_init( contact, hecMESH, cparam, myrank )
153  type(tcontact), intent(inout) :: contact
154  type(hecmwst_local_mesh), pointer :: hecmesh
155  type(tcontactparam), target :: cparam
156  integer(kint),intent(in),optional :: myrank
157 
158  integer :: i, j, is, ie, cgrp, nsurf, nslave, ic, ic_type, iss, nn, ii
159  integer :: count,nodeid
160 
161  fstr_contact_init = .false.
162 
163  contact%cparam => cparam
164 
165  ! master surface
166  cgrp = contact%surf_id2
167  if( cgrp<=0 ) return
168  is= hecmesh%surf_group%grp_index(cgrp-1) + 1
169  ie= hecmesh%surf_group%grp_index(cgrp )
170 
171  if(present(myrank)) then
172  ! PARA_CONTACT
173  count = 0
174  do i=is,ie
175  ic = hecmesh%surf_group%grp_item(2*i-1)
176  if(hecmesh%elem_ID(ic*2) /= myrank) cycle
177  count = count + 1
178  enddo
179  allocate( contact%master(count) )
180  count = 0
181  do i=is,ie
182  ic = hecmesh%surf_group%grp_item(2*i-1)
183  if(hecmesh%elem_ID(ic*2) /= myrank) cycle
184  count = count + 1
185  nsurf = hecmesh%surf_group%grp_item(2*i)
186  ic_type = hecmesh%elem_type(ic)
187  call initialize_surf( ic, ic_type, nsurf, contact%master(count) )
188  iss = hecmesh%elem_node_index(ic-1)
189  do j=1, size( contact%master(count)%nodes )
190  nn = contact%master(count)%nodes(j)
191  contact%master(count)%nodes(j) = hecmesh%elem_node_item( iss+nn )
192  enddo
193  enddo
194 
195  else
196  ! not PARA_CONTACT
197  allocate( contact%master(ie-is+1) )
198  do i=is,ie
199  ic = hecmesh%surf_group%grp_item(2*i-1)
200  nsurf = hecmesh%surf_group%grp_item(2*i)
201  ic_type = hecmesh%elem_type(ic)
202  call initialize_surf( ic, ic_type, nsurf, contact%master(i-is+1) )
203  iss = hecmesh%elem_node_index(ic-1)
204  do j=1, size( contact%master(i-is+1)%nodes )
205  nn = contact%master(i-is+1)%nodes(j)
206  contact%master(i-is+1)%nodes(j) = hecmesh%elem_node_item( iss+nn )
207  enddo
208  enddo
209 
210  endif
211 
212  call update_surface_reflen( contact%master, hecmesh%node )
213 
214  cgrp = contact%surf_id1
215  if( cgrp<=0 ) return
216  is= hecmesh%node_group%grp_index(cgrp-1) + 1
217  ie= hecmesh%node_group%grp_index(cgrp )
218  nslave = 0
219  do i=is,ie
220  nodeid = hecmesh%global_node_ID(hecmesh%node_group%grp_item(i))
221  if(present(myrank)) then
222  ! PARA_CONTACT
223  nslave = nslave + 1
224  else
225  ! not PARA_CONTACT
226  if( hecmesh%node_group%grp_item(i) <= hecmesh%nn_internal) then
227  nslave = nslave + 1
228  endif
229  endif
230  enddo
231  allocate( contact%slave(nslave) )
232  allocate( contact%states(nslave) )
233  ii = 0
234  do i=is,ie
235  if(.not.present(myrank)) then
236  ! not PARA_CONTACT
237  if( hecmesh%node_group%grp_item(i) > hecmesh%nn_internal) cycle
238  endif
239  ii = ii + 1
240  contact%slave(ii) = hecmesh%node_group%grp_item(i)
241  contact%states(ii)%state = -1
242  contact%states(ii)%multiplier(:) = 0.d0
243  contact%states(ii)%tangentForce(:) = 0.d0
244  contact%states(ii)%tangentForce1(:) = 0.d0
245  contact%states(ii)%tangentForce_trial(:) = 0.d0
246  contact%states(ii)%tangentForce_final(:) = 0.d0
247  contact%states(ii)%reldisp(:) = 0.d0
248  enddo
249  ! endif
250 
251  ! contact state
252  ! allocate( contact%states(nslave) )
253  do i=1,nslave
254  call contact_state_init( contact%states(i) )
255  enddo
256 
257  ! neighborhood of surface group
258  call update_surface_box_info( contact%master, hecmesh%node )
259  call bucketdb_init( contact%master_bktDB )
260  call update_surface_bucket_info( contact%master, contact%master_bktDB )
261  call find_surface_neighbor( contact%master, contact%master_bktDB )
262 
263  ! initialize contact communication table
264  call hecmw_contact_comm_init( contact%comm, hecmesh, 1, nslave, contact%slave )
265 
266  contact%symmetric = .true.
267  fstr_contact_init = .true.
268  end function
269 
271  logical function fstr_embed_init( embed, hecMESH, cparam, myrank )
272  type(tcontact), intent(inout) :: embed
273  type(hecmwst_local_mesh), pointer :: hecmesh
274  type(tcontactparam), target :: cparam
275  integer(kint),intent(in),optional :: myrank
276 
277  integer :: i, j, is, ie, cgrp, nsurf, nslave, ic, ic_type, iss, nn, ii
278  integer :: count,nodeid
279 
280  fstr_embed_init = .false.
281 
282  embed%cparam => cparam
283 
284  ! master surface
285  cgrp = embed%surf_id2
286  if( cgrp<=0 ) return
287  is= hecmesh%elem_group%grp_index(cgrp-1) + 1
288  ie= hecmesh%elem_group%grp_index(cgrp )
289 
290  if(present(myrank)) then
291  ! PARA_CONTACT
292  count = 0
293  do i=is,ie
294  ic = hecmesh%elem_group%grp_item(i)
295  if(hecmesh%elem_ID(ic*2) /= myrank) cycle
296  count = count + 1
297  enddo
298  allocate( embed%master(count) )
299  count = 0
300  do i=is,ie
301  ic = hecmesh%elem_group%grp_item(i)
302  if(hecmesh%elem_ID(ic*2) /= myrank) cycle
303  count = count + 1
304  ic_type = hecmesh%elem_type(ic)
305  call initialize_surf( ic, ic_type, 0, embed%master(count) )
306  iss = hecmesh%elem_node_index(ic-1)
307  do j=1, size( embed%master(count)%nodes )
308  nn = embed%master(count)%nodes(j)
309  embed%master(count)%nodes(j) = hecmesh%elem_node_item( iss+nn )
310  enddo
311  enddo
312 
313  else
314  ! not PARA_CONTACT
315  allocate( embed%master(ie-is+1) )
316  do i=is,ie
317  ic = hecmesh%elem_group%grp_item(i)
318  ic_type = hecmesh%elem_type(ic)
319  call initialize_surf( ic, ic_type, 0, embed%master(i-is+1) )
320  iss = hecmesh%elem_node_index(ic-1)
321  do j=1, size( embed%master(i-is+1)%nodes )
322  nn = embed%master(i-is+1)%nodes(j)
323  embed%master(i-is+1)%nodes(j) = hecmesh%elem_node_item( iss+nn )
324  enddo
325  enddo
326 
327  endif
328 
329  !call update_surface_reflen( embed%master, hecMESH%node )
330 
331  ! slave surface
332  ! if( contact%ctype==1 ) then
333  cgrp = embed%surf_id1
334  if( cgrp<=0 ) return
335  is= hecmesh%node_group%grp_index(cgrp-1) + 1
336  ie= hecmesh%node_group%grp_index(cgrp )
337  nslave = 0
338  do i=is,ie
339  nodeid = hecmesh%global_node_ID(hecmesh%node_group%grp_item(i))
340  if(present(myrank)) then
341  ! PARA_CONTACT
342  nslave = nslave + 1
343  else
344  ! not PARA_CONTACT
345  if( hecmesh%node_group%grp_item(i) <= hecmesh%nn_internal) then
346  nslave = nslave + 1
347  endif
348  endif
349  enddo
350  allocate( embed%slave(nslave) )
351  allocate( embed%states(nslave) )
352  ii = 0
353  do i=is,ie
354  if(.not.present(myrank)) then
355  ! not PARA_CONTACT
356  if( hecmesh%node_group%grp_item(i) > hecmesh%nn_internal) cycle
357  endif
358  ii = ii + 1
359  embed%slave(ii) = hecmesh%node_group%grp_item(i)
360  embed%states(ii)%state = -1
361  embed%states(ii)%multiplier(:) = 0.d0
362  embed%states(ii)%tangentForce(:) = 0.d0
363  embed%states(ii)%tangentForce1(:) = 0.d0
364  embed%states(ii)%tangentForce_trial(:) = 0.d0
365  embed%states(ii)%tangentForce_final(:) = 0.d0
366  embed%states(ii)%reldisp(:) = 0.d0
367  enddo
368  ! endif
369 
370  ! contact state
371  ! allocate( contact%states(nslave) )
372  do i=1,nslave
373  call contact_state_init( embed%states(i) )
374  enddo
375 
376  ! neighborhood of surface group
377  call update_surface_box_info( embed%master, hecmesh%node )
378  call bucketdb_init( embed%master_bktDB )
379  call update_surface_bucket_info( embed%master, embed%master_bktDB )
380  call find_surface_neighbor( embed%master, embed%master_bktDB )
381 
382  ! initialize contact communication table
383  call hecmw_contact_comm_init( embed%comm, hecmesh, 1, nslave, embed%slave )
384 
385  embed%symmetric = .true.
386  fstr_embed_init = .true.
387  end function
388 
389 
391  subroutine clear_contact_state( contact )
392  type(tcontact), intent(inout) :: contact
393  integer :: i
394  if( .not. associated(contact%states) ) return
395  do i=1,size( contact%states )
396  contact%states(i)%state = -1
397  enddo
398  end subroutine
399 
401  logical function is_active_contact( acgrp, contact )
402  integer, intent(in) :: acgrp(:)
403  type(tcontact), intent(in) :: contact
404  if( any( acgrp==contact%group ) ) then
405  is_active_contact = .true.
406  else
407  is_active_contact = .false.
408  endif
409  end function
410 
414  subroutine scan_contact_state( flag_ctAlgo, contact, currpos, currdisp, ndforce, infoCTChange, &
415  nodeID, elemID, is_init, active, mu, B )
416  character(len=9), intent(in) :: flag_ctAlgo
417  type( tcontact ), intent(inout) :: contact
418  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
419  real(kind=kreal), intent(in) :: currpos(:)
420  real(kind=kreal), intent(in) :: currdisp(:)
421  real(kind=kreal), intent(in) :: ndforce(:)
422  integer(kind=kint), intent(in) :: nodeID(:)
423  integer(kind=kint), intent(in) :: elemID(:)
424  logical, intent(in) :: is_init
425  logical, intent(out) :: active
426  real(kind=kreal), intent(in) :: mu
427  real(kind=kreal), optional, target :: b(:)
428 
429  real(kind=kreal) :: distclr
430  integer(kind=kint) :: slave, id, etype
431  integer(kind=kint) :: nn, i, j, iSS, nactive
432  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node )
433  real(kind=kreal) :: nlforce, slforce(3)
434  logical :: isin
435  integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
436  !
437  integer, pointer :: indexMaster(:),indexCand(:)
438  integer :: nMaster,idm,nMasterMax,bktID,nCand
439  logical :: is_cand, is_present_B
440  real(kind=kreal), pointer :: bp(:)
441 
442  if( is_init ) then
443  distclr = contact%cparam%DISTCLR_INIT
444  else
445  distclr = contact%cparam%DISTCLR_FREE
446  if( contact%algtype == contacttied ) then
447  active = .false.
448  do i= 1, size(contact%slave)
449  if( contact%states(i)%state==contactstick ) then
450  active = .true.
451  exit
452  endif
453  enddo
454  endif
455  endif
456 
457  allocate(contact_surf(size(nodeid)))
458  allocate(states_prev(size(contact%slave)))
459  contact_surf(:) = huge(0)
460  do i = 1, size(contact%slave)
461  states_prev(i) = contact%states(i)%state
462  enddo
463 
464  call update_surface_box_info( contact%master, currpos )
465  call update_surface_bucket_info( contact%master, contact%master_bktDB )
466 
467  ! for gfortran-10: optional parameter seems not allowed within omp parallel
468  is_present_b = present(b)
469  if( is_present_b ) bp => b
470 
471  !$omp parallel do &
472  !$omp& default(none) &
473  !$omp& private(i,slave,slforce,id,nlforce,coord,indexMaster,nMaster,nn,j,iSS,elem,is_cand,idm,etype,isin, &
474  !$omp& bktID,nCand,indexCand) &
475  !$omp& firstprivate(nMasterMax,is_present_B) &
476  !$omp& shared(contact,ndforce,flag_ctAlgo,infoCTChange,currpos,currdisp,mu,nodeID,elemID,Bp,distclr,contact_surf,is_init) &
477  !$omp& reduction(.or.:active) &
478  !$omp& schedule(dynamic,1)
479  do i= 1, size(contact%slave)
480  slave = contact%slave(i)
481  if( contact%states(i)%state==contactstick .or. contact%states(i)%state==contactslip ) then
482  slforce(1:3)=ndforce(3*slave-2:3*slave)
483  id = contact%states(i)%surface
484  nlforce = contact%states(i)%multiplier(1)
485 
486  ! update direction of TIED contact
487  if( contact%algtype == contacttied ) then
488  call update_direction( i, contact, currpos )
489  if (.not.is_init) cycle
490  endif
491 
492  if( nlforce < contact%cparam%TENSILE_FORCE ) then
493  contact%states(i)%state = contactfree
494  contact%states(i)%multiplier(:) = 0.d0
495  write(*,'(A,i10,A,i10,A,e12.3)') "Node",nodeid(slave)," free from contact with element", &
496  elemid(contact%master(id)%eid), " with tensile force ", nlforce
497  cycle
498  endif
499  if( contact%algtype /= contactfslid .or. (.not. is_present_b) ) then ! small slide problem
500  contact_surf(contact%slave(i)) = -elemid(contact%master(id)%eid)
501  else
502  call track_contact_position( flag_ctalgo, i, contact, currpos, currdisp, mu, infoctchange, nodeid, elemid, bp )
503  if( contact%states(i)%state /= contactfree ) then
504  id = contact%states(i)%surface
505  contact_surf(contact%slave(i)) = -elemid(contact%master(id)%eid)
506  endif
507  endif
508 
509  else if( contact%states(i)%state==contactfree ) then
510  if( contact%algtype == contacttied .and. .not. is_init ) cycle
511  coord(:) = currpos(3*slave-2:3*slave)
512 
513  ! get master candidates from bucketDB
514  bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
515  ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
516  if (ncand == 0) cycle
517  allocate(indexcand(ncand))
518  call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
519 
520  nmastermax = ncand
521  allocate(indexmaster(nmastermax))
522  nmaster = 0
523 
524  ! narrow down candidates
525  do idm= 1, ncand
526  id = indexcand(idm)
527  if (.not. is_in_surface_box( contact%master(id), coord(1:3), contact%cparam%BOX_EXP_RATE )) cycle
528  nmaster = nmaster + 1
529  indexmaster(nmaster) = id
530  enddo
531  deallocate(indexcand)
532 
533  if(nmaster == 0) then
534  deallocate(indexmaster)
535  cycle
536  endif
537 
538  do idm = 1,nmaster
539  id = indexmaster(idm)
540  etype = contact%master(id)%etype
541  nn = size( contact%master(id)%nodes )
542  do j=1,nn
543  iss = contact%master(id)%nodes(j)
544  elem(1:3,j)=currpos(3*iss-2:3*iss)
545  enddo
546  call project_point2element( coord,etype,nn,elem,contact%master(id)%reflen,contact%states(i), &
547  isin,distclr,localclr=contact%cparam%CLEARANCE )
548  if( .not. isin ) cycle
549  contact%states(i)%surface = id
550  contact%states(i)%multiplier(:) = 0.d0
551  iss = isinsideelement( etype, contact%states(i)%lpos(1:2), contact%cparam%CLR_CAL_NORM )
552  if( iss>0 ) &
553  call cal_node_normal( id, iss, contact%master, currpos, contact%states(i)%lpos(1:2), &
554  contact%states(i)%direction(:) )
555  contact_surf(contact%slave(i)) = elemid(contact%master(id)%eid)
556  write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3,A,3f7.3,A,i6)') "Node",nodeid(slave)," contact with element", &
557  elemid(contact%master(id)%eid), &
558  " with distance ", contact%states(i)%distance," at ",contact%states(i)%lpos(1:2), &
559  " along direction ", contact%states(i)%direction," rank=",hecmw_comm_get_rank()
560  exit
561  enddo
562  deallocate(indexmaster)
563  endif
564  enddo
565  !$omp end parallel do
566 
567  if( contact%algtype == contacttied .and. .not. is_init ) then
568  deallocate(contact_surf)
569  deallocate(states_prev)
570  return
571  endif
572 
573  call hecmw_contact_comm_allreduce_i(contact%comm, contact_surf, hecmw_min)
574  nactive = 0
575  do i = 1, size(contact%slave)
576  if (contact%states(i)%state /= contactfree) then ! any slave in contact
577  id = contact%states(i)%surface
578  if (abs(contact_surf(contact%slave(i))) /= elemid(contact%master(id)%eid)) then ! that is in contact with other surface
579  contact%states(i)%state = contactfree ! should be freed
580  write(*,'(A,i10,A,i10,A,i6,A,i6,A)') "Node",nodeid(contact%slave(i))," contact with element", &
581  & elemid(contact%master(id)%eid), " in rank",hecmw_comm_get_rank()," freed due to duplication"
582  else
583  nactive = nactive + 1
584  endif
585  endif
586  if (states_prev(i) == contactfree .and. contact%states(i)%state /= contactfree) then
587  infoctchange%free2contact = infoctchange%free2contact + 1
588  elseif (states_prev(i) /= contactfree .and. contact%states(i)%state == contactfree) then
589  infoctchange%contact2free = infoctchange%contact2free + 1
590  endif
591  enddo
592  active = (nactive > 0)
593  deallocate(contact_surf)
594  deallocate(states_prev)
595  end subroutine scan_contact_state
596 
600  subroutine scan_embed_state( flag_ctAlgo, embed, currpos, currdisp, ndforce, infoCTChange, &
601  nodeID, elemID, is_init, active, mu, B )
602  character(len=9), intent(in) :: flag_ctAlgo
603  type( tcontact ), intent(inout) :: embed
604  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
605  real(kind=kreal), intent(in) :: currpos(:)
606  real(kind=kreal), intent(in) :: currdisp(:)
607  real(kind=kreal), intent(in) :: ndforce(:)
608  integer(kind=kint), intent(in) :: nodeID(:)
609  integer(kind=kint), intent(in) :: elemID(:)
610  logical, intent(in) :: is_init
611  logical, intent(out) :: active
612  real(kind=kreal), intent(in) :: mu
613  real(kind=kreal), optional, target :: b(:)
614 
615  real(kind=kreal) :: distclr
616  integer(kind=kint) :: slave, id, etype
617  integer(kind=kint) :: nn, i, j, iSS, nactive
618  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node )
619  real(kind=kreal) :: nlforce, slforce(3)
620  logical :: isin
621  integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
622  !
623  integer, pointer :: indexMaster(:),indexCand(:)
624  integer :: nMaster,idm,nMasterMax,bktID,nCand
625  logical :: is_cand, is_present_B
626  real(kind=kreal), pointer :: bp(:)
627 
628  if( is_init ) then
629  distclr = embed%cparam%DISTCLR_INIT
630  else
631  distclr = embed%cparam%DISTCLR_FREE
632  active = .false.
633  do i= 1, size(embed%slave)
634  if( embed%states(i)%state==contactstick ) then
635  active = .true.
636  exit
637  endif
638  enddo
639  return
640  endif
641 
642  allocate(contact_surf(size(nodeid)))
643  allocate(states_prev(size(embed%slave)))
644  contact_surf(:) = huge(0)
645  do i = 1, size(embed%slave)
646  states_prev(i) = embed%states(i)%state
647  enddo
648 
649  call update_surface_box_info( embed%master, currpos )
650  call update_surface_bucket_info( embed%master, embed%master_bktDB )
651 
652  ! for gfortran-10: optional parameter seems not allowed within omp parallel
653  is_present_b = present(b)
654  if( is_present_b ) bp => b
655 
656  !$omp parallel do &
657  !$omp& default(none) &
658  !$omp& private(i,slave,slforce,id,nlforce,coord,indexMaster,nMaster,nn,j,iSS,elem,is_cand,idm,etype,isin, &
659  !$omp& bktID,nCand,indexCand) &
660  !$omp& firstprivate(nMasterMax,is_present_B) &
661  !$omp& shared(embed,ndforce,flag_ctAlgo,infoCTChange,currpos,currdisp,mu,nodeID,elemID,Bp,distclr,contact_surf) &
662  !$omp& reduction(.or.:active) &
663  !$omp& schedule(dynamic,1)
664  do i= 1, size(embed%slave)
665  slave = embed%slave(i)
666  if( embed%states(i)%state==contactfree ) then
667  coord(:) = currpos(3*slave-2:3*slave)
668 
669  ! get master candidates from bucketDB
670  bktid = bucketdb_getbucketid(embed%master_bktDB, coord)
671  ncand = bucketdb_getnumcand(embed%master_bktDB, bktid)
672  if (ncand == 0) cycle
673  allocate(indexcand(ncand))
674  call bucketdb_getcand(embed%master_bktDB, bktid, ncand, indexcand)
675 
676  nmastermax = ncand
677  allocate(indexmaster(nmastermax))
678  nmaster = 0
679 
680  ! narrow down candidates
681  do idm= 1, ncand
682  id = indexcand(idm)
683  if (.not. is_in_surface_box( embed%master(id), coord(1:3), embed%cparam%BOX_EXP_RATE )) cycle
684  nmaster = nmaster + 1
685  indexmaster(nmaster) = id
686  enddo
687  deallocate(indexcand)
688 
689  if(nmaster == 0) then
690  deallocate(indexmaster)
691  cycle
692  endif
693 
694  do idm = 1,nmaster
695  id = indexmaster(idm)
696  etype = embed%master(id)%etype
697  if( mod(etype,10) == 2 ) etype = etype - 1 !search by 1st-order shape function
698  nn = getnumberofnodes(etype)
699  do j=1,nn
700  iss = embed%master(id)%nodes(j)
701  elem(1:3,j)=currpos(3*iss-2:3*iss)
702  enddo
703  call project_point2solidelement( coord,etype,nn,elem,embed%master(id)%reflen,embed%states(i), &
704  isin,distclr,localclr=embed%cparam%CLEARANCE )
705  if( .not. isin ) cycle
706  embed%states(i)%surface = id
707  embed%states(i)%multiplier(:) = 0.d0
708  contact_surf(embed%slave(i)) = elemid(embed%master(id)%eid)
709  write(*,'(A,i10,A,i10,A,3f7.3,A,i6)') "Node",nodeid(slave)," embeded to element", &
710  elemid(embed%master(id)%eid), " at ",embed%states(i)%lpos(:)," rank=",hecmw_comm_get_rank()
711  exit
712  enddo
713  deallocate(indexmaster)
714  endif
715  enddo
716  !$omp end parallel do
717 
718  call hecmw_contact_comm_allreduce_i(embed%comm, contact_surf, hecmw_min)
719  nactive = 0
720  do i = 1, size(embed%slave)
721  if (embed%states(i)%state /= contactfree) then ! any slave in contact
722  id = embed%states(i)%surface
723  if (abs(contact_surf(embed%slave(i))) /= elemid(embed%master(id)%eid)) then ! that is in contact with other surface
724  embed%states(i)%state = contactfree ! should be freed
725  write(*,'(A,i10,A,i10,A,i6,A,i6,A)') "Node",nodeid(embed%slave(i))," contact with element", &
726  & elemid(embed%master(id)%eid), " in rank",hecmw_comm_get_rank()," freed due to duplication"
727  else
728  nactive = nactive + 1
729  endif
730  endif
731  if (states_prev(i) == contactfree .and. embed%states(i)%state /= contactfree) then
732  infoctchange%free2contact = infoctchange%free2contact + 1
733  elseif (states_prev(i) /= contactfree .and. embed%states(i)%state == contactfree) then
734  infoctchange%contact2free = infoctchange%contact2free + 1
735  endif
736  enddo
737  active = (nactive > 0)
738  deallocate(contact_surf)
739  deallocate(states_prev)
740 end subroutine scan_embed_state
741 
742 
744  subroutine cal_node_normal( csurf, isin, surf, currpos, lpos, normal )
745  use elementinfo, only:getvertexcoord, surfacenormal
746  integer, intent(in) :: csurf
747  integer, intent(in) :: isin
748  type(tsurfelement), intent(in) :: surf(:)
749  real(kind=kreal), intent(in) :: currpos(:)
750  real(kind=kreal), intent(in) :: lpos(:)
751  real(kind=kreal), intent(out) :: normal(3)
752  integer(kind=kint) :: cnode, i, j, cnt, nd1, gn, etype, iSS, nn,cgn
753  real(kind=kreal) :: cnpos(2), elem(3, l_max_elem_node )
754  integer(kind=kint) :: cnode1, cnode2, gn1, gn2, nsurf, cgn1, cgn2, isin_n
755  real(kind=kreal) :: x=0, normal_n(3), lpos_n(2)
756 
757  if( 1 <= isin .and. isin <= 4 ) then ! corner
758  cnode = isin
759  gn = surf(csurf)%nodes(cnode)
760  etype = surf(csurf)%etype
761  call getvertexcoord( etype, cnode, cnpos )
762  nn = size( surf(csurf)%nodes )
763  do j=1,nn
764  iss = surf(csurf)%nodes(j)
765  elem(1:3,j)=currpos(3*iss-2:3*iss)
766  enddo
767  normal = surfacenormal( etype, nn, cnpos, elem )
768  cnt = 1
769  do i=1,surf(csurf)%n_neighbor
770  nd1 = surf(csurf)%neighbor(i)
771  nn = size( surf(nd1)%nodes )
772  etype = surf(nd1)%etype
773  cgn = 0
774  do j=1,nn
775  iss = surf(nd1)%nodes(j)
776  elem(1:3,j)=currpos(3*iss-2:3*iss)
777  if( iss==gn ) cgn=j
778  enddo
779  if( cgn>0 ) then
780  call getvertexcoord( etype, cgn, cnpos )
781  !normal = normal+SurfaceNormal( etype, nn, cnpos, elem )
782  normal_n = surfacenormal( etype, nn, cnpos, elem )
783  normal = normal+normal_n
784  cnt = cnt+1
785  endif
786  enddo
787  !normal = normal/cnt !!-???
788  elseif( 12 <= isin .and. isin <= 41 ) then ! edge
789  cnode1 = isin / 10
790  cnode2 = mod(isin, 10)
791  gn1 = surf(csurf)%nodes(cnode1)
792  gn2 = surf(csurf)%nodes(cnode2)
793  etype = surf(csurf)%etype
794  nn = size( surf(csurf)%nodes )
795  do j=1,nn
796  iss = surf(csurf)%nodes(j)
797  elem(1:3,j)=currpos(3*iss-2:3*iss)
798  enddo
799  normal = surfacenormal( etype, nn, lpos, elem )
800  select case (etype)
801  case (fe_tri3n, fe_tri6n, fe_tri6nc)
802  if ( isin==12 ) then; x=lpos(2)-lpos(1)
803  elseif( isin==23 ) then; x=1.d0-2.d0*lpos(2)
804  elseif( isin==31 ) then; x=2.d0*lpos(1)-1.d0
805  else; stop "Error: cal_node_normal: invalid isin"
806  endif
807  case (fe_quad4n, fe_quad8n)
808  if ( isin==12 ) then; x=lpos(1)
809  elseif( isin==23 ) then; x=lpos(2)
810  elseif( isin==34 ) then; x=-lpos(1)
811  elseif( isin==41 ) then; x=-lpos(2)
812  else; stop "Error: cal_node_normal: invalid isin"
813  endif
814  end select
815  ! find neighbor surf that includes cnode1 and cnode2
816  nsurf = 0
817  neib_loop: do i=1, surf(csurf)%n_neighbor
818  nd1 = surf(csurf)%neighbor(i)
819  nn = size( surf(nd1)%nodes )
820  etype = surf(nd1)%etype
821  cgn1 = 0
822  cgn2 = 0
823  do j=1,nn
824  iss = surf(nd1)%nodes(j)
825  elem(1:3,j)=currpos(3*iss-2:3*iss)
826  if( iss==gn1 ) cgn1=j
827  if( iss==gn2 ) cgn2=j
828  enddo
829  if( cgn1>0 .and. cgn2>0 ) then
830  nsurf = nd1
831  isin_n = 10*cgn2 + cgn1
832  x = -x
833  select case (etype)
834  case (fe_tri3n, fe_tri6n, fe_tri6nc)
835  if ( isin_n==12 ) then; lpos_n(1)=0.5d0*(1.d0-x); lpos_n(2)=0.5d0*(1.d0+x)
836  elseif( isin_n==23 ) then; lpos_n(1)=0.d0; lpos_n(2)=0.5d0*(1.d0-x)
837  elseif( isin_n==31 ) then; lpos_n(1)=0.5d0*(1.d0+x); lpos_n(2)=0.d0
838  else; stop "Error: cal_node_normal: invalid isin_n"
839  endif
840  case (fe_quad4n, fe_quad8n)
841  if ( isin_n==12 ) then; lpos_n(1)= x; lpos_n(2)=-1.d0
842  elseif( isin_n==23 ) then; lpos_n(1)= 1.d0; lpos_n(2)= x
843  elseif( isin_n==34 ) then; lpos_n(1)=-x; lpos_n(2)= 1.d0
844  elseif( isin_n==41 ) then; lpos_n(1)=-1.d0; lpos_n(2)=-x
845  else; stop "Error: cal_node_normal: invalid isin_n"
846  endif
847  end select
848  !normal = normal + SurfaceNormal( etype, nn, lpos_n, elem )
849  normal_n = surfacenormal( etype, nn, lpos_n, elem )
850  normal = normal+normal_n
851  exit neib_loop
852  endif
853  enddo neib_loop
854  !if( nsurf==0 ) write(0,*) "Warning: cal_node_normal: neighbor surf not found"
855  !normal = normal/2
856  endif
857  normal = normal/ dsqrt( dot_product( normal, normal ) )
858  end subroutine cal_node_normal
859 
861  subroutine track_contact_position( flag_ctAlgo, nslave, contact, currpos, currdisp, mu, infoCTChange, nodeID, elemID, B )
862  character(len=9), intent(in) :: flag_ctAlgo
863  integer, intent(in) :: nslave
864  type( tcontact ), intent(inout) :: contact
865  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
866  real(kind=kreal), intent(in) :: currpos(:)
867  real(kind=kreal), intent(in) :: currdisp(:)
868  real(kind=kreal), intent(in) :: mu
869  integer(kind=kint), intent(in) :: nodeID(:)
870  integer(kind=kint), intent(in) :: elemID(:)
871  real(kind=kreal), intent(inout) :: b(:)
872 
873  integer(kind=kint) :: slave, sid0, sid, etype
874  integer(kind=kint) :: nn, i, j, iSS
875  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node ), elem0(3, l_max_elem_node )
876  logical :: isin
877  real(kind=kreal) :: opos(2), odirec(3)
878  integer(kind=kint) :: bktID, nCand, idm
879  integer(kind=kint), allocatable :: indexCand(:)
880 
881  sid = 0
882 
883  slave = contact%slave(nslave)
884  coord(:) = currpos(3*slave-2:3*slave)
886  sid0 = contact%states(nslave)%surface
887  opos = contact%states(nslave)%lpos(1:2)
888  odirec = contact%states(nslave)%direction
889  etype = contact%master(sid0)%etype
890  nn = getnumberofnodes( etype )
891  do j=1,nn
892  iss = contact%master(sid0)%nodes(j)
893  elem(1:3,j)=currpos(3*iss-2:3*iss)
894  elem0(1:3,j)=currpos(3*iss-2:3*iss)-currdisp(3*iss-2:3*iss)
895  enddo
896  call project_point2element( coord,etype,nn,elem,contact%master(sid0)%reflen,contact%states(nslave), &
897  isin,contact%cparam%DISTCLR_NOCHECK,contact%states(nslave)%lpos(1:2),contact%cparam%CLR_SAME_ELEM )
898  if( .not. isin ) then
899  do i=1, contact%master(sid0)%n_neighbor
900  sid = contact%master(sid0)%neighbor(i)
901  etype = contact%master(sid)%etype
902  nn = getnumberofnodes( etype )
903  do j=1,nn
904  iss = contact%master(sid)%nodes(j)
905  elem(1:3,j)=currpos(3*iss-2:3*iss)
906  enddo
907  call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
908  isin,contact%cparam%DISTCLR_NOCHECK,localclr=contact%cparam%CLEARANCE )
909  if( isin ) then
910  contact%states(nslave)%surface = sid
911  exit
912  endif
913  enddo
914  endif
915 
916  if( .not. isin ) then ! such case is considered to rarely or never occur
917  write(*,*) 'Warning: contact moved beyond neighbor elements'
918  ! get master candidates from bucketDB
919  bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
920  ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
921  if (ncand > 0) then
922  allocate(indexcand(ncand))
923  call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
924  do idm= 1, ncand
925  sid = indexcand(idm)
926  if( sid==sid0 ) cycle
927  if( associated(contact%master(sid0)%neighbor) ) then
928  if( any(sid==contact%master(sid0)%neighbor(:)) ) cycle
929  endif
930  if (.not. is_in_surface_box( contact%master(sid), coord(1:3), contact%cparam%BOX_EXP_RATE )) cycle
931  etype = contact%master(sid)%etype
932  nn = size( contact%master(sid)%nodes )
933  do j=1,nn
934  iss = contact%master(sid)%nodes(j)
935  elem(1:3,j)=currpos(3*iss-2:3*iss)
936  enddo
937  call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
938  isin,contact%cparam%DISTCLR_NOCHECK,localclr=contact%cparam%CLEARANCE )
939  if( isin ) then
940  contact%states(nslave)%surface = sid
941  exit
942  endif
943  enddo
944  deallocate(indexcand)
945  endif
946  endif
947 
948  if( isin ) then
949  if( contact%states(nslave)%surface==sid0 ) then
950  if(any(dabs(contact%states(nslave)%lpos(1:2)-opos(:)) >= contact%cparam%CLR_DIFFLPOS)) then
951  !$omp atomic
952  infoctchange%contact2difflpos = infoctchange%contact2difflpos + 1
953  endif
954  else
955  write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3)') "Node",nodeid(slave)," move to contact with", &
956  elemid(contact%master(sid)%eid), " with distance ", &
957  contact%states(nslave)%distance," at ",contact%states(nslave)%lpos(1:2)
958  !$omp atomic
959  infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
960  if( flag_ctalgo=='ALagrange' ) &
961  call reset_contact_force( contact, currpos, nslave, sid0, opos, odirec, b )
962  endif
963  if( flag_ctalgo=='SLagrange' ) call update_tangentforce(etype,nn,elem0,elem,contact%states(nslave))
964  iss = isinsideelement( etype, contact%states(nslave)%lpos(1:2), contact%cparam%CLR_CAL_NORM )
965  if( iss>0 ) &
966  call cal_node_normal( contact%states(nslave)%surface, iss, contact%master, currpos, &
967  contact%states(nslave)%lpos(1:2), contact%states(nslave)%direction(:) )
968  else if( .not. isin ) then
969  write(*,'(A,i10,A)') "Node",nodeid(slave)," move out of contact"
970  contact%states(nslave)%state = contactfree
971  contact%states(nslave)%multiplier(:) = 0.d0
972  endif
973 
974  end subroutine track_contact_position
975 
977  subroutine update_direction( nslave, contact, currpos )
978  integer, intent(in) :: nslave
979  type( tcontact ), intent(inout) :: contact
980  real(kind=kreal), intent(in) :: currpos(:)
981 
982  integer(kind=kint) :: slave, sid0, sid, etype
983  integer(kind=kint) :: nn, i, j, iSS
984  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node ), elem0(3, l_max_elem_node )
985  logical :: isin
986  real(kind=kreal) :: opos(2), odirec(3)
987  integer(kind=kint) :: bktID, nCand, idm
988  integer(kind=kint), allocatable :: indexCand(:)
989  type(tcontactstate) :: cstate_tmp
990 
991  sid = 0
992 
993  slave = contact%slave(nslave)
994  coord(:) = currpos(3*slave-2:3*slave)
996  sid0 = contact%states(nslave)%surface
997  opos = contact%states(nslave)%lpos(1:2)
998  odirec = contact%states(nslave)%direction
999  etype = contact%master(sid0)%etype
1000  nn = getnumberofnodes( etype )
1001  do j=1,nn
1002  iss = contact%master(sid0)%nodes(j)
1003  elem(1:3,j)=currpos(3*iss-2:3*iss)
1004  enddo
1005 
1006  cstate_tmp%state = contact%states(nslave)%state
1007  cstate_tmp%surface = sid0
1008  call project_point2element( coord,etype,nn,elem,contact%master(sid0)%reflen,cstate_tmp, &
1009  isin,contact%cparam%DISTCLR_NOCHECK,contact%states(nslave)%lpos,contact%cparam%CLR_SAME_ELEM )
1010 
1011  if( isin ) then
1012  iss = isinsideelement( etype, cstate_tmp%lpos, contact%cparam%CLR_CAL_NORM )
1013  if( iss>0 ) &
1014  call cal_node_normal( cstate_tmp%surface, iss, contact%master, currpos, &
1015  cstate_tmp%lpos, cstate_tmp%direction(:) )
1016  endif
1017 
1018  contact%states(nslave)%direction = cstate_tmp%direction
1019 
1020  end subroutine
1021 
1022 
1025  subroutine reset_contact_force( contact, currpos, lslave, omaster, opos, odirec, B )
1026  type( tcontact ), intent(inout) :: contact
1027  real(kind=kreal), intent(in) :: currpos(:)
1028  integer, intent(in) :: lslave
1029  integer, intent(in) :: omaster
1030  real(kind=kreal), intent(in) :: opos(2)
1031  real(kind=kreal), intent(in) :: odirec(3)
1032  real(kind=kreal), intent(inout) :: b(:)
1033 
1034  integer(kind=kint) :: slave, etype, master
1035  integer(kind=kint) :: nn, j, iSS
1036  real(kind=kreal) :: nrlforce, fcoeff, tangent(3,2)
1037  real(kind=kreal) :: elemcrd(3, l_max_elem_node )
1038  real(kind=kreal) :: shapefunc(l_max_surface_node)
1039  real(kind=kreal) :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
1040  real(kind=kreal) :: fric(2), f3(l_max_elem_node*3+3)
1041  integer(kind=kint) :: i, idx0
1042 
1043  slave = contact%slave(lslave)
1044  fcoeff = contact%fcoeff
1045  ! clear contact force in former contacted element
1046  nrlforce = contact%states(lslave)%multiplier(1)
1047  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)+nrlforce*odirec
1048  nn = size( contact%master(omaster)%nodes )
1049  etype = contact%master(omaster)%etype
1050  call getshapefunc( etype, opos(:), shapefunc )
1051  do j=1,nn
1052  iss = contact%master(omaster)%nodes(j)
1053  ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)-nrlforce*shapefunc(j)*odirec
1054  idx0 = 3*(iss-1)
1055  do i=1,3
1056  !$omp atomic
1057  b(idx0+i) = b(idx0+i)-nrlforce*shapefunc(j)*odirec(i)
1058  enddo
1059  enddo
1060  if( fcoeff/=0.d0 ) then
1061  do j=1,nn
1062  iss = contact%master(omaster)%nodes(j)
1063  elemcrd(:,j) = currpos(3*iss-2:3*iss)
1064  enddo
1065  call dispincrematrix( opos(:), etype, nn, elemcrd, tangent, &
1066  metric, dispmat )
1067  fric(1:2) = contact%states(lslave)%multiplier(2:3)
1068  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
1069  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)+f3(1:3)
1070  do j=1,nn
1071  iss = contact%master(omaster)%nodes(j)
1072  ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)+f3(3*j+1:3*j+3)
1073  idx0 = 3*(iss-1)
1074  do i=1,3
1075  !$omp atomic
1076  b(idx0+i) = b(idx0+i)+f3(3*j+i)
1077  enddo
1078  enddo
1079  endif
1080 
1081  ! reset contact force in new contacting element
1082  master = contact%states(lslave)%surface
1083  nn = size( contact%master(master)%nodes )
1084  etype = contact%master(master)%etype
1085  call getshapefunc( etype, contact%states(lslave)%lpos(1:2), shapefunc )
1086  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-nrlforce*contact%states(lslave)%direction
1087  do j=1,nn
1088  iss = contact%master(master)%nodes(j)
1089  ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)+nrlforce* &
1090  ! shapefunc(j)*contact%states(lslave)%direction
1091  idx0 = 3*(iss-1)
1092  do i=1,3
1093  !$omp atomic
1094  b(idx0+i) = b(idx0+i)+nrlforce* &
1095  shapefunc(j)*contact%states(lslave)%direction(i)
1096  enddo
1097  enddo
1098  if( fcoeff/=0.d0 ) then
1099  do j=1,nn
1100  iss = contact%master(master)%nodes(j)
1101  elemcrd(:,j) = currpos(3*iss-2:3*iss)
1102  enddo
1103  call dispincrematrix( contact%states(lslave)%lpos(1:2), etype, nn, elemcrd, tangent, &
1104  metric, dispmat )
1105  fric(1:2) = contact%states(lslave)%multiplier(2:3)
1106  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
1107  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-f3(1:3)
1108  do j=1,nn
1109  iss = contact%master(master)%nodes(j)
1110  ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)-f3(3*j+1:3*j+3)
1111  idx0 = 3*(iss-1)
1112  do i=1,3
1113  !$omp atomic
1114  b(idx0+i) = b(idx0+i)-f3(3*j+i)
1115  enddo
1116  enddo
1117  endif
1118 
1119  end subroutine reset_contact_force
1120 
1124  subroutine calcu_contact_force0( contact, coord, disp, ddisp, fcoeff, mu, &
1125  mut, B )
1126  type( tcontact ), intent(inout) :: contact
1127  real(kind=kreal), intent(in) :: coord(:)
1128  real(kind=kreal), intent(in) :: disp(:)
1129  real(kind=kreal), intent(in) :: ddisp(:)
1130  real(kind=kreal), intent(in) :: fcoeff
1131  real(kind=kreal), intent(in) :: mu, mut
1132  real(kind=kreal), intent(inout) :: b(:)
1133 
1134  integer(kind=kint) :: slave, etype, master
1135  integer(kind=kint) :: nn, i, j, iSS
1136  real(kind=kreal) :: nrlforce, elemdisp(3,l_max_elem_node), tangent(3,2)
1137  real(kind=kreal) :: dg(3), elemg(3), elemcrd(3, l_max_elem_node )
1138  real(kind=kreal) :: dgn, dxi(2), dxy(2), shapefunc(l_max_surface_node)
1139  real(kind=kreal) :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
1140  real(kind=kreal) :: fric(2), f3(3*l_max_elem_node+3), edisp(3*l_max_elem_node+3)
1141 
1142  do i= 1, size(contact%slave)
1143  if( contact%states(i)%state==contactfree ) cycle ! not in contact
1144  slave = contact%slave(i)
1145  edisp(1:3) = ddisp(3*slave-2:3*slave)
1146  master = contact%states(i)%surface
1147 
1148  nn = size( contact%master(master)%nodes )
1149  etype = contact%master(master)%etype
1150  do j=1,nn
1151  iss = contact%master(master)%nodes(j)
1152  elemdisp(:,j) = ddisp(3*iss-2:3*iss)
1153  edisp(3*j+1:3*j+3) = ddisp(3*iss-2:3*iss)
1154  elemcrd(:,j) = coord(3*iss-2:3*iss)+disp(3*iss-2:3*iss)
1155  enddo
1156  call getshapefunc( etype, contact%states(i)%lpos(1:2), shapefunc )
1157 
1158  ! normal component
1159  elemg = 0.d0
1160  do j=1,nn
1161  elemg(:) = elemg(:)+shapefunc(j)*elemdisp(:,j)
1162  enddo
1163  dg(:) = ddisp(3*slave-2:3*slave) - elemg(:)
1164  dgn = dot_product( contact%states(i)%direction, dg )
1165  nrlforce = contact%states(i)%multiplier(1)- mu*(contact%states(i)%wkdist-dgn)
1166  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-nrlforce*contact%states(i)%direction
1167 
1168  do j=1,nn
1169  iss = contact%master(master)%nodes(j)
1170  b(3*iss-2:3*iss) = b(3*iss-2:3*iss)+nrlforce* &
1171  shapefunc(j)*contact%states(i)%direction
1172  enddo
1173 
1174  if( fcoeff==0.d0 ) cycle
1175  ! tangent component
1176  call dispincrematrix( contact%states(i)%lpos(1:2), etype, nn, elemcrd, tangent, &
1177  metric, dispmat )
1178  dxi(1) = dot_product( dispmat(1,1:nn*3+3), edisp(1:nn*3+3) )
1179  dxi(2) = dot_product( dispmat(2,1:nn*3+3), edisp(1:nn*3+3) )
1180  dxy(:) = matmul( metric, dxi )
1181  fric(1:2) = contact%states(i)%multiplier(2:3) + mut*dxy(1:2)
1182  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
1183 
1184  if( contact%states(i)%state==contactslip ) then
1185  dgn = dsqrt( f3(1)*f3(1)+f3(2)*f3(2)+f3(3)*f3(3) )
1186  f3(:) = f3(:)*fcoeff*contact%states(i)%multiplier(1)/dgn
1187  endif
1188  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-f3(1:3)
1189  do j=1,nn
1190  iss = contact%master(master)%nodes(j)
1191  b(3*iss-2:3*iss) = b(3*iss-2:3*iss)-f3(3*j+1:3*j+3)
1192  enddo
1193  enddo
1194  end subroutine calcu_contact_force0
1195 
1199  subroutine calcu_tied_force0( contact, disp, ddisp, mu, B )
1200  type( tcontact ), intent(inout) :: contact
1201  real(kind=kreal), intent(in) :: disp(:)
1202  real(kind=kreal), intent(in) :: ddisp(:)
1203  real(kind=kreal), intent(in) :: mu
1204  real(kind=kreal), intent(inout) :: b(:)
1205 
1206  integer(kind=kint) :: slave, etype, master
1207  integer(kind=kint) :: nn, i, j, iss
1208  real(kind=kreal) :: nrlforce(3)
1209  real(kind=kreal) :: dg(3)
1210  real(kind=kreal) :: shapefunc(l_max_surface_node)
1211  real(kind=kreal) :: edisp(3*l_max_elem_node+3)
1212 
1213  do i= 1, size(contact%slave)
1214  if( contact%states(i)%state==contactfree ) cycle ! not in contact
1215  slave = contact%slave(i)
1216  edisp(1:3) = disp(3*slave-2:3*slave)+ddisp(3*slave-2:3*slave)
1217  master = contact%states(i)%surface
1218 
1219  nn = size( contact%master(master)%nodes )
1220  etype = contact%master(master)%etype
1221  do j=1,nn
1222  iss = contact%master(master)%nodes(j)
1223  edisp(3*j+1:3*j+3) = disp(3*iss-2:3*iss)+ddisp(3*iss-2:3*iss)
1224  enddo
1225  call getshapefunc( etype, contact%states(i)%lpos(1:2), shapefunc )
1226 
1227  ! normal component
1228  dg(1:3) = edisp(1:3)
1229  do j=1,nn
1230  dg(1:3) = dg(1:3)-shapefunc(j)*edisp(3*j+1:3*j+3)
1231  enddo
1232 
1233  nrlforce(1:3) = contact%states(i)%multiplier(1:3)+mu*dg(1:3)
1234 
1235  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-nrlforce(1:3)
1236  do j=1,nn
1237  iss = contact%master(master)%nodes(j)
1238  b(3*iss-2:3*iss) = b(3*iss-2:3*iss) + shapefunc(j)*nrlforce(1:3)
1239  enddo
1240  enddo
1241 
1242  end subroutine
1243 
1246  subroutine update_contact_multiplier( contact, coord, disp, ddisp, fcoeff, mu, mut, &
1247  gnt, ctchanged )
1248  type( tcontact ), intent(inout) :: contact
1249  real(kind=kreal), intent(in) :: coord(:)
1250  real(kind=kreal), intent(in) :: disp(:)
1251  real(kind=kreal), intent(in) :: ddisp(:)
1252  real(kind=kreal), intent(in) :: fcoeff
1253  real(kind=kreal), intent(in) :: mu, mut
1254  real(kind=kreal), intent(out) :: gnt(2)
1255  logical, intent(inout) :: ctchanged
1256 
1257  integer(kind=kint) :: slave, etype, master
1258  integer(kind=kint) :: nn, i, j, iss, cnt
1259  real(kind=kreal) :: elemdisp(3,l_max_elem_node), tangent(3,2)
1260  real(kind=kreal) :: dg(3), elemg(3), elemcrd(3, l_max_elem_node )
1261  real(kind=kreal) :: dgn, dxi(2), dxy(2), shapefunc(l_max_surface_node)
1262  real(kind=kreal) :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
1263  real(kind=kreal) :: lgnt(2), fric(2), f3(3*l_max_elem_node+3), edisp(3*l_max_elem_node+3)
1264 
1265  cnt =0; lgnt(:)=0.d0
1266  do i= 1, size(contact%slave)
1267  if( contact%states(i)%state==contactfree ) cycle ! not in contact
1268  slave = contact%slave(i)
1269  edisp(1:3) = ddisp(3*slave-2:3*slave)
1270  master = contact%states(i)%surface
1271 
1272  nn = size( contact%master(master)%nodes )
1273  etype = contact%master(master)%etype
1274  do j=1,nn
1275  iss = contact%master(master)%nodes(j)
1276  elemdisp(:,j) = ddisp(3*iss-2:3*iss)
1277  edisp(3*j+1:3*j+3) = ddisp(3*iss-2:3*iss)
1278  elemcrd(:,j) = coord(3*iss-2:3*iss)+disp(3*iss-2:3*iss)
1279  enddo
1280  call getshapefunc( etype, contact%states(i)%lpos(1:2), shapefunc )
1281 
1282  ! normal component
1283  elemg = 0.d0
1284  do j=1,nn
1285  elemg(:) = elemg(:)+shapefunc(j)*elemdisp(:,j)
1286  enddo
1287  dg(:) = ddisp(3*slave-2:3*slave) - elemg(:)
1288  dgn = dot_product( contact%states(i)%direction, dg )
1289  contact%states(i)%wkdist = contact%states(i)%wkdist-dgn
1290  contact%states(i)%multiplier(1) = contact%states(i)%multiplier(1)- mu*contact%states(i)%wkdist
1291  contact%states(i)%distance = contact%states(i)%distance - dgn
1292  cnt = cnt+1
1293  lgnt(1) = lgnt(1)- contact%states(i)%wkdist
1294 
1295  if( fcoeff==0.d0 ) cycle
1296  ! tangent component
1297  call dispincrematrix( contact%states(i)%lpos(1:2), etype, nn, elemcrd, tangent, &
1298  metric, dispmat )
1299  dxi(1) = dot_product( dispmat(1,1:nn*3+3), edisp(1:nn*3+3) )
1300  dxi(2) = dot_product( dispmat(2,1:nn*3+3), edisp(1:nn*3+3) )
1301  dxy(:) = matmul( metric, dxi )
1302  fric(1:2) = contact%states(i)%multiplier(2:3) + mut*dxy(1:2)
1303  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
1304  dgn = dsqrt( f3(1)*f3(1)+f3(2)*f3(2)+f3(3)*f3(3) )
1305  if( contact%states(i)%multiplier(1)>0.d0 ) then
1306  if( dgn > fcoeff*contact%states(i)%multiplier(1) ) then
1307  if( contact%states(i)%state==contactstick ) then
1308  ctchanged= .true.
1309  print *, "Node", slave, "to slip state",dgn, fcoeff*contact%states(i)%multiplier(1)
1310  endif
1311  contact%states(i)%state = contactslip
1312  fric(:) = fric(:)*fcoeff*contact%states(i)%multiplier(1)/dgn
1313  else
1314  if( contact%states(i)%state==contactslip ) then
1315  ctchanged= .true.
1316  print *, "Node", slave, "to stick state",dgn, fcoeff*contact%states(i)%multiplier(1)
1317  endif
1318  contact%states(i)%state = contactstick
1319  endif
1320  endif
1321  contact%states(i)%multiplier(2:3) = fric(:)
1322  dxy(:) = matmul( dg, tangent )
1323  lgnt(2) = lgnt(2)+dsqrt( dxy(1)*dxy(1)+dxy(2)*dxy(2) )
1324  enddo
1325  if(cnt>0) lgnt(:) = lgnt(:)/cnt
1326  gnt = gnt + lgnt
1327  end subroutine update_contact_multiplier
1328 
1331  subroutine update_tied_multiplier( contact, disp, ddisp, mu, ctchanged )
1332  type( tcontact ), intent(inout) :: contact
1333  real(kind=kreal), intent(in) :: disp(:)
1334  real(kind=kreal), intent(in) :: ddisp(:)
1335  real(kind=kreal), intent(in) :: mu
1336  logical, intent(inout) :: ctchanged
1337 
1338  integer(kind=kint) :: slave, etype, master
1339  integer(kind=kint) :: nn, i, j, iss, cnt
1340  real(kind=kreal) :: dg(3), dgmax
1341  real(kind=kreal) :: shapefunc(l_max_surface_node)
1342  real(kind=kreal) :: edisp(3*l_max_elem_node+3)
1343 
1344  do i= 1, size(contact%slave)
1345  if( contact%states(i)%state==contactfree ) cycle ! not in contact
1346  slave = contact%slave(i)
1347  edisp(1:3) = disp(3*slave-2:3*slave)+ddisp(3*slave-2:3*slave)
1348  master = contact%states(i)%surface
1349 
1350  nn = size( contact%master(master)%nodes )
1351  etype = contact%master(master)%etype
1352  do j=1,nn
1353  iss = contact%master(master)%nodes(j)
1354  edisp(3*j+1:3*j+3) = disp(3*iss-2:3*iss)+ddisp(3*iss-2:3*iss)
1355  enddo
1356  call getshapefunc( etype, contact%states(i)%lpos(1:2), shapefunc )
1357 
1358  ! normal component
1359  dg(1:3) = edisp(1:3)
1360  do j=1,nn
1361  dg(1:3) = dg(1:3)-shapefunc(j)*edisp(3*j+1:3*j+3)
1362  enddo
1363 
1364  contact%states(i)%multiplier(1:3) = contact%states(i)%multiplier(1:3) + mu*dg(1:3)
1365 
1366  ! check if tied constraint converged
1367  dgmax = 0.d0
1368  do j=1,(nn+1)*3
1369  dgmax = dgmax + dabs(edisp(j))
1370  enddo
1371  dgmax = dgmax/dble((nn+1)*3)
1372  do j=1,3
1373  if( dabs(dg(j))/dmax1(1.d0,dgmax) > 1.d-3 ) ctchanged = .true.
1374  enddo
1375 
1376  enddo
1377  end subroutine
1378 
1380  subroutine ass_contact_force( contact, coord, disp, B )
1381  type( tcontact ), intent(in) :: contact
1382  real(kind=kreal), intent(in) :: coord(:)
1383  real(kind=kreal), intent(in) :: disp(:)
1384  real(kind=kreal), intent(inout) :: b(:)
1385 
1386  integer(kind=kint) :: slave, etype, master
1387  integer(kind=kint) :: nn, i, j, iss
1388  real(kind=kreal) :: fcoeff, nrlforce, tangent(3,2)
1389  real(kind=kreal) :: elemcrd(3, l_max_elem_node )
1390  real(kind=kreal) :: shapefunc(l_max_surface_node)
1391  real(kind=kreal) :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
1392  real(kind=kreal) :: fric(2), f3(l_max_elem_node*3+3)
1393  fcoeff = contact%fcoeff
1394  do i= 1, size(contact%slave)
1395  if( contact%states(i)%state==contactfree ) cycle ! not in contact
1396  slave = contact%slave(i)
1397  master = contact%states(i)%surface
1398 
1399  nn = size( contact%master(master)%nodes )
1400  etype = contact%master(master)%etype
1401  call getshapefunc( etype, contact%states(i)%lpos(1:2), shapefunc )
1402 
1403  ! normal component
1404  nrlforce = contact%states(i)%multiplier(1)
1405  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-nrlforce*contact%states(i)%direction
1406 
1407  do j=1,nn
1408  iss = contact%master(master)%nodes(j)
1409  b(3*iss-2:3*iss) = b(3*iss-2:3*iss)+nrlforce* &
1410  shapefunc(j)*contact%states(i)%direction
1411  enddo
1412 
1413  if( fcoeff==0.d0 ) cycle
1414  ! tangent component
1415  do j=1,nn
1416  iss = contact%master(master)%nodes(j)
1417  elemcrd(:,j) = coord(3*iss-2:3*iss)+disp(3*iss-2:3*iss)
1418  enddo
1419  call dispincrematrix( contact%states(i)%lpos(1:2), etype, nn, elemcrd, tangent, &
1420  metric, dispmat )
1421 
1422  fric(1:2) = contact%states(i)%multiplier(2:3)
1423  f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
1424  b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-f3(1:3)
1425  do j=1,nn
1426  iss = contact%master(master)%nodes(j)
1427  b(3*iss-2:3*iss) = b(3*iss-2:3*iss)-f3(3*j+1:3*j+3)
1428  enddo
1429  enddo
1430 
1431  end subroutine ass_contact_force
1432 
1434  subroutine set_contact_state_vector( contact, dt, relvel_vec, state_vec )
1435  type( tcontact ), intent(in) :: contact
1436  real(kind=kreal), intent(in) :: dt
1437  real(kind=kreal), intent(inout) :: relvel_vec(:)
1438  real(kind=kreal), intent(inout) :: state_vec(:)
1439 
1440  integer(kind=kint) :: i, slave
1441 
1442  do i= 1, size(contact%slave)
1443  slave = contact%slave(i)
1444  if( state_vec(slave) < 0.1d0 .or. contact%states(i)%state > 0 ) &
1445  & state_vec(slave) = dble(contact%states(i)%state)
1446 
1447  if( contact%states(i)%state==contactfree ) cycle ! not in contact
1448  if( dt < 1.d-16 ) cycle ! too small delta t
1449  relvel_vec(3*slave-2:3*slave) = contact%states(i)%reldisp(1:3)/dt
1450  enddo
1451 
1452  end subroutine set_contact_state_vector
1453 
1454  subroutine update_contact_tangentforce( contact )
1455  type( tcontact ), intent(inout) :: contact
1456 
1457  integer(kind=kint) :: i
1458 
1459  do i= 1, size(contact%slave)
1460  if( contact%states(i)%state==contactfree ) then
1461  contact%states(i)%tangentForce(1:3) = 0.d0
1462  contact%states(i)%tangentForce_trial(1:3) = 0.d0
1463  contact%states(i)%tangentForce_final(1:3) = 0.d0
1464  else
1465  contact%states(i)%tangentForce(1:3) = contact%states(i)%tangentForce_final(1:3)
1466  end if
1467  contact%states(i)%tangentForce1(1:3) = contact%states(i)%tangentForce(1:3)
1468  enddo
1469  end subroutine update_contact_tangentforce
1470 
1472  subroutine track_contact_position_exp( nslave, contact, currpos, currdisp, infoCTChange, nodeID, elemID )
1473  integer, intent(in) :: nslave
1474  type( tcontact ), intent(inout) :: contact
1475  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
1476  real(kind=kreal), intent(in) :: currpos(:)
1477  real(kind=kreal), intent(in) :: currdisp(:)
1478  integer(kind=kint), intent(in) :: nodeID(:)
1479  integer(kind=kint), intent(in) :: elemID(:)
1480 
1481  integer(kind=kint) :: slave, sid0, sid, etype
1482  integer(kind=kint) :: nn, i, j, iSS
1483  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node ), elem0(3, l_max_elem_node )
1484  logical :: isin
1485  real(kind=kreal) :: opos(2), odirec(3)
1486  integer(kind=kint) :: bktID, nCand, idm
1487  integer(kind=kint), allocatable :: indexCand(:)
1488 
1489  sid = 0
1490 
1491  slave = contact%slave(nslave)
1492  coord(:) = currpos(3*slave-2:3*slave)
1494  sid0 = contact%states(nslave)%surface
1495  opos = contact%states(nslave)%lpos(1:2)
1496  odirec = contact%states(nslave)%direction
1497  etype = contact%master(sid0)%etype
1498  nn = getnumberofnodes( etype )
1499  do j=1,nn
1500  iss = contact%master(sid0)%nodes(j)
1501  elem(1:3,j)=currpos(3*iss-2:3*iss)
1502  elem0(1:3,j)=currpos(3*iss-2:3*iss)-currdisp(3*iss-2:3*iss)
1503  enddo
1504  call project_point2element( coord,etype,nn,elem,contact%master(sid0)%reflen,contact%states(nslave), &
1505  isin,contact%cparam%DISTCLR_NOCHECK,contact%states(nslave)%lpos(1:2),contact%cparam%CLR_SAME_ELEM )
1506  if( .not. isin ) then
1507  do i=1, contact%master(sid0)%n_neighbor
1508  sid = contact%master(sid0)%neighbor(i)
1509  etype = contact%master(sid)%etype
1510  nn = getnumberofnodes( etype )
1511  do j=1,nn
1512  iss = contact%master(sid)%nodes(j)
1513  elem(1:3,j)=currpos(3*iss-2:3*iss)
1514  enddo
1515  call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
1516  isin,contact%cparam%DISTCLR_NOCHECK,localclr=contact%cparam%CLEARANCE )
1517  if( isin ) then
1518  contact%states(nslave)%surface = sid
1519  exit
1520  endif
1521  enddo
1522  endif
1523 
1524  if( .not. isin ) then ! such case is considered to rarely or never occur
1525  write(*,*) 'Warning: contact moved beyond neighbor elements'
1526  ! get master candidates from bucketDB
1527  bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
1528  ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
1529  if (ncand > 0) then
1530  allocate(indexcand(ncand))
1531  call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
1532  do idm= 1, ncand
1533  sid = indexcand(idm)
1534  if( sid==sid0 ) cycle
1535  if( any(sid==contact%master(sid0)%neighbor(:)) ) cycle
1536  if (.not. is_in_surface_box( contact%master(sid), coord(1:3), contact%cparam%BOX_EXP_RATE )) cycle
1537  etype = contact%master(sid)%etype
1538  nn = size( contact%master(sid)%nodes )
1539  do j=1,nn
1540  iss = contact%master(sid)%nodes(j)
1541  elem(1:3,j)=currpos(3*iss-2:3*iss)
1542  enddo
1543  call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
1544  isin,contact%cparam%DISTCLR_NOCHECK,localclr=contact%cparam%CLEARANCE )
1545  if( isin ) then
1546  contact%states(nslave)%surface = sid
1547  exit
1548  endif
1549  enddo
1550  deallocate(indexcand)
1551  endif
1552  endif
1553 
1554  if( isin ) then
1555  if( contact%states(nslave)%surface==sid0 ) then
1556  if(any(dabs(contact%states(nslave)%lpos(1:2)-opos(:)) >= contact%cparam%CLR_DIFFLPOS)) then
1557  !$omp atomic
1558  infoctchange%contact2difflpos = infoctchange%contact2difflpos + 1
1559  endif
1560  else
1561  write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3)') "Node",nodeid(slave)," move to contact with", &
1562  elemid(contact%master(sid)%eid), " with distance ", &
1563  contact%states(nslave)%distance," at ",contact%states(nslave)%lpos(1:2)
1564  !$omp atomic
1565  infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
1566  endif
1567  iss = isinsideelement( etype, contact%states(nslave)%lpos(1:2), contact%cparam%CLR_CAL_NORM )
1568  if( iss>0 ) then
1569  call cal_node_normal( contact%states(nslave)%surface, iss, contact%master, currpos, &
1570  contact%states(nslave)%lpos(1:2), contact%states(nslave)%direction(:) )
1571  endif
1572  else if( .not. isin ) then
1573  write(*,'(A,i10,A)') "Node",nodeid(slave)," move out of contact"
1574  contact%states(nslave)%state = contactfree
1575  contact%states(nslave)%multiplier(:) = 0.d0
1576  endif
1577 
1578  end subroutine track_contact_position_exp
1579 
1583  subroutine scan_contact_state_exp( contact, currpos, currdisp, infoCTChange, &
1584  nodeID, elemID, is_init, active )
1585  type( tcontact ), intent(inout) :: contact
1586  type( fstr_info_contactchange ), intent(inout) :: infoCTChange
1587  real(kind=kreal), intent(in) :: currpos(:)
1588  real(kind=kreal), intent(in) :: currdisp(:)
1589  integer(kind=kint), intent(in) :: nodeID(:)
1590  integer(kind=kint), intent(in) :: elemID(:)
1591  logical, intent(in) :: is_init
1592  logical, intent(out) :: active
1593 
1594  real(kind=kreal) :: distclr
1595  integer(kind=kint) :: slave, id, etype
1596  integer(kind=kint) :: nn, i, j, iSS, nactive
1597  real(kind=kreal) :: coord(3), elem(3, l_max_elem_node )
1598  real(kind=kreal) :: nlforce
1599  logical :: isin
1600  integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
1601  !
1602  integer, pointer :: indexMaster(:),indexCand(:)
1603  integer :: nMaster,idm,nMasterMax,bktID,nCand
1604  logical :: is_cand
1605 
1606  if( is_init ) then
1607  distclr = contact%cparam%DISTCLR_INIT
1608  else
1609  distclr = contact%cparam%DISTCLR_FREE
1610  if( contact%algtype == contacttied ) then
1611  active = .false.
1612  do i= 1, size(contact%slave)
1613  if( contact%states(i)%state==contactstick ) then
1614  active = .true.
1615  exit
1616  endif
1617  enddo
1618  return
1619  endif
1620  endif
1621 
1622  allocate(contact_surf(size(nodeid)))
1623  allocate(states_prev(size(contact%slave)))
1624  contact_surf(:) = size(elemid)+1
1625  do i = 1, size(contact%slave)
1626  states_prev(i) = contact%states(i)%state
1627  enddo
1628 
1629  call update_surface_box_info( contact%master, currpos )
1630  call update_surface_bucket_info( contact%master, contact%master_bktDB )
1631 
1632  !$omp parallel do &
1633  !$omp& default(none) &
1634  !$omp& private(i,slave,id,nlforce,coord,indexMaster,nMaster,nn,j,iSS,elem,is_cand,idm,etype,isin, &
1635  !$omp& bktID,nCand,indexCand) &
1636  !$omp& firstprivate(nMasterMax) &
1637  !$omp& shared(contact,infoCTChange,currpos,currdisp,nodeID,elemID,distclr,contact_surf) &
1638  !$omp& reduction(.or.:active) &
1639  !$omp& schedule(dynamic,1)
1640  do i= 1, size(contact%slave)
1641  slave = contact%slave(i)
1642  if( contact%states(i)%state==contactstick .or. contact%states(i)%state==contactslip ) then
1643  call track_contact_position_exp( i, contact, currpos, currdisp, infoctchange, nodeid, elemid )
1644  if( contact%states(i)%state /= contactfree ) then
1645  contact_surf(contact%slave(i)) = -contact%states(i)%surface
1646  endif
1647  else if( contact%states(i)%state==contactfree ) then
1648  coord(:) = currpos(3*slave-2:3*slave)
1649 
1650  ! get master candidates from bucketDB
1651  bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
1652  ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
1653  if (ncand == 0) cycle
1654  allocate(indexcand(ncand))
1655  call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
1656 
1657  nmastermax = ncand
1658  allocate(indexmaster(nmastermax))
1659  nmaster = 0
1660 
1661  ! narrow down candidates
1662  do idm= 1, ncand
1663  id = indexcand(idm)
1664  if (.not. is_in_surface_box( contact%master(id), coord(1:3), contact%cparam%BOX_EXP_RATE )) cycle
1665  nmaster = nmaster + 1
1666  indexmaster(nmaster) = id
1667  enddo
1668  deallocate(indexcand)
1669 
1670  if(nmaster == 0) then
1671  deallocate(indexmaster)
1672  cycle
1673  endif
1674 
1675  do idm = 1,nmaster
1676  id = indexmaster(idm)
1677  etype = contact%master(id)%etype
1678  nn = size( contact%master(id)%nodes )
1679  do j=1,nn
1680  iss = contact%master(id)%nodes(j)
1681  elem(1:3,j)=currpos(3*iss-2:3*iss)
1682  enddo
1683  call project_point2element( coord,etype,nn,elem,contact%master(id)%reflen,contact%states(i), &
1684  isin,distclr,localclr=contact%cparam%CLEARANCE )
1685  if( .not. isin ) cycle
1686  contact%states(i)%surface = id
1687  contact%states(i)%multiplier(:) = 0.d0
1688  iss = isinsideelement( etype, contact%states(i)%lpos(1:2), contact%cparam%CLR_CAL_NORM )
1689  if( iss>0 ) &
1690  call cal_node_normal( id, iss, contact%master, currpos, contact%states(i)%lpos(1:2), &
1691  contact%states(i)%direction(:) )
1692  contact_surf(contact%slave(i)) = id
1693  write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3,A,3f7.3)') "Node",nodeid(slave)," contact with element", &
1694  elemid(contact%master(id)%eid), &
1695  " with distance ", contact%states(i)%distance," at ",contact%states(i)%lpos(1:2), &
1696  " along direction ", contact%states(i)%direction
1697  exit
1698  enddo
1699  deallocate(indexmaster)
1700  endif
1701  enddo
1702  !$omp end parallel do
1703 
1704  call hecmw_contact_comm_allreduce_i(contact%comm, contact_surf, hecmw_min)
1705  nactive = 0
1706  do i = 1, size(contact%slave)
1707  if (contact%states(i)%state /= contactfree) then ! any slave in contact
1708  if (abs(contact_surf(contact%slave(i))) /= contact%states(i)%surface) then ! that is in contact with other surface
1709  contact%states(i)%state = contactfree ! should be freed
1710  write(*,'(A,i10,A,i6,A,i6,A)') "Node",nodeid(contact%slave(i)), &
1711  " in rank",hecmw_comm_get_rank()," freed due to duplication"
1712  else
1713  nactive = nactive + 1
1714  endif
1715  endif
1716  if (states_prev(i) == contactfree .and. contact%states(i)%state /= contactfree) then
1717  infoctchange%free2contact = infoctchange%free2contact + 1
1718  elseif (states_prev(i) /= contactfree .and. contact%states(i)%state == contactfree) then
1719  infoctchange%contact2free = infoctchange%contact2free + 1
1720  endif
1721  enddo
1722  active = (nactive > 0)
1723  deallocate(contact_surf)
1724  deallocate(states_prev)
1725  end subroutine scan_contact_state_exp
1726 
1727 end module mcontactdef
mcontact::mu
real(kind=kreal), save mu
penalty, default value
Definition: fstr_contact.f90:17
mcontactdef::tcontact
Structure to includes all info needed by contact calculation.
Definition: fstr_contact_def.F90:26
mcontactdef::update_contact_multiplier
subroutine update_contact_multiplier(contact, coord, disp, ddisp, fcoeff, mu, mut, gnt, ctchanged)
This subroutine update lagrangian multiplier and the distance between contacting nodes.
Definition: fstr_contact_def.F90:1248
mcontactdef::track_contact_position_exp
subroutine track_contact_position_exp(nslave, contact, currpos, currdisp, infoCTChange, nodeID, elemID)
This subroutine tracks down next contact position after a finite slide.
Definition: fstr_contact_def.F90:1473
m_hecmw_contact_comm
Definition: hecmw_comm_contact_f.f90:6
mcontactdef::fstr_contact_finalize
subroutine fstr_contact_finalize(contact)
Finalizer.
Definition: fstr_contact_def.F90:108
m_fstr::myrank
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:96
msurfelement::finalize_surf
subroutine finalize_surf(surf)
Memory management subroutine.
Definition: surf_ele.f90:69
m_contact_lib::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: contact_lib.f90:365
mcontactdef
This module manages the data structure for contact calculation.
Definition: fstr_contact_def.F90:12
mcontactdef::scan_embed_state
subroutine scan_embed_state(flag_ctAlgo, embed, currpos, currdisp, ndforce, infoCTChange, nodeID, elemID, is_init, active, mu, B)
This subroutine update contact states, which include.
Definition: fstr_contact_def.F90:602
mcontactdef::calcu_tied_force0
subroutine calcu_tied_force0(contact, disp, ddisp, mu, B)
This subroutine updates contact condition as follows:
Definition: fstr_contact_def.F90:1200
mcontactdef::ass_contact_force
subroutine ass_contact_force(contact, coord, disp, B)
This subroutine assemble contact force into contacing nodes.
Definition: fstr_contact_def.F90:1381
elementinfo::surfacenormal
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:870
m_contact_lib::contacttied
integer, parameter contacttied
contact type or algorithm definition
Definition: contact_lib.f90:21
mcontactdef::scan_contact_state
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.
Definition: fstr_contact_def.F90:416
bucket_search
This module provides bucket-search functionality It provides definition of bucket info and its access...
Definition: bucket_search.f90:7
msurfelement::update_surface_bucket_info
subroutine update_surface_bucket_info(surf, bktDB)
Update bucket info for searching surface elements.
Definition: surf_ele.f90:280
mcontactparam
This module manage the parameters for contact calculation.
Definition: fstr_contact_param.f90:7
mcontactdef::fstr_contact_init
logical function fstr_contact_init(contact, hecMESH, cparam, myrank)
Initializer of tContactState.
Definition: fstr_contact_def.F90:153
msurfelement::write_surf
subroutine write_surf(file, surf)
Write out elemental surface.
Definition: surf_ele.f90:76
m_hecmw_contact_comm::hecmw_contact_comm_finalize
subroutine, public hecmw_contact_comm_finalize(conComm)
Definition: hecmw_comm_contact_f.f90:160
mcontactdef::update_tied_multiplier
subroutine update_tied_multiplier(contact, disp, ddisp, mu, ctchanged)
This subroutine update lagrangian multiplier and the distance between contacting nodes.
Definition: fstr_contact_def.F90:1332
msurfelement::is_in_surface_box
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:265
m_contact_lib::contactslip
integer, parameter contactslip
Definition: contact_lib.f90:18
mcontactdef::fstr_write_contact
subroutine fstr_write_contact(file, contact)
Write out contact definition.
Definition: fstr_contact_def.F90:86
m_contact_lib::update_tangentforce
subroutine update_tangentforce(etype, nn, elemt0, elemt, cstate)
This subroutine find the projection of a slave point onto master surface.
Definition: contact_lib.f90:452
m_contact_lib::contactfree
integer, parameter contactfree
contact state definition
Definition: contact_lib.f90:16
m_fstr::dt
real(kind=kreal) dt
ANALYSIS CONTROL for NLGEOM and HEAT.
Definition: m_fstr.f90:139
mcontactdef::fstr_info_contactchange
Definition: fstr_contact_def.F90:59
msurfelement::tsurfelement
Structure to define surface group.
Definition: surf_ele.f90:19
msurfelement::update_surface_box_info
subroutine update_surface_box_info(surf, currpos)
Update info of cubic box including surface elements.
Definition: surf_ele.f90:240
mcontactdef::scan_contact_state_exp
subroutine scan_contact_state_exp(contact, currpos, currdisp, infoCTChange, nodeID, elemID, is_init, active)
This subroutine update contact states, which include.
Definition: fstr_contact_def.F90:1585
mcontact::gnt
real(kind=kreal), dimension(2), save gnt
1:current average penetration; 2:current relative tangent displacement
Definition: fstr_contact.f90:24
msurfelement::find_surface_neighbor
subroutine find_surface_neighbor(surf, bktDB)
Find neighboring surface elements.
Definition: surf_ele.f90:88
mcontactdef::fstr_embed_init
logical function fstr_embed_init(embed, hecMESH, cparam, myrank)
Initializer of tContactState for embed case.
Definition: fstr_contact_def.F90:272
m_contact_lib::contactstick
integer, parameter contactstick
Definition: contact_lib.f90:17
msurfelement::update_surface_reflen
subroutine update_surface_reflen(surf, coord)
Compute reference length of surface elements.
Definition: surf_ele.f90:220
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
bucket_search::bucketdb_getnumcand
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...
Definition: bucket_search.f90:313
bucket_search::bucketdb_getcand
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...
Definition: bucket_search.f90:344
hecmw
Definition: hecmw.f90:6
mcontactdef::update_contact_tangentforce
subroutine update_contact_tangentforce(contact)
Definition: fstr_contact_def.F90:1455
bucket_search::bucketdb_finalize
subroutine, public bucketdb_finalize(bktdb)
Finalizer.
Definition: bucket_search.f90:151
m_contact_lib
This module provide functions of contact stiffness calculation.
Definition: contact_lib.f90:6
mcontactdef::set_contact_state_vector
subroutine set_contact_state_vector(contact, dt, relvel_vec, state_vec)
This subroutine setup contact output nodal vectors.
Definition: fstr_contact_def.F90:1435
m_hecmw_contact_comm::hecmw_contact_comm_allreduce_i
subroutine, public hecmw_contact_comm_allreduce_i(conComm, vec, op)
Definition: hecmw_comm_contact_f.f90:224
bucket_search::bucketdb_getbucketid
integer(kind=kint) function, public bucketdb_getbucketid(bktdb, x)
Get bucket ID that includes given point.
Definition: bucket_search.f90:244
m_contact_lib::contact_state_init
subroutine contact_state_init(cstate)
Initializer.
Definition: contact_lib.f90:48
elementinfo::getvertexcoord
subroutine getvertexcoord(fetype, cnode, localcoord)
Get the natural coord of a vertex node.
Definition: element.f90:1136
m_contact_lib::dispincrematrix
subroutine dispincrematrix(pos, etype, nnode, ele, tangent, tensor, matrix)
This subroutine calculate the relation between global disp and displacement along natural coordinate ...
Definition: contact_lib.f90:216
mcontactdef::update_direction
subroutine update_direction(nslave, contact, currpos)
This subroutine tracks down next contact position after a finite slide.
Definition: fstr_contact_def.F90:978
mcontactdef::clear_contact_state
subroutine clear_contact_state(contact)
Reset contact state all to free.
Definition: fstr_contact_def.F90:392
mcontactdef::calcu_contact_force0
subroutine calcu_contact_force0(contact, coord, disp, ddisp, fcoeff, mu, mut, B)
This subroutine updates contact condition as follows:
Definition: fstr_contact_def.F90:1126
m_hecmw_contact_comm::hecmw_contact_comm_init
subroutine, public hecmw_contact_comm_init(conComm, hecMESH, ndof, n_contact_dof, contact_dofs)
Definition: hecmw_comm_contact_f.f90:42
mcontact::mut
real(kind=kreal), save mut
penalty along tangent direction
Definition: fstr_contact.f90:18
m_contact_lib::project_point2element
subroutine project_point2element(xyz, etype, nn, elemt, reflen, cstate, isin, distclr, ctpos, localclr)
This subroutine find the projection of a slave point onto master surface.
Definition: contact_lib.f90:257
bucket_search::bucketdb_init
subroutine, public bucketdb_init(bktdb)
Initializer.
Definition: bucket_search.f90:138
mcontactdef::fstr_contact_check
logical function fstr_contact_check(contact, hecMESH)
Check the consistency with given mesh of contact definition.
Definition: fstr_contact_def.F90:124
msurfelement::initialize_surf
subroutine initialize_surf(eid, etype, nsurf, surf)
Initializer.
Definition: surf_ele.f90:38
m_contact_lib::contactfslid
integer, parameter contactfslid
Definition: contact_lib.f90:24
msurfelement
This module manages surface elements in 3D It provides basic definition of surface elements (triangla...
Definition: surf_ele.f90:8
m_contact_lib::tcontactstate
This structure records contact status.
Definition: contact_lib.f90:27