FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
fstr_contact.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!-------------------------------------------------------------------------------
6module mcontact
7
9 use hecmw
10 use m_fstr
11 implicit none
12
13 private :: l_contact2mpc, l_tied2mpc
14 integer(kind=kint), save :: n_contact_mpc
15 logical, private :: active
16
17 real(kind=kreal), save :: mu=1.d10
18 real(kind=kreal), save :: mut=1.d6
19 real(kind=kreal), save :: cdotp=1.d3
20
21 real(kind=kreal), save :: cgn=1.d-5
22 real(kind=kreal), save :: cgt=1.d-3
23
24 real(kind=kreal), save :: gnt(2)
26 real(kind=kreal), save :: bakgnt(2)
28
29contains
30
32 subroutine print_contatct_pair( file, pair )
33 integer(kind=kint), intent(in) :: file
34 type( hecmwst_contact_pair ), intent(in) :: pair
35
36 integer(kind=kint) :: i
37 write(file,*) "Number of contact pair", pair%n_pair
38 do i=1,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)
41 enddo
42 end subroutine
43
44 subroutine fstr_set_contact_penalty( maxv )
45 real(kind=kreal), intent(in) :: maxv
46 mu = cdotp * maxv
47 if( gnt(1)<1.d-3 ) mu=cdotp*10.d0* maxv
48 bakgnt = 0.d0
49 end subroutine
50
51 logical function fstr_is_contact_active()
53 end function
54
56 logical, intent(in) :: a
57 active = a
58 end subroutine
59
60 logical function fstr_is_contact_conv(ctAlgo,infoCTChange,hecMESH)
61 integer(kind=kint), intent(in) :: ctalgo
62 type (fstr_info_contactchange), intent(in) :: infoctchange
63 type (hecmwst_local_mesh), intent(in) :: hecmesh
64 integer(kind=kint) :: is_conv
65 fstr_is_contact_conv = .false.
66 if( ctalgo== kcaslagrange ) then
67 if( infoctchange%contact2free+infoctchange%contact2neighbor+ &
68 infoctchange%contact2difflpos+infoctchange%free2contact == 0 ) &
70 elseif( ctalgo==kcaalagrange ) then
71 if( gnt(1)<cgn .and. gnt(2)<cgt ) fstr_is_contact_conv = .true.
72 write(*,'(a,2e15.7)') "--Relative displacement in contact surface",gnt
73 ! if( dabs( bakgnt(1)-gnt(1) )<cgn*1.d-2 .and. &
74 ! dabs( bakgnt(2)-gnt(2) )<cgn ) fstr_is_contact_conv = .true.
75 bakgnt = gnt
76 endif
77 is_conv = 0
78 if (fstr_is_contact_conv) is_conv = 1
79 call hecmw_allreduce_i1(hecmesh, is_conv, hecmw_min)
80 if (is_conv == 0) then
81 fstr_is_contact_conv = .false.
82 else
84 endif
85 end function
86
87 logical function fstr_is_matrixstructure_changed(infoCTChange)
88 type (fstr_info_contactchange) :: infoctchange
90 if( infoctchange%contact2free+infoctchange%contact2neighbor+infoctchange%free2contact > 0 ) &
92 end function
93
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 ! 3D problem only, currently
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)
105 nmpc=0
106 do i=1,size(contact%states)
107 if( contact%states(i)%state == -1 ) cycle ! in free
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) )
114 tdof = 0
115 do j=1,ndof
116 if( dabs(values(j))<tol ) cycle
117 tdof = tdof+1
118 nodes(tdof) = contact%slave(i)
119 dofs(tdof) = j
120 val(tdof) = values(j)
121 enddo
122 do j=1,nenode
123 nn = contact%master(csurf)%nodes(j)
124 nodes( j*ndof+1:j*ndof+ndof ) = nn
125 do k=1,ndof
126 if( dabs(values(j*ndof+k)) < tol ) cycle
127 tdof=tdof+1
128 nodes(tdof)=nn
129 dofs(tdof ) = k
130 val(tdof)=values(j*ndof+k)
131 enddo
132 enddo
133 val(tdof+1) = values(nenode*ndof+ndof+1)
134
135 call fstr_append_mpc( tdof, nodes(1:tdof), dofs(1:tdof), val(1:tdof+1), mpcs )
136 nmpc=nmpc+1
137 enddo
138 end subroutine l_contact2mpc
139
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)
149 nmpc=0
150 do i=1,size(contact%slave)
151 csurf = contact%states(i)%surface
152 if( csurf<=0 ) cycle ! contactor not exists
153 nenode = size(contact%master(csurf)%nodes)
154 tdof = nenode+1
155 nodes(1) = contact%slave(i)
156 nodes( 2:tdof ) = contact%master(csurf)%nodes(:)
157 values(1) = -1.d0
158 values(2:tdof) = 1.d0
159 values(tdof+1) = 0.d0
160 etype = contact%master(csurf)%etype
161 do j=1,3
162 dofs(1:tdof) = j
163 call fstr_append_mpc( tdof, nodes(1:tdof), dofs(1:tdof), values(1:tdof+1), mpcs )
164 nmpc=nmpc+1
165 enddo
166 enddo
167 end subroutine l_tied2mpc
168
170 subroutine fstr_contact2mpc( contacts, mpcs )
171 type( tcontact ), intent(in) :: contacts(:)
172 type( hecmwst_mpc ), intent(inout) :: mpcs
173 integer(kind=kint) :: i, nmpc
174 n_contact_mpc = 0
175 do i=1,size(contacts)
176 if( contacts(i)%algtype == contactunknown ) cycle ! not initialized
177 if( contacts(i)%algtype == contactfslid ) then
178 print *, "Cannot deal with finit slip problems by MPC!"
179 cycle
180 endif
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 )
187 endif
188 enddo
189 end subroutine
190
192 subroutine fstr_del_contactmpc( mpcs )
194 type( hecmwst_mpc ), intent(inout) :: mpcs
195 call fstr_delete_mpc( n_contact_mpc, mpcs )
196 end subroutine
197
199 subroutine fstr_write_mpc( file, mpcs )
200 integer(kind=kint), intent(in) :: file
201 type( hecmwst_mpc ), intent(in) :: mpcs
202
203 integer(kind=kint) :: i,j,n0,n1
204 write(file, *) "Number of equation", mpcs%n_mpc
205 do i=1,mpcs%n_mpc
206 write(file,*) "--Equation",i
207 n0=mpcs%mpc_index(i-1)+1
208 n1=mpcs%mpc_index(i)
209 write(file, *) n0,n1
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)
213 enddo
214 end subroutine
215
217 subroutine fstr_scan_contact_state( cstep, sub_step, cont_step, dt, ctAlgo, hecMESH, fstrSOLID, infoCTChange, B )
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
224 type(fstr_solid), intent(inout) :: fstrSOLID
225 type(fstr_info_contactchange), intent(inout):: infoCTChange
226 ! logical, intent(inout) :: changed !< if contact state changed
227 real(kind=kreal), optional :: b(:)
228 character(len=9) :: flag_ctAlgo
229 integer(kind=kint) :: i, grpid
230 logical :: iactive, is_init
231
232 fstrsolid%CONT_RELVEL(:) = 0.d0
233 fstrsolid%CONT_STATE(:) = 0.d0
234
235 if( ctalgo == kcaslagrange ) then
236 flag_ctalgo = 'SLagrange'
237 elseif( ctalgo == kcaalagrange ) then
238 flag_ctalgo = 'ALagrange'
239 endif
240
241 ! P.A. We redefine fstrSOLID%ddunode as current coordinate of every nodes
242 ! fstrSOLID%ddunode(:) = fstrSOLID%unode(:) + fstrSOLID%dunode(:)
243 do i = 1, size(fstrsolid%unode)
244 fstrsolid%ddunode(i) = hecmesh%node(i) + fstrsolid%unode(i) + fstrsolid%dunode(i)
245 enddo
246 active = .false.
247
248 infoctchange%contact2free = 0
249 infoctchange%contact2neighbor = 0
250 infoctchange%contact2diffLpos = 0
251 infoctchange%free2contact = 0
252 infoctchange%contactNode_current = 0
253
254 is_init = ( cstep == 1 .and. sub_step == 1 .and. cont_step == 0 )
255
256 do i=1,size(fstrsolid%contacts)
257 grpid = fstrsolid%contacts(i)%group
258 if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) then
259 call clear_contact_state(fstrsolid%contacts(i)); cycle
260 endif
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 )
264 else
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 )
267 endif
268 if( .not. active ) active = iactive
269
270 !for output contact state
271 call set_contact_state_vector( fstrsolid%contacts(i), dt, fstrsolid%CONT_RELVEL, fstrsolid%CONT_STATE )
272 enddo
273
274 infoctchange%contactNode_current = infoctchange%contactNode_previous+infoctchange%free2contact-infoctchange%contact2free
275 infoctchange%contactNode_previous = infoctchange%contactNode_current
276
277 if( .not. active ) then
278 fstrsolid%CONT_NFORCE(:) = 0.d0
279 fstrsolid%CONT_FRIC(:) = 0.d0
280 end if
281
282 end subroutine
283
285 subroutine fstr_scan_contact_state_exp( cstep, hecMESH, fstrSOLID, infoCTChange )
286 integer(kind=kint), intent(in) :: cstep
287 type( hecmwst_local_mesh ), intent(in) :: hecMESH
288 type(fstr_solid), intent(inout) :: fstrSOLID
289 type(fstr_info_contactchange), intent(inout) :: infoCTChange
290
291 integer(kind=kint) :: i
292 logical :: iactive, is_init
293
294
295 ! P.A. We redefine fstrSOLID%ddunode as current coordinate of every nodes
296 ! fstrSOLID%ddunode(:) = fstrSOLID%unode(:) + fstrSOLID%dunode(:)
297 do i = 1, size(fstrsolid%unode)
298 fstrsolid%ddunode(i) = hecmesh%node(i) + fstrsolid%unode(i) + fstrsolid%dunode(i)
299 enddo
300 infoctchange%active = .false.
301
302 infoctchange%contact2free = 0
303 infoctchange%contact2neighbor = 0
304 infoctchange%contact2diffLpos = 0
305 infoctchange%free2contact = 0
306 infoctchange%contactNode_current = 0
307
308 is_init = ( cstep == 1 )
309
310 do i=1,size(fstrsolid%contacts)
311 ! grpid = fstrSOLID%contacts(i)%group
312 ! if( .not. fstr_isContactActive( fstrSOLID, grpid, cstep ) ) then
313 ! call clear_contact_state(fstrSOLID%contacts(i)); cycle
314 ! endif
315
316 call scan_contact_state_exp( fstrsolid%contacts(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
317 & infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive )
318
319 infoctchange%active = infoctchange%active .or. iactive
320 enddo
321
322 infoctchange%contactNode_current = infoctchange%contactNode_previous+infoctchange%free2contact-infoctchange%contact2free
323 infoctchange%contactNode_previous = infoctchange%contactNode_current
324 fstrsolid%ddunode = 0.d0
325 end subroutine
326
328 subroutine fstr_update_contact0( hecMESH, fstrSOLID, B )
329 type( hecmwst_local_mesh ), intent(in) :: hecMESH
330 type(fstr_solid), intent(inout) :: fstrSOLID
331 real(kind=kreal), intent(inout) :: b(:)
332
333 integer(kind=kint) :: i
334 do i=1, size(fstrsolid%contacts)
335 ! if( contacts(i)%mpced ) cycle
336 call calcu_contact_force0( fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:) &
337 , fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff, mu, mut, b )
338 enddo
339 end subroutine
340
342 subroutine fstr_update_contact_multiplier( hecMESH, fstrSOLID, ctchanged )
343 type( hecmwst_local_mesh ), intent(in) :: hecMESH
344 type(fstr_solid), intent(inout) :: fstrSOLID
345 logical, intent(out) :: ctchanged
346
347 integer(kind=kint) :: i, nc
348 gnt = 0.d0; ctchanged = .false.
349 nc = size(fstrsolid%contacts)
350 do i=1, nc
351 ! if( contacts(i)%mpced ) cycle
352 call update_contact_multiplier( fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:) &
353 , fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff, mu, mut, gnt, ctchanged )
354 enddo
355 if( nc>0 ) gnt = gnt/nc
356 end subroutine
357
359 subroutine fstr_update_contact_tangentforce( fstrSOLID )
360 type(fstr_solid), intent(inout) :: fstrSOLID
361
362 integer(kind=kint) :: i, nc
363
364 nc = size(fstrsolid%contacts)
365 do i=1, nc
366 call update_contact_tangentforce( fstrsolid%contacts(i) )
367 enddo
368 end subroutine
369
371 subroutine fstr_contactbc( iter, hecMESH, hecMAT, fstrSOLID )
373 integer(kind=kint), intent(in) :: iter
374 type (hecmwST_local_mesh), intent(inout) :: hecMESH
375 type (hecmwST_matrix), intent(inout) :: hecMAT
376 type(fstr_solid), intent(inout) :: fstrSOLID
377
378 integer(kind=kint), parameter :: NDOF=3
379
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)
385 ! if( n_contact_mpc>0 ) call fstr_delete_mpc( n_contact_mpc, hecMESH%mpc )
386 ! call fstr_contact2mpc( fstrSOLID%contacts, hecMESH%mpc )
387 ! temp. need modification
388 ! do i=1,size(fstrSOLID%contacts)
389 ! fstrSOLID%contacts(i)%mpced = .true.
390 ! enddo
391 ! call fstr_write_mpc( 6, hecMESH%mpc )
392 !print *,"Contact to mpc ok!",n_contact_mpc,"equations generated"
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 ! free
398 ctsurf = fstrsolid%contacts(i)%states(j)%surface ! contacting 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)
402 do k=1,nnode
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))
407 enddo
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 )
412
413 if( iter>1 ) cycle
414 ! if( fstrSOLID%contacts(i)%states(j)%multiplier(1)/=0.d0 ) cycle
415 ! In case of new contact nodes, add enforced disp constraint
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
419 do m=1,nnode+1
420 nd = ndlocal(m)
421 do k=1,ndof
422 hecmat%B(ndof*(nd-1)+k)=hecmat%B(ndof*(nd-1)+k)-force((m-1)*ndof+k)
423 enddo
424 enddo
425
426 enddo
427 enddo
428
429 end subroutine
430
431 subroutine initialize_contact_output_vectors(fstrSOLID,hecMAT)
432 type(fstr_solid) :: fstrsolid
433 type(hecmwst_matrix) :: hecMAT
434
435 if( .not. associated(fstrsolid%CONT_NFORCE) ) then
436 allocate( fstrsolid%CONT_NFORCE(hecmat%NP*3) )
437 fstrsolid%CONT_NFORCE(:) = 0.d0
438 end if
439
440 if( .not. associated(fstrsolid%CONT_FRIC) ) then
441 allocate( fstrsolid%CONT_FRIC(hecmat%NP*3) )
442 fstrsolid%CONT_FRIC(:) = 0.d0
443 end if
444
445 if( .not. associated(fstrsolid%CONT_RELVEL) ) then
446 allocate( fstrsolid%CONT_RELVEL(hecmat%NP*3) )
447 fstrsolid%CONT_RELVEL(:) = 0.d0
448 end if
449
450 if( .not. associated(fstrsolid%CONT_STATE) ) then
451 allocate( fstrsolid%CONT_STATE(hecmat%NP*1) )
452 fstrsolid%CONT_STATE(:) = 0.d0
453 end if
454
455 if( .not. associated(fstrsolid%CONT_AREA) ) then
456 allocate( fstrsolid%CONT_AREA(hecmat%NP) )
457 fstrsolid%CONT_AREA(:) = 0.d0
458 end if
459
460 if( .not. associated(fstrsolid%CONT_NTRAC) ) then
461 allocate( fstrsolid%CONT_NTRAC(hecmat%NP*3) )
462 fstrsolid%CONT_NTRAC(:) = 0.d0
463 end if
464
465 if( .not. associated(fstrsolid%CONT_FTRAC) ) then
466 allocate( fstrsolid%CONT_FTRAC(hecmat%NP*3) )
467 fstrsolid%CONT_FTRAC(:) = 0.d0
468 end if
469
470 end subroutine
471
472 subroutine setup_contact_elesurf_for_area( cstep, hecMESH, fstrSOLID )
473 integer(kind=kint), intent(in) :: cstep
474 type( hecmwst_local_mesh ), intent(in) :: hecMESH
475 type(fstr_solid), intent(inout) :: fstrSOLID
476
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(:)
480
481 if( associated(fstrsolid%CONT_SGRP_ID) ) deallocate(fstrsolid%CONT_SGRP_ID)
482
483 allocate(cdef_surf(l_max_elem_surf,hecmesh%n_elem))
484 cdef_surf(:,:) = .false.
485
486 ! label contact defined surfaces
487 n_cdsurfs = 0
488 do i=1, size(fstrsolid%contacts)
489 !grpid = fstrSOLID%contacts(i)%group
490 !if( .not. fstr_isContactActive( fstrSOLID, grpid, cstep ) ) then
491 ! call clear_contact_state(fstrSOLID%contacts(i)); cycle
492 !endif
493
494 do k=1,2 !slave,master
495 if( k==1 ) then !slave
496 sgrp_id = fstrsolid%contacts(i)%surf_id1_sgrp
497 else if( k==2 ) then !master
498 sgrp_id = fstrsolid%contacts(i)%surf_id2
499 end if
500
501 if( sgrp_id <= 0 ) cycle
502
503 is = hecmesh%surf_group%grp_index(sgrp_id-1) + 1
504 ie = hecmesh%surf_group%grp_index(sgrp_id )
505 do j=is,ie
506 eid = hecmesh%surf_group%grp_item(2*j-1)
507 sid = hecmesh%surf_group%grp_item(2*j)
508 ! only internal and boundary element should be added
509 if( .not. cdef_surf(sid,eid) ) n_cdsurfs = n_cdsurfs + 1
510 cdef_surf(sid,eid) = .true.
511 enddo
512 end do
513 enddo
514
515 !gather info of contact defined surfaces
516 allocate(fstrsolid%CONT_SGRP_ID(2*n_cdsurfs))
517 n_cdsurfs = 0
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
524 endif
525 end do
526 end do
527 deallocate(cdef_surf)
528
529 end subroutine
530
531 subroutine calc_contact_area( hecMESH, fstrSOLID, flag )
532 type( hecmwst_local_mesh ), intent(in) :: hecmesh
533 type(fstr_solid), intent(inout) :: fstrSOLID
534 integer(kind=kint), intent(in) :: flag
535
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)
541
542 fstrsolid%CONT_AREA(:) = 0.d0
543
544 if( .not. associated(fstrsolid%CONT_SGRP_ID) ) return
545
546 allocate(coord(ndof*hecmesh%n_node))
547 do i=1,ndof*hecmesh%n_node
548 coord(i) = hecmesh%node(i)+fstrsolid%unode(i)
549 end do
550 if( flag == 1 ) then
551 do i=1,ndof*hecmesh%n_node
552 coord(i) = coord(i)+fstrsolid%dunode(i)
553 end do
554 end if
555
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 )
559
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)
564
565 do i=1,nn
566 idx = ndof*(ndlocal(i)-1)
567 ecoord(1:ndof,i) = coord(idx+1:idx+ndof)
568 end do
569
570 call calc_nodalarea_surfelement( etype, nn, ecoord, sid, vect )
571
572 do i=1,nn
573 idx = ndlocal(i)
574 fstrsolid%CONT_AREA(idx) = fstrsolid%CONT_AREA(idx) + vect(i)
575 end do
576
577 end do
578
579 deallocate(coord)
580 end subroutine
581
582 subroutine calc_nodalarea_surfelement( etype, nn, ecoord, sid, vect )
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(:)
588
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)
594
595 vect(:) = 0.d0
596
597 call getsubface( etype, sid, stype, nod )
598 nsur = getnumberofnodes( stype )
599 do i=1,nsur
600 scoord(1:ndof,i) = ecoord(1:ndof,nod(i))
601 end do
602
603 area = 0.d0
604 do ig0=1,numofquadpoints( stype )
605 call getquadpoint( stype, ig0, localcoord(1:2) )
606 call getshapefunc( stype, localcoord(1:2), h(1:nsur) )
607
608 wg=getweight( stype, ig0 )
609 ! normal = dx/dr_1 \times dx/dr_2
610 normal(1:3) = surfacenormal( stype, nsur, localcoord(1:2), scoord(1:ndof,1:nsur) )
611 area = area + dsqrt(dot_product(normal,normal))*wg
612 enddo
613 area = area/dble(nsur)
614 do i=1,nsur
615 vect(nod(i)) = area
616 end do
617
618 end subroutine
619
620end module mcontact
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.
Definition: hecmw.f90:6
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint), parameter kcaslagrange
contact analysis algorithm
Definition: m_fstr.f90:53
real(kind=kreal) dt
ANALYSIS CONTROL for NLGEOM and HEAT.
Definition: m_fstr.f90:123
integer(kind=kint), parameter kcaalagrange
Definition: m_fstr.f90:54
logical function fstr_iscontactactive(fstrsolid, nbc, cstep)
Definition: m_fstr.f90:1021
This module provides functions to calcualte contact stiff matrix.
Definition: fstr_contact.f90:6
logical function fstr_is_matrixstructure_changed(infoctchange)
subroutine fstr_update_contact_tangentforce(fstrsolid)
Update tangent force.
subroutine fstr_set_contact_active(a)
subroutine fstr_write_mpc(file, mpcs)
Print out mpc conditions.
subroutine fstr_scan_contact_state(cstep, sub_step, cont_step, dt, ctalgo, hecmesh, fstrsolid, infoctchange, b)
Scanning contact state.
subroutine fstr_scan_contact_state_exp(cstep, hecmesh, fstrsolid, infoctchange)
Scanning contact state.
subroutine fstr_contactbc(iter, hecmesh, hecmat, fstrsolid)
Introduce contact stiff into global stiff matrix or mpc conditions into hecMESH.
subroutine fstr_contact2mpc(contacts, mpcs)
Contact states to equation conditions.
real(kind=kreal), dimension(2), save bakgnt
1:current avarage penetration; 2:current relative tangent displacement!
subroutine setup_contact_elesurf_for_area(cstep, hecmesh, fstrsolid)
logical function fstr_is_contact_conv(ctalgo, infoctchange, hecmesh)
real(kind=kreal), save cdotp
mu=cdotp*maxval
subroutine fstr_del_contactmpc(mpcs)
Delete mpcs derived from contact conditions.
subroutine calc_contact_area(hecmesh, fstrsolid, flag)
integer(kind=kint), save n_contact_mpc
subroutine fstr_set_contact_penalty(maxv)
real(kind=kreal), save mut
penalty along tangent direction
subroutine print_contatct_pair(file, pair)
Write out the contact definition read from mesh file.
real(kind=kreal), dimension(2), save gnt
1:current avarage penetration; 2:current relative tangent displacement
real(kind=kreal), save mu
penalty, default value
subroutine fstr_update_contact_multiplier(hecmesh, fstrsolid, ctchanged)
Update lagrangian multiplier.
real(kind=kreal), save cgt
convergent condition of relative tangent disp
subroutine fstr_update_contact0(hecmesh, fstrsolid, b)
Update lagrangian multiplier.
logical function fstr_is_contact_active()
subroutine calc_nodalarea_surfelement(etype, nn, ecoord, sid, vect)
subroutine initialize_contact_output_vectors(fstrsolid, hecmat)
real(kind=kreal), save cgn
convergent condition of penetration
This module manage the data structure for contact calculation.
subroutine set_contact_state_vector(contact, dt, relvel_vec, state_vec)
This subroutine setup contact output nodal vectors.
subroutine scan_contact_state(flag_ctalgo, contact, currpos, currdisp, ndforce, infoctchange, nodeid, elemid, is_init, active, mu, b)
This subroutine update contact states, which include.
subroutine update_contact_multiplier(contact, coord, disp, ddisp, fcoeff, mu, mut, gnt, ctchanged)
This subroutine update lagrangian multiplier and the distance between contacting nodes.
subroutine clear_contact_state(contact)
Reset contact state all to free.
subroutine update_contact_tangentforce(contact)
subroutine calcu_contact_force0(contact, coord, disp, ddisp, fcoeff, mu, mut, b)
This subroutine update contact condition as follows:
subroutine scan_contact_state_exp(contact, currpos, currdisp, infoctchange, nodeid, elemid, is_init, active)
This subroutine update contact states, which include.
Structure to includes all info needed by contact calculation.