22 include
'fstr_ctrl_util_f.inc'
27 character(len=HECMW_NAME_LEN) :: name
30 character(len=HECMW_NAME_LEN) :: pair_name
31 integer :: surf_id1, surf_id2
32 integer :: surf_id1_sgrp
34 integer,
pointer :: slave(:)=>null()
35 real(kind=kreal) :: fcoeff
36 real(kind=kreal) :: tpenalty
52 type(fstrst_contact_comm) :: comm
53 type(bucketdb) :: master_bktdb
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
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
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
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
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
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.
103 integer(kind=kint),
intent(in) :: file
104 type(
tcontact),
intent(in) :: contact
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)
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)
125 type(
tcontact),
intent(inout) :: contact
127 if(
associated( contact%slave ) )
deallocate(contact%slave)
128 if(
associated( contact%master ) )
then
129 do i=1,
size( contact%master )
132 deallocate(contact%master)
134 if(
associated(contact%states) )
deallocate(contact%states)
141 type(
tcontact),
intent(inout) :: contact
142 type(hecmwst_local_mesh),
pointer :: hecmesh
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)
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
170 type(
tcontact),
intent(inout) :: contact
171 type(hecmwst_local_mesh),
pointer :: hecmesh
172 integer(kint),
intent(in),
optional :: myrank
174 integer :: i, j, is, ie, cgrp, nsurf, nslave, ic, ic_type, iss, nn, ii
175 integer :: count,nodeid
179 cgrp = contact%surf_id2
181 is= hecmesh%surf_group%grp_index(cgrp-1) + 1
182 ie= hecmesh%surf_group%grp_index(cgrp )
184 if(
present(myrank))
then
188 ic = hecmesh%surf_group%grp_item(2*i-1)
189 if(hecmesh%elem_ID(ic*2) /= myrank) cycle
192 allocate( contact%master(count) )
195 ic = hecmesh%surf_group%grp_item(2*i-1)
196 if(hecmesh%elem_ID(ic*2) /= myrank) cycle
198 nsurf = hecmesh%surf_group%grp_item(2*i)
199 ic_type = hecmesh%elem_type(ic)
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 )
210 allocate( contact%master(ie-is+1) )
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)
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 )
229 cgrp = contact%surf_id1
231 is= hecmesh%node_group%grp_index(cgrp-1) + 1
232 ie= hecmesh%node_group%grp_index(cgrp )
235 nodeid = hecmesh%global_node_ID(hecmesh%node_group%grp_item(i))
236 if(
present(myrank))
then
241 if( hecmesh%node_group%grp_item(i) <= hecmesh%nn_internal)
then
246 allocate( contact%slave(nslave) )
247 allocate( contact%states(nslave) )
250 if(.not.
present(myrank))
then
252 if( hecmesh%node_group%grp_item(i) > hecmesh%nn_internal) cycle
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
281 contact%symmetric = .true.
287 type(
tcontact),
intent(inout) :: contact
289 if( .not.
associated(contact%states) )
return
290 do i=1,
size( contact%states )
291 contact%states(i)%state = -1
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.
302 is_active_contact = .false.
310 nodeID, elemID, is_init, active, mu, B )
311 character(len=9),
intent(in) :: flag_ctAlgo
312 type(
tcontact ),
intent(inout) :: contact
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(:)
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)
330 integer(kind=kint),
allocatable :: contact_surf(:), states_prev(:)
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(:)
337 if( contact%algtype<=2 )
return
342 distclr = distclr_free
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
356 is_present_b =
present(b)
357 if( is_present_b ) bp => b
367 do i= 1,
size(contact%slave)
368 slave = contact%slave(i)
370 slforce(1:3)=ndforce(3*slave-2:3*slave)
371 id = contact%states(i)%surface
372 nlforce = contact%states(i)%multiplier(1)
374 if( nlforce < tensile_force )
then
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
381 if( contact%algtype /=
contactfslid .or. (.not. is_present_b) )
then
382 contact_surf(contact%slave(i)) = -id
384 call track_contact_position( flag_ctalgo, i, contact, currpos, currdisp, mu, infoctchange, nodeid, elemid, bp )
386 contact_surf(contact%slave(i)) = -contact%states(i)%surface
390 else if( contact%states(i)%state==
contactfree )
then
391 coord(:) = currpos(3*slave-2:3*slave)
396 if (ncand == 0) cycle
397 allocate(indexcand(ncand))
401 allocate(indexmaster(nmastermax))
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
411 deallocate(indexcand)
413 if(nmaster == 0)
then
414 deallocate(indexmaster)
419 id = indexmaster(idm)
420 etype = contact%master(id)%etype
421 nn =
size( contact%master(id)%nodes )
423 iss = contact%master(id)%nodes(j)
424 elem(1:3,j)=currpos(3*iss-2:3*iss)
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 )
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
442 deallocate(indexmaster)
449 do i = 1,
size(contact%slave)
451 if (abs(contact_surf(contact%slave(i))) /= contact%states(i)%surface)
then
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"
456 nactive = nactive + 1
460 infoctchange%free2contact = infoctchange%free2contact + 1
462 infoctchange%contact2free = infoctchange%contact2free + 1
465 active = (nactive > 0)
466 deallocate(contact_surf)
467 deallocate(states_prev)
471 subroutine cal_node_normal( csurf, isin, surf, currpos, lpos, normal )
473 integer,
intent(in) :: csurf
474 integer,
intent(in) :: isin
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)
484 if( 1 <= isin .and. isin <= 4 )
then
486 gn = surf(csurf)%nodes(cnode)
487 etype = surf(csurf)%etype
489 nn =
size( surf(csurf)%nodes )
491 iss = surf(csurf)%nodes(j)
492 elem(1:3,j)=currpos(3*iss-2:3*iss)
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
502 iss = surf(nd1)%nodes(j)
503 elem(1:3,j)=currpos(3*iss-2:3*iss)
510 normal = normal+normal_n
515 elseif( 12 <= isin .and. isin <= 41 )
then
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 )
523 iss = surf(csurf)%nodes(j)
524 elem(1:3,j)=currpos(3*iss-2:3*iss)
528 case (fe_tri3n, fe_tri6n, fe_tri6nc)
529 if ( isin==12 ) then; x=lpos(2)-lpos(1)
530 elseif( isin==23 ) then; x=1.d0-2.d0*lpos(2)
531 elseif( isin==31 ) then; x=2.d0*lpos(1)-1.d0
532 else; stop
"Error: cal_node_normal: invalid isin"
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"
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
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
556 if( cgn1>0 .and. cgn2>0 )
then
558 isin_n = 10*cgn2 + cgn1
561 case (fe_tri3n, fe_tri6n, fe_tri6nc)
562 if ( isin_n==12 ) then; lpos_n(1)=0.5d0*(1.d0-x); lpos_n(2)=0.5d0*(1.d0+x)
563 elseif( isin_n==23 ) then; lpos_n(1)=0.d0; lpos_n(2)=0.5d0*(1.d0-x)
564 elseif( isin_n==31 ) then; lpos_n(1)=0.5d0*(1.d0+x); lpos_n(2)=0.d0
565 else; stop
"Error: cal_node_normal: invalid isin_n"
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"
577 normal = normal+normal_n
584 normal = normal/ dsqrt( dot_product( normal, normal ) )
585 end subroutine cal_node_normal
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
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(:)
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 )
604 real(kind=kreal) :: opos(2), odirec(3)
605 integer(kind=kint) :: bktID, nCand, idm
606 integer(kind=kint),
allocatable :: indexCand(:)
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 )
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)
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 )
631 iss = contact%master(sid)%nodes(j)
632 elem(1:3,j)=currpos(3*iss-2:3*iss)
635 isin,distclr_nocheck,localclr=clearance )
637 contact%states(nslave)%surface = sid
643 if( .not. isin )
then
644 write(*,*)
'Warning: contact moved beyond neighbor elements'
649 allocate(indexcand(ncand))
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 )
659 iss = contact%master(sid)%nodes(j)
660 elem(1:3,j)=currpos(3*iss-2:3*iss)
663 isin,distclr_nocheck,localclr=clearance )
665 contact%states(nslave)%surface = sid
669 deallocate(indexcand)
674 if( contact%states(nslave)%surface==sid0 )
then
675 if(any(dabs(contact%states(nslave)%lpos(:)-opos(:)) >= clr_difflpos))
then
677 infoctchange%contact2difflpos = infoctchange%contact2difflpos + 1
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(:)
684 infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
685 if( flag_ctalgo==
'ALagrange' ) &
686 call reset_contact_force( contact, currpos, nslave, sid0, opos, odirec, b )
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 )
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"
696 contact%states(nslave)%multiplier(:) = 0.d0
699 end subroutine track_contact_position
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(:)
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
721 slave = contact%slave(lslave)
722 fcoeff = contact%fcoeff
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 )
730 iss = contact%master(omaster)%nodes(j)
735 b(idx0+i) = b(idx0+i)-nrlforce*shapefunc(j)*odirec(i)
738 if( fcoeff/=0.d0 )
then
740 iss = contact%master(omaster)%nodes(j)
741 elemcrd(:,j) = currpos(3*iss-2:3*iss)
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)
749 iss = contact%master(omaster)%nodes(j)
754 b(idx0+i) = b(idx0+i)+f3(3*j+i)
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
766 iss = contact%master(master)%nodes(j)
772 b(idx0+i) = b(idx0+i)+nrlforce* &
773 shapefunc(j)*contact%states(lslave)%direction(i)
776 if( fcoeff/=0.d0 )
then
778 iss = contact%master(master)%nodes(j)
779 elemcrd(:,j) = currpos(3*iss-2:3*iss)
781 call dispincrematrix( contact%states(lslave)%lpos, etype, nn, elemcrd, tangent, &
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)
787 iss = contact%master(master)%nodes(j)
792 b(idx0+i) = b(idx0+i)-f3(3*j+i)
797 end subroutine reset_contact_force
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(:)
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)
820 do i= 1,
size(contact%slave)
822 slave = contact%slave(i)
823 edisp(1:3) = ddisp(3*slave-2:3*slave)
824 master = contact%states(i)%surface
826 nn =
size( contact%master(master)%nodes )
827 etype = contact%master(master)%etype
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)
834 call getshapefunc( etype, contact%states(i)%lpos(:), shapefunc )
839 elemg(:) = elemg(:)+shapefunc(j)*elemdisp(:,j)
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
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
852 if( fcoeff==0.d0 ) cycle
854 call dispincrematrix( contact%states(i)%lpos, etype, nn, elemcrd, tangent, &
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,:)
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
866 b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-f3(1:3)
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)
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
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)
897 do i= 1,
size(contact%slave)
899 slave = contact%slave(i)
900 edisp(1:3) = ddisp(3*slave-2:3*slave)
901 master = contact%states(i)%surface
903 nn =
size( contact%master(master)%nodes )
904 etype = contact%master(master)%etype
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)
911 call getshapefunc( etype, contact%states(i)%lpos(:), shapefunc )
916 elemg(:) = elemg(:)+shapefunc(j)*elemdisp(:,j)
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
924 lgnt(1) = lgnt(1)- contact%states(i)%wkdist
926 if( fcoeff==0.d0 ) cycle
928 call dispincrematrix( contact%states(i)%lpos, etype, nn, elemcrd, tangent, &
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
940 print *,
"Node", slave,
"to slip state",dgn, fcoeff*contact%states(i)%multiplier(1)
943 fric(:) = fric(:)*fcoeff*contact%states(i)%multiplier(1)/dgn
947 print *,
"Node", slave,
"to stick state",dgn, fcoeff*contact%states(i)%multiplier(1)
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) )
956 if(cnt>0) lgnt(:) = lgnt(:)/cnt
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(:)
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)
977 slave = contact%slave(i)
978 master = contact%states(i)%surface
980 nn =
size( contact%master(master)%nodes )
981 etype = contact%master(master)%etype
982 call getshapefunc( etype, contact%states(i)%lpos(:), shapefunc )
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
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
994 if( fcoeff==0.d0 ) cycle
997 iss = contact%master(master)%nodes(j)
998 elemcrd(:,j) = coord(3*iss-2:3*iss)+disp(3*iss-2:3*iss)
1000 call dispincrematrix( contact%states(i)%lpos, etype, nn, elemcrd, tangent, &
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)
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)
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(:)
1021 integer(kind=kint) :: i, slave
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)
1029 if( dt < 1.d-16 ) cycle
1030 relvel_vec(3*slave-2:3*slave) = contact%states(i)%reldisp(1:3)/dt
1036 type(
tcontact ),
intent(inout) :: contact
1038 integer(kind=kint) :: i
1040 do i= 1,
size(contact%slave)
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
1046 contact%states(i)%tangentForce(1:3) = contact%states(i)%tangentForce_final(1:3)
1048 contact%states(i)%tangentForce1(1:3) = contact%states(i)%tangentForce(1:3)
1054 integer,
intent(in) :: nslave
1055 type(
tcontact ),
intent(inout) :: contact
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(:)
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 )
1066 real(kind=kreal) :: opos(2), odirec(3)
1067 integer(kind=kint) :: bktID, nCand, idm
1068 integer(kind=kint),
allocatable :: indexCand(:)
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 )
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)
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 )
1093 iss = contact%master(sid)%nodes(j)
1094 elem(1:3,j)=currpos(3*iss-2:3*iss)
1096 call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
1097 isin,distclr_nocheck,localclr=clearance )
1099 contact%states(nslave)%surface = sid
1105 if( .not. isin )
then
1106 write(*,*)
'Warning: contact moved beyond neighbor elements'
1111 allocate(indexcand(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 )
1121 iss = contact%master(sid)%nodes(j)
1122 elem(1:3,j)=currpos(3*iss-2:3*iss)
1124 call project_point2element( coord,etype,nn,elem,contact%master(sid)%reflen,contact%states(nslave), &
1125 isin,distclr_nocheck,localclr=clearance )
1127 contact%states(nslave)%surface = sid
1131 deallocate(indexcand)
1136 if( contact%states(nslave)%surface==sid0 )
then
1137 if(any(dabs(contact%states(nslave)%lpos(:)-opos(:)) >= clr_difflpos))
then
1139 infoctchange%contact2difflpos = infoctchange%contact2difflpos + 1
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(:)
1146 infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
1148 iss = isinsideelement( etype, contact%states(nslave)%lpos, clr_cal_norm )
1150 call cal_node_normal( contact%states(nslave)%surface, iss, contact%master, currpos, &
1151 contact%states(nslave)%lpos, contact%states(nslave)%direction(:) )
1153 else if( .not. isin )
then
1154 write(*,
'(A,i10,A)')
"Node",nodeid(slave),
" move out of contact"
1156 contact%states(nslave)%multiplier(:) = 0.d0
1165 nodeID, elemID, is_init, active )
1166 type(
tcontact ),
intent(inout) :: contact
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
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
1181 integer(kind=kint),
allocatable :: contact_surf(:), states_prev(:)
1183 integer,
pointer :: indexMaster(:),indexCand(:)
1184 integer :: nMaster,idm,nMasterMax,bktID,nCand
1190 distclr = distclr_free
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
1211 do i= 1,
size(contact%slave)
1212 slave = contact%slave(i)
1216 contact_surf(contact%slave(i)) = -contact%states(i)%surface
1218 else if( contact%states(i)%state==
contactfree )
then
1219 coord(:) = currpos(3*slave-2:3*slave)
1224 if (ncand == 0) cycle
1225 allocate(indexcand(ncand))
1229 allocate(indexmaster(nmastermax))
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
1239 deallocate(indexcand)
1241 if(nmaster == 0)
then
1242 deallocate(indexmaster)
1247 id = indexmaster(idm)
1248 etype = contact%master(id)%etype
1249 nn =
size( contact%master(id)%nodes )
1251 iss = contact%master(id)%nodes(j)
1252 elem(1:3,j)=currpos(3*iss-2:3*iss)
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 )
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
1270 deallocate(indexmaster)
1277 do i = 1,
size(contact%slave)
1279 if (abs(contact_surf(contact%slave(i))) /= contact%states(i)%surface)
then
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"
1284 nactive = nactive + 1
1288 infoctchange%free2contact = infoctchange%free2contact + 1
1290 infoctchange%contact2free = infoctchange%contact2free + 1
1293 active = (nactive > 0)
1294 deallocate(contact_surf)
1295 deallocate(states_prev)
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.
subroutine getvertexcoord(fetype, cnode, localcoord)
Get the natural coord of a vertex node.
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
This module manage surface elements in 3D It provides basic definition of surface elements (trianglar...
subroutine initialize_surf(eid, etype, nsurf, surf)
Initializer.
subroutine update_surface_reflen(surf, coord)
Compute reference length of surface elements.
subroutine write_surf(file, surf)
Write out elemental surface.
subroutine update_surface_box_info(surf, currpos)
Update info of cubic box including surface elements.
subroutine find_surface_neighbor(surf, bktDB)
Find neighboring surface elements.
subroutine update_surface_bucket_info(surf, bktDB)
Update bucket info for searching surface elements.
subroutine finalize_surf(surf)
Memeory management subroutine.
logical function is_in_surface_box(surf, x0, exp_rate)
Check if given point is within cubic box including surface element.
Structure to define surface group.