13 private :: l_contact2mpc, l_tied2mpc
15 logical,
private :: active
17 real(kind=kreal),
save ::
mu=1.d10
18 real(kind=kreal),
save ::
mut=1.d6
19 real(kind=kreal),
save ::
cdotp=1.d3
21 real(kind=kreal),
save ::
cgn=1.d-5
22 real(kind=kreal),
save ::
cgt=1.d-3
24 real(kind=kreal),
save ::
gnt(2)
33 integer(kind=kint),
intent(in) :: file
34 type( hecmwst_contact_pair ),
intent(in) :: pair
36 integer(kind=kint) :: i
37 write(file,*)
"Number of contact pair", pair%n_pair
39 write(file,*) trim(pair%name(i)), pair%type(i), pair%slave_grp_id(i) &
40 ,pair%master_grp_id(i), pair%slave_orisgrp_id(i)
45 real(kind=kreal),
intent(in) :: maxv
56 logical,
intent(in) :: a
61 integer(kind=kint),
intent(in) :: ctalgo
63 type (hecmwst_local_mesh),
intent(in) :: hecmesh
64 integer(kind=kint) :: is_conv
67 if( infoctchange%contact2free+infoctchange%contact2neighbor+ &
68 infoctchange%contact2difflpos+infoctchange%free2contact == 0 ) &
72 write(*,
'(a,2e15.7)')
"--Relative displacement in contact surface",
gnt
79 call hecmw_allreduce_i1(hecmesh, is_conv, hecmw_min)
80 if (is_conv == 0)
then
90 if( infoctchange%contact2free+infoctchange%contact2neighbor+infoctchange%free2contact > 0 ) &
95 subroutine l_contact2mpc( contact, mpcs, nmpc )
97 type(
tcontact ),
intent(in) :: contact
98 type( hecmwst_mpc ),
intent(inout) :: mpcs
99 integer(kind=kint),
intent(out) :: nmpc
100 integer(kind=kint),
parameter :: ndof = 3
101 real(kind=kreal),
parameter :: tol =1.d-10
102 integer(kind=kint) :: i, j, k, nn, csurf, nenode, etype, tdof
103 integer(kind=kint) :: nodes(l_max_surface_node*ndof+ndof), dofs(l_max_surface_node*ndof+ndof)
104 real(kind=kreal) :: values(l_max_surface_node*ndof+ndof+1),val(l_max_surface_node*ndof+ndof+1)
106 do i=1,
size(contact%states)
107 if( contact%states(i)%state == -1 ) cycle
108 csurf = contact%states(i)%surface
109 if( csurf<=0 ) stop
"error in contact state"
110 etype = contact%master(csurf)%etype
111 nenode =
size(contact%master(csurf)%nodes)
112 tdof = nenode*ndof+ndof
113 call contact2mpcval( contact%states(i), etype, nenode, values(1:tdof+1) )
116 if( dabs(values(j))<tol ) cycle
118 nodes(tdof) = contact%slave(i)
120 val(tdof) = values(j)
123 nn = contact%master(csurf)%nodes(j)
124 nodes( j*ndof+1:j*ndof+ndof ) = nn
126 if( dabs(values(j*ndof+k)) < tol ) cycle
130 val(tdof)=values(j*ndof+k)
133 val(tdof+1) = values(nenode*ndof+ndof+1)
135 call fstr_append_mpc( tdof, nodes(1:tdof), dofs(1:tdof), val(1:tdof+1), mpcs )
138 end subroutine l_contact2mpc
141 subroutine l_tied2mpc( contact, mpcs, nmpc )
143 type(
tcontact ),
intent(in) :: contact
144 type( hecmwst_mpc ),
intent(inout) :: mpcs
145 integer(kind=kint),
intent(out) :: nmpc
146 integer(kind=kint) :: i, j, csurf, nenode, etype, tdof
147 integer(kind=kint) :: nodes(l_max_surface_node+1), dofs(l_max_surface_node+1)
148 real(kind=kreal) :: values(l_max_surface_node+2)
150 do i=1,
size(contact%slave)
151 csurf = contact%states(i)%surface
153 nenode =
size(contact%master(csurf)%nodes)
155 nodes(1) = contact%slave(i)
156 nodes( 2:tdof ) = contact%master(csurf)%nodes(:)
158 values(2:tdof) = 1.d0
159 values(tdof+1) = 0.d0
160 etype = contact%master(csurf)%etype
163 call fstr_append_mpc( tdof, nodes(1:tdof), dofs(1:tdof), values(1:tdof+1), mpcs )
167 end subroutine l_tied2mpc
171 type(
tcontact ),
intent(in) :: contacts(:)
172 type( hecmwst_mpc ),
intent(inout) :: mpcs
173 integer(kind=kint) :: i, nmpc
175 do i=1,
size(contacts)
176 if( contacts(i)%algtype == contactunknown ) cycle
177 if( contacts(i)%algtype == contactfslid )
then
178 print *,
"Cannot deal with finit slip problems by MPC!"
181 if( contacts(i)%algtype == contactsslid )
then
182 call l_contact2mpc( contacts(i), mpcs, nmpc )
184 elseif( contacts(i)%algtype == contacttied )
then
185 call l_tied2mpc( contacts(i), mpcs, nmpc )
194 type( hecmwst_mpc ),
intent(inout) :: mpcs
200 integer(kind=kint),
intent(in) :: file
201 type( hecmwst_mpc ),
intent(in) :: mpcs
203 integer(kind=kint) :: i,j,n0,n1
204 write(file, *)
"Number of equation", mpcs%n_mpc
206 write(file,*)
"--Equation",i
207 n0=mpcs%mpc_index(i-1)+1
210 write(file,
'(30i5)') (mpcs%mpc_item(j),j=n0,n1)
211 write(file,
'(30i5)') (mpcs%mpc_dof(j),j=n0,n1)
212 write(file,
'(30f7.2)') (mpcs%mpc_val(j),j=n0,n1),mpcs%mpc_const(i)
218 integer(kind=kint),
intent(in) :: cstep
219 integer(kind=kint),
intent(in) :: sub_step
220 integer(kind=kint),
intent(in) :: cont_step
221 real(kind=kreal),
intent(in) ::
dt
222 integer(kind=kint),
intent(in) :: ctAlgo
223 type( hecmwst_local_mesh ),
intent(in) :: hecMESH
227 real(kind=kreal),
optional :: b(:)
228 character(len=9) :: flag_ctAlgo
229 integer(kind=kint) :: i, grpid
230 logical :: iactive, is_init
232 fstrsolid%CONT_RELVEL(:) = 0.d0
233 fstrsolid%CONT_STATE(:) = 0.d0
236 flag_ctalgo =
'SLagrange'
238 flag_ctalgo =
'ALagrange'
243 do i = 1,
size(fstrsolid%unode)
244 fstrsolid%ddunode(i) = hecmesh%node(i) + fstrsolid%unode(i) + fstrsolid%dunode(i)
248 infoctchange%contact2free = 0
249 infoctchange%contact2neighbor = 0
250 infoctchange%contact2diffLpos = 0
251 infoctchange%free2contact = 0
252 infoctchange%contactNode_current = 0
254 is_init = ( cstep == 1 .and. sub_step == 1 .and. cont_step == 0 )
256 do i=1,
size(fstrsolid%contacts)
257 grpid = fstrsolid%contacts(i)%group
261 if(
present(b) )
then
262 call scan_contact_state( flag_ctalgo, fstrsolid%contacts(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
263 & fstrsolid%QFORCE(:), infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive,
mu, b )
265 call scan_contact_state( flag_ctalgo, fstrsolid%contacts(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
266 & fstrsolid%QFORCE(:), infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive,
mu )
268 if( .not. active ) active = iactive
274 infoctchange%contactNode_current = infoctchange%contactNode_previous+infoctchange%free2contact-infoctchange%contact2free
275 infoctchange%contactNode_previous = infoctchange%contactNode_current
277 if( .not. active )
then
278 fstrsolid%CONT_NFORCE(:) = 0.d0
279 fstrsolid%CONT_FRIC(:) = 0.d0
286 integer(kind=kint),
intent(in) :: cstep
287 type( hecmwst_local_mesh ),
intent(in) :: hecMESH
291 integer(kind=kint) :: i
292 logical :: iactive, is_init
297 do i = 1,
size(fstrsolid%unode)
298 fstrsolid%ddunode(i) = hecmesh%node(i) + fstrsolid%unode(i) + fstrsolid%dunode(i)
300 infoctchange%active = .false.
302 infoctchange%contact2free = 0
303 infoctchange%contact2neighbor = 0
304 infoctchange%contact2diffLpos = 0
305 infoctchange%free2contact = 0
306 infoctchange%contactNode_current = 0
308 is_init = ( cstep == 1 )
310 do i=1,
size(fstrsolid%contacts)
317 & infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive )
319 infoctchange%active = infoctchange%active .or. iactive
322 infoctchange%contactNode_current = infoctchange%contactNode_previous+infoctchange%free2contact-infoctchange%contact2free
323 infoctchange%contactNode_previous = infoctchange%contactNode_current
324 fstrsolid%ddunode = 0.d0
329 type( hecmwst_local_mesh ),
intent(in) :: hecMESH
331 real(kind=kreal),
intent(inout) :: b(:)
333 integer(kind=kint) :: i
334 do i=1,
size(fstrsolid%contacts)
337 , fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff,
mu,
mut, b )
343 type( hecmwst_local_mesh ),
intent(in) :: hecMESH
345 logical,
intent(out) :: ctchanged
347 integer(kind=kint) :: i, nc
348 gnt = 0.d0; ctchanged = .false.
349 nc =
size(fstrsolid%contacts)
353 , fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff,
mu,
mut,
gnt, ctchanged )
362 integer(kind=kint) :: i, nc
364 nc =
size(fstrsolid%contacts)
373 integer(kind=kint),
intent(in) :: iter
374 type (hecmwST_local_mesh),
intent(inout) :: hecMESH
375 type (hecmwST_matrix),
intent(inout) :: hecMAT
378 integer(kind=kint),
parameter :: NDOF=3
380 integer(kind=kint) :: i, j, k, m, nnode, nd, etype
381 integer(kind=kint) :: ctsurf, ndLocal(l_max_surface_node+1)
382 real(kind=kreal) :: factor, elecoord( 3, l_max_surface_node)
383 real(kind=kreal) :: stiff(l_max_surface_node*3+3, l_max_surface_node*3+3)
384 real(kind=kreal) :: nrlforce, force(l_max_surface_node*3+3)
393 factor = fstrsolid%FACTOR(2)
394 call hecmw_cmat_clear( hecmat%cmat )
395 do i=1,
size(fstrsolid%contacts)
396 do j=1,
size(fstrsolid%contacts(i)%slave)
397 if( fstrsolid%contacts(i)%states(j)%state==contactfree ) cycle
398 ctsurf = fstrsolid%contacts(i)%states(j)%surface
399 etype = fstrsolid%contacts(i)%master(ctsurf)%etype
400 nnode =
size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
401 ndlocal(1) = fstrsolid%contacts(i)%slave(j)
403 ndlocal(k+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(k)
404 elecoord(1,k)=hecmesh%node(3*ndlocal(k+1)-2)
405 elecoord(2,k)=hecmesh%node(3*ndlocal(k+1)-1)
406 elecoord(3,k)=hecmesh%node(3*ndlocal(k+1))
408 call contact2stiff( fstrsolid%contacts(i)%algtype, fstrsolid%contacts(i)%states(j), &
409 etype, nnode, elecoord(:,:),
mu,
mut, fstrsolid%contacts(i)%fcoeff, &
410 fstrsolid%contacts(i)%symmetric, stiff(:,:), force(:) )
411 call hecmw_mat_ass_contact( hecmat,nnode+1,ndlocal,stiff )
416 fstrsolid%contacts(i)%states(j)%wkdist = fstrsolid%contacts(i)%states(j)%distance
417 nrlforce = -
mu*fstrsolid%contacts(i)%states(j)%distance
418 force(1:nnode*ndof+ndof) = force(1:nnode*ndof+ndof)*nrlforce
422 hecmat%B(ndof*(nd-1)+k)=hecmat%B(ndof*(nd-1)+k)-force((m-1)*ndof+k)
433 type(hecmwst_matrix) :: hecMAT
435 if( .not.
associated(fstrsolid%CONT_NFORCE) )
then
436 allocate( fstrsolid%CONT_NFORCE(hecmat%NP*3) )
437 fstrsolid%CONT_NFORCE(:) = 0.d0
440 if( .not.
associated(fstrsolid%CONT_FRIC) )
then
441 allocate( fstrsolid%CONT_FRIC(hecmat%NP*3) )
442 fstrsolid%CONT_FRIC(:) = 0.d0
445 if( .not.
associated(fstrsolid%CONT_RELVEL) )
then
446 allocate( fstrsolid%CONT_RELVEL(hecmat%NP*3) )
447 fstrsolid%CONT_RELVEL(:) = 0.d0
450 if( .not.
associated(fstrsolid%CONT_STATE) )
then
451 allocate( fstrsolid%CONT_STATE(hecmat%NP*1) )
452 fstrsolid%CONT_STATE(:) = 0.d0
455 if( .not.
associated(fstrsolid%CONT_AREA) )
then
456 allocate( fstrsolid%CONT_AREA(hecmat%NP) )
457 fstrsolid%CONT_AREA(:) = 0.d0
460 if( .not.
associated(fstrsolid%CONT_NTRAC) )
then
461 allocate( fstrsolid%CONT_NTRAC(hecmat%NP*3) )
462 fstrsolid%CONT_NTRAC(:) = 0.d0
465 if( .not.
associated(fstrsolid%CONT_FTRAC) )
then
466 allocate( fstrsolid%CONT_FTRAC(hecmat%NP*3) )
467 fstrsolid%CONT_FTRAC(:) = 0.d0
473 integer(kind=kint),
intent(in) :: cstep
474 type( hecmwst_local_mesh ),
intent(in) :: hecMESH
477 integer(kind=kint) :: i, j, k, sgrp_id, iS, iE, eid, sid, n_cdsurfs
478 logical,
pointer :: cdef_surf(:,:)
479 real(kind=kreal),
pointer :: coord(:)
481 if(
associated(fstrsolid%CONT_SGRP_ID) )
deallocate(fstrsolid%CONT_SGRP_ID)
483 allocate(cdef_surf(l_max_elem_surf,hecmesh%n_elem))
484 cdef_surf(:,:) = .false.
488 do i=1,
size(fstrsolid%contacts)
496 sgrp_id = fstrsolid%contacts(i)%surf_id1_sgrp
498 sgrp_id = fstrsolid%contacts(i)%surf_id2
501 if( sgrp_id <= 0 ) cycle
503 is = hecmesh%surf_group%grp_index(sgrp_id-1) + 1
504 ie = hecmesh%surf_group%grp_index(sgrp_id )
506 eid = hecmesh%surf_group%grp_item(2*j-1)
507 sid = hecmesh%surf_group%grp_item(2*j)
509 if( .not. cdef_surf(sid,eid) ) n_cdsurfs = n_cdsurfs + 1
510 cdef_surf(sid,eid) = .true.
516 allocate(fstrsolid%CONT_SGRP_ID(2*n_cdsurfs))
518 do i=1,hecmesh%n_elem
519 do j=1,l_max_elem_surf
520 if( cdef_surf(j,i) )
then
521 n_cdsurfs = n_cdsurfs + 1
522 fstrsolid%CONT_SGRP_ID(2*n_cdsurfs-1) = i
523 fstrsolid%CONT_SGRP_ID(2*n_cdsurfs ) = j
527 deallocate(cdef_surf)
532 type( hecmwst_local_mesh ),
intent(in) :: hecmesh
534 integer(kind=kint),
intent(in) :: flag
536 integer(kind=kint),
parameter :: NDOF=3
537 integer(kind=kint) :: i, isuf, icel, sid, etype, nn, iS, stype, idx
538 integer(kind=kint) :: ndlocal(l_max_elem_node)
539 real(kind=kreal),
pointer :: coord(:)
540 real(kind=kreal) :: ecoord(ndof,l_max_elem_node), vect(l_max_elem_node)
542 fstrsolid%CONT_AREA(:) = 0.d0
544 if( .not.
associated(fstrsolid%CONT_SGRP_ID) )
return
546 allocate(coord(ndof*hecmesh%n_node))
547 do i=1,ndof*hecmesh%n_node
548 coord(i) = hecmesh%node(i)+fstrsolid%unode(i)
551 do i=1,ndof*hecmesh%n_node
552 coord(i) = coord(i)+fstrsolid%dunode(i)
556 do isuf=1,
size(fstrsolid%CONT_SGRP_ID)/2
557 icel = fstrsolid%CONT_SGRP_ID(2*isuf-1)
558 sid = fstrsolid%CONT_SGRP_ID(2*isuf )
560 etype = hecmesh%elem_type(icel)
561 nn = hecmw_get_max_node(etype)
562 is = hecmesh%elem_node_index(icel-1)
563 ndlocal(1:nn) = hecmesh%elem_node_item (is+1:is+nn)
566 idx = ndof*(ndlocal(i)-1)
567 ecoord(1:ndof,i) = coord(idx+1:idx+ndof)
574 fstrsolid%CONT_AREA(idx) = fstrsolid%CONT_AREA(idx) + vect(i)
583 integer(kind=kint),
intent(in) :: etype
584 integer(kind=kint),
intent(in) :: nn
585 real(kind=kreal),
intent(in) :: ecoord(:,:)
586 integer(kind=kint),
intent(in) :: sid
587 real(kind=kreal),
intent(out) :: vect(:)
589 integer(kind=kint),
parameter :: NDOF=3
590 integer(kind=kint) :: nod(l_max_surface_node)
591 integer(kind=kint) :: nsur, stype, ig0, i
592 real(kind=kreal) :: localcoord(2), normal(3), area, wg
593 real(kind=kreal) :: scoord(ndof,l_max_surface_node), h(l_max_surface_node)
597 call getsubface( etype, sid, stype, nod )
598 nsur = getnumberofnodes( stype )
600 scoord(1:ndof,i) = ecoord(1:ndof,nod(i))
604 do ig0=1,numofquadpoints( stype )
605 call getquadpoint( stype, ig0, localcoord(1:2) )
606 call getshapefunc( stype, localcoord(1:2), h(1:nsur) )
608 wg=getweight( stype, ig0 )
610 normal(1:3) = surfacenormal( stype, nsur, localcoord(1:2), scoord(1:ndof,1:nsur) )
611 area = area + dsqrt(dot_product(normal,normal))*wg
613 area = area/dble(nsur)
This module provides functions to modify MPC conditions.
subroutine fstr_delete_mpc(np, mpcs)
Delete last n equation conditions from current mpc condition.
subroutine fstr_append_mpc(np, nodes, dofs, values, mpcs)
Append new equation condition at end of existing mpc conditions.
This module defined coomon data and basic structures for analysis.
integer(kind=kint), parameter kcaslagrange
contact analysis algorithm
real(kind=kreal) dt
ANALYSIS CONTROL for NLGEOM and HEAT.
integer(kind=kint), parameter kcaalagrange
logical function fstr_iscontactactive(fstrsolid, nbc, cstep)