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