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