FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
fstr_AddContactStiff.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
12
14
15 use m_fstr
16 use elementinfo
21
22 implicit none
23
27 public :: fstr_ass_load_contact
29
30 private :: getcontactstiffness
31 private :: fstr_mat_ass_contact
32 private :: getcontactnodalforce
33 private :: gettrialfricforceandcheckfricstate
34
35contains
36
39 subroutine fstr_addcontactstiffness(cstep,iter,hecMAT,fstrMAT,fstrSOLID)
40
41 integer(kind=kint) :: cstep
42 integer(kind=kint) :: iter
43 type(hecmwst_matrix) :: hecmat
44 type(fstr_solid) :: fstrsolid
45 type(fstrst_matrix_contact_lagrange) :: fstrmat
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)
50
51 fstrmat%AL_lagrange = 0.0d0
52 fstrmat%AU_lagrange = 0.0d0
53
54 id_lagrange = 0
55
56 do i = 1, size(fstrsolid%contacts)
57
58 grpid = fstrsolid%contacts(i)%group
59 if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
60
61 do j = 1, size(fstrsolid%contacts(i)%slave)
62
63 if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
64
65 id_lagrange = id_lagrange + 1
66 lagrange = fstrmat%Lagrange(id_lagrange)
67
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)
73
74 ! Obtain contact stiffness matrix of contact pair
75 call getcontactstiffness(iter,etype,nnode,fstrsolid%contacts(i)%states(j), &
76 fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,stiffness)
77
78 ! Assemble contact stiffness matrix of contact pair into global stiffness matrix
79 call fstr_mat_ass_contact(nnode,ndlocal,id_lagrange,fstrsolid%contacts(i)%fcoeff,stiffness,hecmat,fstrmat)
80
81 enddo
82
83 enddo
84
85
86 end subroutine fstr_addcontactstiffness
87
88
90 subroutine getcontactstiffness(iter,etype,nnode,ctState,tPenalty,fcoeff,lagrange,stiffness)
91
92 type(tcontactstate) :: ctstate
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)
103
104 stiffness = 0.0d0
105
106 call getshapefunc( etype, ctstate%lpos(:), shapefunc )
107
108 normal(1:3) = ctstate%direction(1:3)
109
110 ntm(1:3) = normal(1:3)
111 do i = 1, nnode
112 ntm(i*3+1:i*3+3) = -shapefunc(i)*normal(1:3)
113 enddo
114
115 i = (nnode+1)*3 + 1
116 do j = 1, (nnode+1)*3
117 stiffness(i,j) = ntm(j); stiffness(j,i) = ntm(j)
118 enddo
119
120
121 if( fcoeff /= 0.0d0 ) then
122 if( lagrange>0.0d0 .or. iter==1 ) then
123
124 do i = 1, nnode+1
125 do j = 1, i
126 do k = 1, 3
127 do l = 1, 3
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)
129 if( k==l ) then
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)
136 endif
137 endif
138 if(i==j) cycle
139 stiffness((j-1)*3+l,(i-1)*3+k) = stiffness((i-1)*3+k,(j-1)*3+l)
140 enddo
141 enddo
142 enddo
143 enddo
144
145 if( ctstate%state == contactslip ) then
146
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
150
151 ttm(1:3) = -tangent(1:3)
152 do i = 1, nnode
153 ttm(i*3+1:i*3+3) = shapefunc(i)*tangent(1:3)
154 enddo
155
156 do i = 1, nnode+1
157 do j = 1, nnode+1
158 do k = 1, 3
159 do l = 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))
163 enddo
164 enddo
165 enddo
166 enddo
167 stiffness(1:(nnode+1)*3,1:(nnode+1)*3) = (fcoeff*lagrange/length_tft)*stiffness(1:(nnode+1)*3,1:(nnode+1)*3)
168
169 !
170 ! do i = 1, (nnode + 1)*3
171 ! do j = 1, i
172 ! do k = 1, 3
173 ! do l = 1, k
174 ! stiffness((i-1)*3+k,(j-1)*3+l) = stiffness((i-1)*3+k,(j-1)*3+l) - mut*tTm((i-1)*3+k)*tTm((j-1)*3+l))
175 ! if( k/=l )stiffness((j-1)*3+l,(i-1)*3+k) = stiffness((i-1)*3+k,(j-1)*3+l)
176 ! enddo !- l
177 ! enddo !- k
178 ! enddo !- j
179 ! enddo !- i
180 ! stiffness(1:(nnode+1)*3,1:(nnode+1)*3) = (fcoeff*dabs(lagrange)/length_tft)*stiffness(1:(nnode+1)*3,1:(nnode+1)*3)
181
182 ! j = (nnode+1)*3 + 1
183 ! do i = 1, (nnode+1)*3
184 ! stiffness(i,j) = stiffness(i,j) - fcoeff*tTm(i)
185 ! enddo
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)
187
188 endif
189 endif
190 endif
191
192 end subroutine getcontactstiffness
193
194
196 subroutine fstr_mat_ass_contact(nnode,ndLocal,id_lagrange,fcoeff,stiffness,hecMAT,fstrMAT)
197
198 type(hecmwst_matrix) :: hecmat
199 type(fstrst_matrix_contact_lagrange) :: fstrmat
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)
208
209 i = nnode + 1 + 1
210 inod = id_lagrange
211 isl = fstrmat%indexL_lagrange(inod-1)+1
212 iel = fstrmat%indexL_lagrange(inod)
213
214 do j = 1, nnode + 1
215 jnod = ndlocal(j)
216 isu = fstrmat%indexU_lagrange(jnod-1)+1
217 ieu = fstrmat%indexU_lagrange(jnod)
218
219 kl = hecmw_array_search_i(fstrmat%itemL_lagrange,isl,iel,jnod)
220 if( kl<isl .or. kl>iel ) then
221 write(*,*) '###ERROR### : cannot find connectivity (Lagrange1)'
222 stop
223 endif
224 ku = hecmw_array_search_i(fstrmat%itemU_lagrange,isu,ieu,inod)
225 if( ku<isu .or. ku>ieu ) then
226 write(*,*) '###ERROR### : cannot find connectivity (Lagrange2)'
227 stop
228 endif
229
230 idxl_base = (kl-1)*3
231 idxu_base = (ku-1)*3
232
233 do l = 1, 3
234 idxl = idxl_base + l
235 fstrmat%AL_lagrange(idxl) = fstrmat%AL_lagrange(idxl) + stiffness((i-1)*3+1,(j-1)*3+l)
236 idxu = idxu_base + l
237 fstrmat%AU_lagrange(idxu) = fstrmat%AU_lagrange(idxu) + stiffness((j-1)*3+l,(i-1)*3+1)
238 enddo
239 enddo
240
241
242 if(fcoeff /= 0.0d0)then
243
244 do i = 1, nnode + 1
245 inod = ndlocal(i)
246 do j = 1, nnode + 1
247 jnod = ndlocal(j)
248 call stf_get_block(stiffness(1:(nnode+1)*3,1:(nnode+1)*3), 3, i, j, a)
249 call hecmw_mat_add_node(hecmat, inod, jnod, a)
250 enddo
251 enddo
252
253 endif
254
255 end subroutine fstr_mat_ass_contact
256
257
260 subroutine fstr_update_ndforce_contact(cstep,hecMESH,hecMAT,fstrMAT,fstrSOLID,conMAT)
261
262 type(hecmwst_local_mesh) :: hecmesh
263 type(hecmwst_matrix) :: hecmat
264 type(fstr_solid) :: fstrsolid
265 type(fstrst_matrix_contact_lagrange) :: fstrmat
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)
273
274 integer(kind=kint) :: cstep
275 integer(kind=kint) :: grpid
276
277 id_lagrange = 0
278 if( associated(fstrsolid%CONT_NFORCE) ) fstrsolid%CONT_NFORCE(:) = 0.d0
279 if( associated(fstrsolid%CONT_FRIC) ) fstrsolid%CONT_FRIC(:) = 0.d0
280
281 do i = 1, size(fstrsolid%contacts)
282
283 grpid = fstrsolid%contacts(i)%group
284 if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
285
286 do j = 1, size(fstrsolid%contacts(i)%slave)
287
288 if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
289
290 id_lagrange = id_lagrange + 1
291 lagrange = fstrmat%Lagrange(id_lagrange)
292
293 fstrsolid%contacts(i)%states(j)%multiplier(1) = fstrmat%Lagrange(id_lagrange)
294
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)
300 do k = 1, nnode+1
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)
305 enddo
306 ! Obtain contact nodal force vector of contact pair
307 call getcontactnodalforce(etype,nnode,ndcoord,nddu,fstrsolid%contacts(i)%states(j), &
308 fstrsolid%contacts(i)%tPenalty,fstrsolid%contacts(i)%fcoeff,lagrange,ctnforce,cttforce)
309 ! Update non-eqilibrited force vector
310 if(present(conmat)) then
311 call update_ndforce_contact(nnode,ndlocal,id_lagrange,lagrange,ctnforce,cttforce,fstrsolid,conmat)
312 else
313 call update_ndforce_contact(nnode,ndlocal,id_lagrange,lagrange,ctnforce,cttforce,fstrsolid,hecmat)
314 endif
315
316 enddo
317
318 enddo
319
320 ! Consider SPC condition
321 call fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, hecmat%B)
322 if(present(conmat)) call fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, conmat%B)
323
324 end subroutine fstr_update_ndforce_contact
325
327 subroutine getcontactnodalforce(etype,nnode,ndCoord,ndDu,ctState,tPenalty,fcoeff,lagrange,ctNForce,ctTForce)
328
329 type(tcontactstate) :: ctstate
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)
340
341 ctnforce = 0.0d0
342 cttforce = 0.0d0
343
344 call getshapefunc( etype, ctstate%lpos(:), shapefunc )
345
346 normal(1:3) = ctstate%direction(1:3)
347
348 ntm(1:3) = -normal(1:3)
349 do i = 1, nnode
350 ntm(i*3+1:i*3+3) = shapefunc(i)*normal(1:3)
351 enddo
352
353 do j = 1, (nnode+1)*3
354 ctnforce(j) = lagrange*ntm(j)
355 enddo
356 j = (nnode+1)*3 + 1
357 ctnforce(j) = dot_product(ntm,ndcoord)
358
359 if(fcoeff /= 0.0d0 .and. lagrange > 0.0d0)then
360
361 call gettrialfricforceandcheckfricstate(nnode,tpenalty,fcoeff,lagrange,normal,shapefunc,ntm,nddu,ctstate)
362
363 if( ctstate%state == contactstick ) then
364 tf_final(1:3) = ctstate%tangentForce_trial(1:3)
365 elseif( ctstate%state == contactslip ) then
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)
370 endif
371
372 cttforce(1:3) = -tf_final(1:3)
373 do j = 1, nnode
374 cttforce(j*3+1:j*3+3) = shapefunc(j)*tf_final(1:3)
375 enddo
376
377 ctstate%tangentForce_final(1:3) = tf_final(1:3)
378
379 endif
380
381 end subroutine getcontactnodalforce
382
383
385 subroutine gettrialfricforceandcheckfricstate(nnode,tPenalty,fcoeff,lagrange,normal,shapefunc,nTm,ndDu,ctstate)
386
387 type(tcontactstate) :: 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
398
399 relativedisp = 0.0d0
400
401 dotp = dot_product(ntm,nddu)
402 do i = 1, 3
403 relativedisp(i) = - nddu(i)
404 do j = 1, nnode
405 relativedisp(i) = relativedisp(i) + shapefunc(j)*nddu(j*3+i)
406 enddo
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)
410 enddo
411
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
415 ctstate%state = contactstick
416 else
417 ctstate%state = contactslip
418 endif
419
420 end subroutine gettrialfricforceandcheckfricstate
421
422
425 subroutine update_ndforce_contact(nnode,ndLocal,id_lagrange,lagrange,ctNForce,ctTForce,fstrSOLID,hecMAT)
426
427 type(fstr_solid) :: fstrsolid
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)
437
438 np = hecmat%NP; ndof = hecmat%NDOF
439
440 do i = 1, nnode + 1
441 inod = ndlocal(i)
442 idx = (inod-1)*3+1
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)
447 enddo
448
449 hecmat%B(np*ndof+id_lagrange) = ctnforce((nnode+1)*3+1)+cttforce((nnode+1)*3+1)
450
451 end subroutine update_ndforce_contact
452
453
456 subroutine fstr_ass_load_contact(cstep, hecMESH, hecMAT, fstrSOLID, fstrMAT)
457
458 type(hecmwst_local_mesh) :: hecmesh
459 type(hecmwst_matrix) :: hecmat
460 type(fstr_solid) :: fstrsolid
461 type(fstrst_matrix_contact_lagrange) :: fstrmat
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)
471
472 np = hecmat%NP ; ndof = hecmat%NDOF
473
474 id_lagrange = 0
475
476 do i = 1, size(fstrsolid%contacts)
477
478 grpid = fstrsolid%contacts(i)%group
479 if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
480
481 do j = 1, size(fstrsolid%contacts(i)%slave)
482
483 if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
484
485 id_lagrange = id_lagrange + 1
486 lagrange = fstrsolid%contacts(i)%states(j)%multiplier(1)
487
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)
493 do k = 1, nnode+1
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)
496 enddo
497
498 ctforce = 0.0d0
499
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)
503 do k = 1, nnode
504 ntm(k*3+1:k*3+3) = shapefunc(k)*normal(1:3)
505 enddo
506 do l = 1, (nnode+1)*3
507 ctforce(l) = lagrange*ntm(l)
508 enddo
509 l = (nnode+1)*3 + 1
510 ctforce(l) = dot_product(ntm(1:(nnode+1)*3),ndcoord(1:(nnode+1)*3))
511
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)
515 do l = 1, nnode
516 ctforce(l*3+1:l*3+3) = ctforce(l*3+1:l*3+3) + shapefunc(l)*tf_final(1:3)
517 enddo
518 endif
519
520 do l = 1, nnode + 1
521 lnod = ndlocal(l)
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)
523 enddo
524 hecmat%B(np*ndof+id_lagrange) = ctforce((nnode+1)*3+1)
525
526 enddo
527
528 enddo
529
530 end subroutine fstr_ass_load_contact
531
532
535 subroutine fstr_mat_ass_bc_contact(hecMAT,fstrMAT,inode,idof,RHS)
536
537 type(hecmwst_matrix) :: hecmat
538 type(fstrst_matrix_contact_lagrange) :: fstrmat
539 integer(kind=kint) :: inode, idof
540 integer(kind=kint) :: isu, ieu, isl, iel, i, l, k
541 real(kind=kreal) :: rhs
542
543 isu = fstrmat%indexU_lagrange(inode-1)+1
544 ieu = fstrmat%indexU_lagrange(inode)
545 do i = isu, ieu
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)
550 k = hecmw_array_search_i(fstrmat%itemL_lagrange,isl,iel,inode)
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
554 enddo
555
556 end subroutine fstr_mat_ass_bc_contact
557
558end module m_addcontactstiffness
559
560
561
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
This module provides functions of reconstructing.
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 functions: 1) obtain contact stiffness matrix of each contact pair and assemble ...
subroutine, public fstr_ass_load_contact(cstep, hecmesh, hecmat, fstrsolid, fstrmat)
This subroutine adds initial contact force vector to the right-hand side vector \at the beginning of ...
subroutine, public fstr_mat_ass_bc_contact(hecmat, fstrmat, inode, idof, rhs)
Modify Lagrange multiplier-related part of stiffness matrix and right-hand side vector for dealing wi...
subroutine, public fstr_update_ndforce_contact(cstep, hecmesh, hecmat, fstrmat, fstrsolid, conmat)
This subroutine obtains contact nodal force vector of each contact pair and assembles it into right-h...
subroutine, public update_ndforce_contact(nnode, ndlocal, id_lagrange, lagrange, ctnforce, cttforce, fstrsolid, hecmat)
This subroutine assembles contact nodal force vector into right-hand side vector to update non-equili...
subroutine, public fstr_addcontactstiffness(cstep, iter, hecmat, fstrmat, fstrsolid)
This subroutine obtains contact stiffness matrix of each contact pair and assembles it into global st...
This module provide functions of contact stiffness calculation.
Definition: contact_lib.f90:6
integer, parameter contactslip
Definition: contact_lib.f90:18
integer, parameter contactfree
contact state definition
Definition: contact_lib.f90:16
integer, parameter contactstick
Definition: contact_lib.f90:17
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.
Definition: m_fstr.f90:15
logical function fstr_iscontactactive(fstrsolid, nbc, cstep)
Definition: m_fstr.f90:1021
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...
This structure records contact status.
Definition: contact_lib.f90:27