30 private :: getcontactstiffness
31 private :: fstr_mat_ass_contact
32 private :: getcontactnodalforce
33 private :: gettrialfricforceandcheckfricstate
41 integer(kind=kint) :: cstep
42 integer(kind=kint) :: iter
43 type(hecmwst_matrix) :: hecmat
46 integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(9)
47 integer(kind=kint) :: i, j, id_lagrange, grpid
48 real(kind=kreal) :: lagrange
49 real(kind=kreal) :: stiffness(9*3 + 1, 9*3 + 1)
51 fstrmat%AL_lagrange = 0.0d0
52 fstrmat%AU_lagrange = 0.0d0
56 do i = 1,
size(fstrsolid%contacts)
58 grpid = fstrsolid%contacts(i)%group
61 do j = 1,
size(fstrsolid%contacts(i)%slave)
63 if( fstrsolid%contacts(i)%states(j)%state ==
contactfree ) cycle
65 id_lagrange = id_lagrange + 1
66 lagrange = fstrmat%Lagrange(id_lagrange)
68 ctsurf = fstrsolid%contacts(i)%states(j)%surface
69 etype = fstrsolid%contacts(i)%master(ctsurf)%etype
70 nnode =
size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
71 ndlocal(1) = fstrsolid%contacts(i)%slave(j)
72 ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
75 call getcontactstiffness(iter,etype,nnode,fstrsolid%contacts(i)%states(j), &
76 fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,stiffness)
79 call fstr_mat_ass_contact(nnode,ndlocal,id_lagrange,fstrsolid%contacts(i)%fcoeff,stiffness,hecmat,fstrmat)
90 subroutine getcontactstiffness(iter,etype,nnode,ctState,tPenalty,fcoeff,lagrange,stiffness)
93 integer(kind=kint) :: iter
94 integer(kind=kint) :: etype, nnode
95 integer(kind=kint) :: i, j, k, l
96 real(kind=kreal) :: normal(3), shapefunc(nnode)
97 real(kind=kreal) :: ntm((nnode + 1)*3)
98 real(kind=kreal) :: fcoeff, tpenalty
99 real(kind=kreal) :: lagrange
100 real(kind=kreal) :: tf_trial(3), length_tft
101 real(kind=kreal) :: tangent(3), ttm((nnode + 1)*3)
102 real(kind=kreal) :: stiffness(9*3 + 1, 9*3 + 1)
108 normal(1:3) = ctstate%direction(1:3)
110 ntm(1:3) = normal(1:3)
112 ntm(i*3+1:i*3+3) = -shapefunc(i)*normal(1:3)
116 do j = 1, (nnode+1)*3
117 stiffness(i,j) = ntm(j); stiffness(j,i) = ntm(j)
121 if( fcoeff /= 0.0d0 )
then
122 if( lagrange>0.0d0 .or. iter==1 )
then
128 stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) - tpenalty*ntm((i-1)*3+k)*ntm((j-1)*3+l)
130 if(i==1 .and. j==1)
then
131 stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) + tpenalty
132 elseif(i>1 .and. j==1)
then
133 stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) - tpenalty*shapefunc(i-1)
134 elseif(i>1 .and. j>1)
then
135 stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) + tpenalty*shapefunc(i-1)*shapefunc(j-1)
139 stiffness((j-1)*3+l,(i-1)*3+k) = stiffness((i-1)*3+k,(j-1)*3+l)
147 tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
148 length_tft = dsqrt(dot_product(tf_trial,tf_trial))
149 tangent(1:3) = tf_trial(1:3)/length_tft
151 ttm(1:3) = -tangent(1:3)
153 ttm(i*3+1:i*3+3) = shapefunc(i)*tangent(1:3)
160 stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) &
161 + tpenalty*(-ttm((i-1)*3+k)*ttm((j-1)*3+l) &
162 +ttm((i-1)*3+k)*ntm((j-1)*3+l)*dot_product(tangent,normal))
167 stiffness(1:(nnode+1)*3,1:(nnode+1)*3) = (fcoeff*lagrange/length_tft)*stiffness(1:(nnode+1)*3,1:(nnode+1)*3)
186 stiffness(1:(nnode+1)*3,(nnode+1)*3+1) = stiffness(1:(nnode+1)*3,(nnode+1)*3+1) - fcoeff*ttm(1:(nnode+1)*3)
192 end subroutine getcontactstiffness
196 subroutine fstr_mat_ass_contact(nnode,ndLocal,id_lagrange,fcoeff,stiffness,hecMAT,fstrMAT)
198 type(hecmwst_matrix) :: hecmat
200 integer(kind=kint) :: nnode, ndlocal(nnode + 1), id_lagrange
203 integer(kind=kint) :: i, j, inod, jnod, l
204 integer(kind=kint) :: isl, iel, idxl_base, kl, idxl, isu, ieu, idxu_base, ku, idxu
205 real(kind=kreal) :: fcoeff
206 real(kind=kreal) :: stiffness(9*3 + 1, 9*3 + 1)
207 real(kind=kreal) :: a(3, 3)
211 isl = fstrmat%indexL_lagrange(inod-1)+1
212 iel = fstrmat%indexL_lagrange(inod)
216 isu = fstrmat%indexU_lagrange(jnod-1)+1
217 ieu = fstrmat%indexU_lagrange(jnod)
220 if( kl<isl .or. kl>iel )
then
221 write(*,*)
'###ERROR### : cannot find connectivity (Lagrange1)'
225 if( ku<isu .or. ku>ieu )
then
226 write(*,*)
'###ERROR### : cannot find connectivity (Lagrange2)'
235 fstrmat%AL_lagrange(idxl) = fstrmat%AL_lagrange(idxl) + stiffness((i-1)*3+1,(j-1)*3+l)
237 fstrmat%AU_lagrange(idxu) = fstrmat%AU_lagrange(idxu) + stiffness((j-1)*3+l,(i-1)*3+1)
242 if(fcoeff /= 0.0d0)
then
248 call stf_get_block(stiffness(1:(nnode+1)*3,1:(nnode+1)*3), 3, i, j, a)
255 end subroutine fstr_mat_ass_contact
262 type(hecmwst_local_mesh) :: hecmesh
263 type(hecmwst_matrix) :: hecmat
266 type(hecmwst_matrix),
optional :: conmat
267 integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(9)
268 integer(kind=kint) :: i, j, k, id_lagrange
269 real(kind=kreal) :: ndcoord(9*3), nddu(9*3)
270 real(kind=kreal) :: lagrange
271 real(kind=kreal) :: ctnforce(9*3+1)
272 real(kind=kreal) :: cttforce(9*3+1)
274 integer(kind=kint) :: cstep
275 integer(kind=kint) :: grpid
278 if(
associated(fstrsolid%CONT_NFORCE) ) fstrsolid%CONT_NFORCE(:) = 0.d0
279 if(
associated(fstrsolid%CONT_FRIC) ) fstrsolid%CONT_FRIC(:) = 0.d0
281 do i = 1,
size(fstrsolid%contacts)
283 grpid = fstrsolid%contacts(i)%group
286 do j = 1,
size(fstrsolid%contacts(i)%slave)
288 if( fstrsolid%contacts(i)%states(j)%state ==
contactfree ) cycle
290 id_lagrange = id_lagrange + 1
291 lagrange = fstrmat%Lagrange(id_lagrange)
293 fstrsolid%contacts(i)%states(j)%multiplier(1) = fstrmat%Lagrange(id_lagrange)
295 ctsurf = fstrsolid%contacts(i)%states(j)%surface
296 etype = fstrsolid%contacts(i)%master(ctsurf)%etype
297 nnode =
size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
298 ndlocal(1) = fstrsolid%contacts(i)%slave(j)
299 ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
301 ndcoord((k-1)*3+1:(k-1)*3+3) = hecmesh%node((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) &
302 + fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) &
303 + fstrsolid%dunode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
304 nddu((k-1)*3+1:(k-1)*3+3) = fstrsolid%dunode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
307 call getcontactnodalforce(etype,nnode,ndcoord,nddu,fstrsolid%contacts(i)%states(j), &
308 fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,ctnforce,cttforce)
310 if(
present(conmat))
then
327 subroutine getcontactnodalforce(etype,nnode,ndCoord,ndDu,ctState,tPenalty,fcoeff,lagrange,ctNForce,ctTForce)
330 integer(kind=kint) :: etype, nnode
331 integer(kind=kint) :: i, j
332 real(kind=kreal) :: fcoeff, tpenalty
333 real(kind=kreal) :: lagrange
334 real(kind=kreal) :: ndcoord((nnode + 1)*3), nddu((nnode + 1)*3)
335 real(kind=kreal) :: normal(3), shapefunc(nnode)
336 real(kind=kreal) :: ntm((nnode + 1)*3)
337 real(kind=kreal) :: tf_trial(3), length_tft, tangent(3), tf_final(3)
338 real(kind=kreal) :: ctnforce((nnode+1)*3+1)
339 real(kind=kreal) :: cttforce((nnode+1)*3+1)
346 normal(1:3) = ctstate%direction(1:3)
348 ntm(1:3) = -normal(1:3)
350 ntm(i*3+1:i*3+3) = shapefunc(i)*normal(1:3)
353 do j = 1, (nnode+1)*3
354 ctnforce(j) = lagrange*ntm(j)
357 ctnforce(j) = dot_product(ntm,ndcoord)
359 if(fcoeff /= 0.0d0 .and. lagrange > 0.0d0)
then
361 call gettrialfricforceandcheckfricstate(nnode,tpenalty,fcoeff,lagrange,normal,shapefunc,ntm,nddu,ctstate)
364 tf_final(1:3) = ctstate%tangentForce_trial(1:3)
366 tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
367 length_tft = dsqrt(dot_product(tf_trial,tf_trial))
368 tangent(1:3) = tf_trial(1:3)/length_tft
369 tf_final(1:3) = fcoeff*dabs(lagrange)*tangent(1:3)
372 cttforce(1:3) = -tf_final(1:3)
374 cttforce(j*3+1:j*3+3) = shapefunc(j)*tf_final(1:3)
377 ctstate%tangentForce_final(1:3) = tf_final(1:3)
381 end subroutine getcontactnodalforce
385 subroutine gettrialfricforceandcheckfricstate(nnode,tPenalty,fcoeff,lagrange,normal,shapefunc,nTm,ndDu,ctstate)
388 integer(kind=kint) :: nnode
389 integer(kind=kint) :: i, j
390 real(kind=kreal) :: fcoeff, tpenalty
391 real(kind=kreal) :: lagrange
392 real(kind=kreal) :: nddu((nnode + 1)*3)
393 real(kind=kreal) :: normal(3), shapefunc(nnode)
394 real(kind=kreal) :: ntm((nnode + 1)*3)
395 real(kind=kreal) :: dotp
396 real(kind=kreal) :: relativedisp(3)
397 real(kind=kreal) :: tf_yield
401 dotp = dot_product(ntm,nddu)
403 relativedisp(i) = - nddu(i)
405 relativedisp(i) = relativedisp(i) + shapefunc(j)*nddu(j*3+i)
407 relativedisp(i) = relativedisp(i) - dotp*normal(i)
408 ctstate%reldisp(i) = -relativedisp(i)
409 ctstate%tangentForce_trial(i) = ctstate%tangentForce1(i) -tpenalty*relativedisp(i)
412 tf_yield = fcoeff*dabs(lagrange)
413 if(ctstate%state ==
contactslip) tf_yield =0.99d0*tf_yield
414 if( dsqrt(dot_product(ctstate%tangentForce_trial,ctstate%tangentForce_trial)) <= tf_yield )
then
420 end subroutine gettrialfricforceandcheckfricstate
428 type(hecmwst_matrix) :: hecmat
429 integer(kind=kint) :: nnode, ndlocal(nnode + 1)
431 integer(kind=kint) :: id_lagrange
432 real(kind=kreal) :: lagrange
433 integer(kind=kint) :: np, ndof
434 integer (kind=kint) :: i, inod, idx
435 real(kind=kreal) :: ctnforce((nnode+1)*3+1)
436 real(kind=kreal) :: cttforce((nnode+1)*3+1)
438 np = hecmat%NP; ndof = hecmat%NDOF
443 hecmat%B(idx:idx+2) = hecmat%B(idx:idx+2) + ctnforce((i-1)*3+1:(i-1)*3+3) + cttforce((i-1)*3+1:(i-1)*3+3)
444 if( lagrange < 0.d0 ) cycle
445 fstrsolid%CONT_NFORCE(idx:idx+2) = fstrsolid%CONT_NFORCE(idx:idx+2) + ctnforce((i-1)*3+1:(i-1)*3+3)
446 fstrsolid%CONT_FRIC(idx:idx+2) = fstrsolid%CONT_FRIC(idx:idx+2) + cttforce((i-1)*3+1:(i-1)*3+3)
449 hecmat%B(np*ndof+id_lagrange) = ctnforce((nnode+1)*3+1)+cttforce((nnode+1)*3+1)
458 type(hecmwst_local_mesh) :: hecmesh
459 type(hecmwst_matrix) :: hecmat
462 integer(kind=kint) :: cstep
463 integer(kind=kint) :: np, ndof
464 integer(kind=kint) :: i, j, k, l, id_lagrange, lnod, grpid
465 integer(kind=kint) :: ctsurf, etype, nnode, ndlocal(9)
466 real(kind=kreal) :: ndcoord(9*3), lagrange
467 real(kind=kreal) :: normal(3), shapefunc(9)
468 real(kind=kreal) :: ntm(10*3)
469 real(kind=kreal) :: tf_final(3)
470 real(kind=kreal) :: ctforce(9*3 + 1)
472 np = hecmat%NP ; ndof = hecmat%NDOF
476 do i = 1,
size(fstrsolid%contacts)
478 grpid = fstrsolid%contacts(i)%group
481 do j = 1,
size(fstrsolid%contacts(i)%slave)
483 if( fstrsolid%contacts(i)%states(j)%state ==
contactfree ) cycle
485 id_lagrange = id_lagrange + 1
486 lagrange = fstrsolid%contacts(i)%states(j)%multiplier(1)
488 ctsurf = fstrsolid%contacts(i)%states(j)%surface
489 etype = fstrsolid%contacts(i)%master(ctsurf)%etype
490 nnode =
size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
491 ndlocal(1) = fstrsolid%contacts(i)%slave(j)
492 ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
494 ndcoord((k-1)*3+1:(k-1)*3+3) = hecmesh%node((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3) &
495 + fstrsolid%unode((ndlocal(k)-1)*3+1:(ndlocal(k)-1)*3+3)
500 call getshapefunc( etype, fstrsolid%contacts(i)%states(j)%lpos(:), shapefunc )
501 normal(1:3) = fstrsolid%contacts(i)%states(j)%direction(1:3)
502 ntm(1:3) = -normal(1:3)
504 ntm(k*3+1:k*3+3) = shapefunc(k)*normal(1:3)
506 do l = 1, (nnode+1)*3
507 ctforce(l) = lagrange*ntm(l)
510 ctforce(l) = dot_product(ntm(1:(nnode+1)*3),ndcoord(1:(nnode+1)*3))
512 if( fstrsolid%contacts(i)%fcoeff/=0.0d0 .and. lagrange>0.0d0 )
then
513 tf_final(1:3) = fstrsolid%contacts(i)%states(j)%tangentForce_final(1:3)
514 ctforce(1:3) = ctforce(1:3) - tf_final(1:3)
516 ctforce(l*3+1:l*3+3) = ctforce(l*3+1:l*3+3) + shapefunc(l)*tf_final(1:3)
522 hecmat%B((lnod-1)*3+1:(lnod-1)*3+3) = hecmat%B((lnod-1)*3+1:(lnod-1)*3+3) + ctforce((l-1)*3+1:(l-1)*3+3)
524 hecmat%B(np*ndof+id_lagrange) = ctforce((nnode+1)*3+1)
537 type(hecmwst_matrix) :: hecmat
539 integer(kind=kint) :: inode, idof
540 integer(kind=kint) :: isu, ieu, isl, iel, i, l, k
541 real(kind=kreal) :: rhs
543 isu = fstrmat%indexU_lagrange(inode-1)+1
544 ieu = fstrmat%indexU_lagrange(inode)
546 fstrmat%AU_lagrange((i-1)*3+idof) = 0.0d0
547 l = fstrmat%itemU_lagrange(i)
548 isl = fstrmat%indexL_lagrange(l-1)+1
549 iel = fstrmat%indexL_lagrange(l)
551 if(k < isl .or. k > iel) cycle
552 hecmat%B(hecmat%NP*hecmat%NDOF+l) = hecmat%B(hecmat%NP*hecmat%NDOF+l) - fstrmat%AL_lagrange((k-1)*3+idof)*rhs
553 fstrmat%AL_lagrange((k-1)*3+idof) = 0.0d0
This module encapsulate the basic functions of all elements provide by this software.
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
subroutine, public hecmw_mat_add_node(hecmat, inod, jnod, a)
integer(kind=kint) function, public hecmw_array_search_i(array, is, ie, ival)
subroutine, public stf_get_block(stiffness, ndof, inod, jnod, a)
This module provides function to calcualte residual of nodal force.
subroutine, public fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, b)
This module defined coomon data and basic structures for analysis.
logical function fstr_iscontactactive(fstrsolid, nbc, cstep)