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
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)
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"
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
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"
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.
integer, parameter fe_tri6n
subroutine getvertexcoord(fetype, cnode, localcoord)
Get the natural coord of a vertex node.
integer, parameter fe_tri3n
integer, parameter fe_quad4n
integer, parameter fe_tri6nc
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
integer, parameter fe_quad8n
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 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.
subroutine update_surface_bucket_info(surf, bktdb)
Update bucket info for searching surface elements.
Structure for bucket search.
Structure to define surface group.