FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 !-------------------------------------------------------------------------------
6 module mcontact
7 
8  use mcontactdef
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 
29 contains
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()
52  fstr_is_contact_active = active
53  end function
54 
55  subroutine fstr_set_contact_active( a )
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 ) &
69  fstr_is_contact_conv = .true.
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
83  fstr_is_contact_conv = .true.
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 
620 end 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
logical function fstr_iscontactactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1021
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
This module provides functions to calcualte contact stiff matrix.
Definition: fstr_contact.f90:6
subroutine fstr_set_contact_active(a)
subroutine fstr_write_mpc(file, mpcs)
Print out mpc conditions.
subroutine fstr_update_contact_multiplier(hecMESH, fstrSOLID, ctchanged)
Update lagrangian multiplier.
logical function fstr_is_matrixstructure_changed(infoCTChange)
logical function fstr_is_contact_conv(ctAlgo, infoCTChange, hecMESH)
subroutine fstr_contact2mpc(contacts, mpcs)
Contact states to equation conditions.
subroutine fstr_contactbc(iter, hecMESH, hecMAT, fstrSOLID)
Introduce contact stiff into global stiff matrix or mpc conditions into hecMESH.
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)
real(kind=kreal), save cdotp
mu=cdotp*maxval
subroutine fstr_del_contactmpc(mpcs)
Delete mpcs derived from contact conditions.
subroutine fstr_scan_contact_state(cstep, sub_step, cont_step, dt, ctAlgo, hecMESH, fstrSOLID, infoCTChange, B)
Scanning contact state.
integer(kind=kint), save n_contact_mpc
subroutine initialize_contact_output_vectors(fstrSOLID, hecMAT)
subroutine fstr_set_contact_penalty(maxv)
subroutine fstr_scan_contact_state_exp(cstep, hecMESH, fstrSOLID, infoCTChange)
Scanning contact state.
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
real(kind=kreal), save cgt
convergent condition of relative tangent disp
subroutine calc_contact_area(hecMESH, fstrSOLID, flag)
subroutine fstr_update_contact_tangentforce(fstrSOLID)
Update tangent force.
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)
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 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 scan_contact_state_exp(contact, currpos, currdisp, infoCTChange, nodeID, elemID, is_init, active)
This subroutine update contact states, which include.
subroutine calcu_contact_force0(contact, coord, disp, ddisp, fcoeff, mu, mut, B)
This subroutine update contact condition as follows:
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.
Structure to includes all info needed by contact calculation.