1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5!> \brief This module provides functions of reconstructing
6!         stiffness matrix structure for the contact analysis
7!         employing standard Lagrange multiplier algorithm
8
9module fstr_matrix_con_contact
10
11  use m_fstr
12  use elementInfo
13  use mContact
14
15  implicit none
16
17  !> Structure for defining stiffness matrix structure
18  type nodeRelated
19    integer(kind=kint)          :: num_node = 0, num_lagrange = 0 !< total number of related nodes and Lagrange multipliers
20    integer(kind=kint), pointer :: id_node(:) => null() !< list of related nodes
21    integer(kind=kint), pointer :: id_lagrange(:) => null() !< list of related Lagrange multipliers
22  end type
23
24  !> Structure for Lagrange multiplier-related part of stiffness matrix
25  !> (Lagrange multiplier-related matrix)
26  type fstrST_matrix_contact_lagrange
27    integer(kind=kint) :: num_lagrange = 0 !< total number of Lagrange multipliers
28    integer(kind=kint) :: numL_lagrange = 0, numU_lagrange = 0 !< node-based number of non-zero items in lower triangular half of matrix
29    !< node-based number of non-zero items in upper triangular half of matrix
30    integer(kind=kint), pointer  :: indexL_lagrange(:) => null(), &
31      indexU_lagrange(:) => null() !< node-based index of first non-zero item of each row
32    integer(kind=kint), pointer  :: itemL_lagrange(:) => null(), &
33      itemU_lagrange(:) => null() !< node-based column number of non-zero items
34    real(kind=kreal),    pointer  :: AL_lagrange(:) => null(), &
35      AU_lagrange(:) => null() !< values of non-zero items
36    real(kind=kreal),    pointer  :: Lagrange(:) => null() !< values of Lagrange multipliers
37  end type fstrST_matrix_contact_lagrange
38
39  integer(kind=kint), save         :: NPL_org, NPU_org !< original number of non-zero items
40  type(nodeRelated), pointer, save :: list_nodeRelated_org(:) => null() !< original structure of matrix
41
42  type(nodeRelated), pointer       :: list_nodeRelated(:) => null() !< current structure of matrix
43
44  logical                          :: permission = .false.
45
46  private                          :: insert_lagrange, insert_node, find_locationINarray, reallocate_memory
47
48contains
49
50
51  !> \brief This subroutine saves original matrix structure constructed originally by hecMW_matrix
52  subroutine fstr_save_originalMatrixStructure(hecMAT)
53
54    type(hecmwST_matrix) :: hecMAT !< type hecmwST_matrix
55
56    integer(kind=kint)   :: numL, numU, num_nodeRelated !< original number of nodes related to each node
57    integer(kind=kint)   :: i, j
58    integer(kind=kint)   :: ierr
59
60    NPL_org = hecMAT%NPL; NPU_org = hecMAT%NPU
61
62    if( associated(list_nodeRelated_org) ) return
63    allocate(list_nodeRelated_org(hecMAT%NP), stat=ierr)
64    if( ierr /= 0) stop " Allocation error, list_nodeRelated_org "
65
66    do i = 1, hecMAT%NP
67
68      numL = hecMAT%indexL(i) - hecMAT%indexL(i-1)
69      numU = hecMAT%indexU(i) - hecMAT%indexU(i-1)
70
71      num_nodeRelated = numL + numU + 1
72
73      allocate(list_nodeRelated_org(i)%id_node(num_nodeRelated), stat=ierr)
74      if( ierr /= 0) stop " Allocation error, list_nodeRelated_org%id_node "
75
76      list_nodeRelated_org(i)%num_node = num_nodeRelated
77
78      do j = 1, numL
79        list_nodeRelated_org(i)%id_node(j) = hecMAT%itemL(hecMAT%indexL(i-1)+j)
80      enddo
81      list_nodeRelated_org(i)%id_node(numL+1) = i
82      do j = 1, numU
83        list_nodeRelated_org(i)%id_node(numL+1+j) = hecMAT%itemU(hecMAT%indexU(i-1)+j)
84      enddo
85
86    enddo
87
88  end subroutine fstr_save_originalMatrixStructure
89
90  !> \brief this subroutine reconstructs node-based (stiffness) matrix structure
91  !> \corresponding to contact state
92  subroutine fstr_mat_con_contact(cstep,hecMAT,fstrSOLID,fstrMAT,infoCTChange,conMAT)
93
94    integer(kind=kint)                   :: cstep !< current loading step
95    type(hecmwST_matrix)                 :: hecMAT !< type hecmwST_matrix
96    type(fstr_solid)                     :: fstrSOLID !< type fstr_solid
97    type(fstrST_matrix_contact_lagrange) :: fstrMAT !< type fstrST_matrix_contact_lagrange
98    type(fstr_info_contactChange)        :: infoCTChange !< type fstr_contactChange
99
100    integer(kind=kint)                   :: num_lagrange !< number of Lagrange multipliers
101    integer(kind=kint)                   :: countNon0LU_node, countNon0LU_lagrange !< counter of node-based number of non-zero items
102    integer(kind=kint)                   :: numNon0_node, numNon0_lagrange !< node-based number of displacement-related non-zero items in half of the matrix
103    !< node-based number of Lagrange multiplier-related non-zero items in half of the matrix
104    type (hecmwST_matrix),optional          :: conMAT
105
106    num_lagrange = infoCTChange%contactNode_current
107    fstrMAT%num_lagrange = num_lagrange
108
109    ! Get original list of related nodes
110    call getOriginalListOFrelatedNodes(hecMAT%NP,num_lagrange)
111
112    ! Construct new list of related nodes and Lagrange multipliers
113    countNon0LU_node = NPL_org + NPU_org
114    countNon0LU_lagrange = 0
115    if( fstr_is_contact_active() ) &
116      call getNewListOFrelatednodesANDLagrangeMultipliers(cstep,hecMAT%NP,fstrSOLID,countNon0LU_node,countNon0LU_lagrange)
117
118    ! Construct new matrix structure(hecMAT&fstrMAT)
119    numNon0_node = countNon0LU_node/2
120    numNon0_lagrange = countNon0LU_lagrange/2
121    !     ----  For Parallel Contact with Multi-Partition Domains
122    if(paraContactFlag.and.present(conMAT)) then
123      call constructNewMatrixStructure(hecMAT,fstrMAT,numNon0_node,numNon0_lagrange,conMAT)
124    else
125      call constructNewMatrixStructure(hecMAT,fstrMAT,numNon0_node,numNon0_lagrange)
126    endif
127
128    ! Copy Lagrange multipliers
129    if( fstr_is_contact_active() ) &
130      call fstr_copy_lagrange_contact(fstrSOLID,fstrMAT)
131
132  end subroutine fstr_mat_con_contact
133
134  !> Get original list of related nodes
135  subroutine getOriginalListOFrelatedNodes(np,num_lagrange)
136
137    integer(kind=kint)            :: np, num_lagrange !< total number of nodes
138    integer(kind=kint)            :: num_nodeRelated_org !< original number of nodes related to each node
139    integer(kind=kint)            :: i, ierr
140
141    allocate(list_nodeRelated(np+num_lagrange), stat=ierr)
142    if( ierr /= 0) stop " Allocation error, list_nodeRelated "
143
144    do i = 1, np  !hecMAT%NP
145      num_nodeRelated_org = list_nodeRelated_org(i)%num_node
146      allocate(list_nodeRelated(i)%id_node(num_nodeRelated_org), stat=ierr)
147      if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_node "
148      list_nodeRelated(i)%num_node = num_nodeRelated_org
149      list_nodeRelated(i)%id_node(1:num_nodeRelated_org) = list_nodeRelated_org(i)%id_node(1:num_nodeRelated_org)
150    enddo
151
152    if( fstr_is_contact_active() ) then
153      do i = np+1, np+num_lagrange  !hecMAT%NP+1, hecMAT%NP+num_lagrange
154        allocate(list_nodeRelated(i)%id_lagrange(5), stat=ierr)
155        if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_lagrange "
156        list_nodeRelated(i)%num_lagrange = 0
157        list_nodeRelated(i)%id_lagrange = 0
158      enddo
159    endif
160
161  end subroutine getOriginalListOFrelatedNodes
162
163
164  !> Construct new list of related nodes and Lagrange multipliers. Here, a procedure similar to HEC_MW is used.
165  subroutine getNewListOFrelatednodesANDLagrangeMultipliers(cstep,np,fstrSOLID,countNon0LU_node,countNon0LU_lagrange)
166
167    type(fstr_solid)              :: fstrSOLID !< type fstr_solid
168
169    integer(kind=kint)            :: cstep !< current loading step
170    integer(kind=kint)            :: np !< total number of nodes
171    integer(kind=kint)            :: countNon0LU_node, countNon0LU_lagrange !< counters of node-based number of non-zero items
172    integer(kind=kint)            :: grpid !< contact pairs group ID
173    integer(kind=kint)            :: count_lagrange !< counter of Lagrange multiplier
174    integer(kind=kint)            :: ctsurf, etype, nnode, ndLocal(l_max_surface_node + 1) !< contants of type tContact
175    integer(kind=kint)            :: i, j, k, l, num, num_nodeRelated_org, ierr
176    real(kind=kreal)              :: fcoeff !< friction coefficient
177
178    count_lagrange = 0
179    do i = 1, size(fstrSOLID%contacts)
180
181      grpid = fstrSOLID%contacts(i)%group
182      if( .not. fstr_isContactActive( fstrSOLID, grpid, cstep ) ) cycle
183
184      fcoeff = fstrSOLID%contacts(i)%fcoeff
185
186      do j = 1, size(fstrSOLID%contacts(i)%slave)
187
188        if( fstrSOLID%contacts(i)%states(j)%state == CONTACTFREE ) cycle
189        ctsurf = fstrSOLID%contacts(i)%states(j)%surface
190        etype = fstrSOLID%contacts(i)%master(ctsurf)%etype
191        if( etype/=fe_tri3n .and. etype/=fe_quad4n ) &
192          stop " ##Error: This element type is not supported in contact analysis !!! "
193        nnode = size(fstrSOLID%contacts(i)%master(ctsurf)%nodes)
194        ndLocal(1) = fstrSOLID%contacts(i)%slave(j)
195        ndLocal(2:nnode+1) = fstrSOLID%contacts(i)%master(ctsurf)%nodes(1:nnode)
196
197        count_lagrange = count_lagrange + 1
198
199        do k = 1, nnode+1
200
201          if( .not. associated(list_nodeRelated(ndLocal(k))%id_lagrange) )then
202            num = 10
203            !             if( k == 1 ) num = 1
204            allocate(list_nodeRelated(ndLocal(k))%id_lagrange(num),stat=ierr)
205            if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_lagrange "
206            list_nodeRelated(ndLocal(k))%num_lagrange = 0
207            list_nodeRelated(ndLocal(k))%id_lagrange = 0
208          endif
209
210          if( fcoeff /= 0.0d0 ) then
211            num_nodeRelated_org = list_nodeRelated_org(ndLocal(k))%num_node
212            if( list_nodeRelated(ndLocal(k))%num_node == num_nodeRelated_org )then
213              num = 10
214              if(k==1) num = 4
215              call reallocate_memory(num,list_nodeRelated(ndLocal(k)))
216            endif
217          endif
218
219          call insert_lagrange(k,count_lagrange,list_nodeRelated(ndLocal(k)),countNon0LU_lagrange)
220
221          do l = k, nnode+1
222            if( fcoeff /= 0.0d0 ) then
223              if( k /= l) then
224                num_nodeRelated_org = list_nodeRelated_org(ndLocal(k))%num_node
225                call insert_node(ndLocal(l),list_nodeRelated(ndLocal(k)),countNon0LU_node)
226                num_nodeRelated_org = list_nodeRelated_org(ndLocal(l))%num_node
227                if( list_nodeRelated(ndLocal(l))%num_node == num_nodeRelated_org )then
228                  num = 10
229                  call reallocate_memory(num,list_nodeRelated(ndLocal(l)))
230                endif
231                call insert_node(ndLocal(k),list_nodeRelated(ndLocal(l)),countNon0LU_node)
232              endif
233            endif
234
235            if(k == 1) &
236              call insert_lagrange(0,ndLocal(l),list_nodeRelated(np+count_lagrange),countNon0LU_lagrange)
237
238          enddo
239
240        enddo
241
242      enddo
243
244    enddo
245
246  end subroutine getNewListOFrelatednodesANDLagrangeMultipliers
247
248
249  !> Construct new stiffness matrix structure
250  subroutine constructNewMatrixStructure(hecMAT,fstrMAT,numNon0_node,numNon0_lagrange,conMAT)
251
252    type(hecmwST_matrix)                 :: hecMAT !< type hecmwST_matrix
253    type(fstrST_matrix_contact_lagrange) :: fstrMAT !< type fstrST_matrix_contact_lagrange
254
255    integer(kind=kint)                   :: numNon0_node, numNon0_lagrange !< node-based number of non-zero items in half of the matrix
256    integer(kind=kint)                   :: countNon0L_node, countNon0U_node, countNon0U_lagrange, countNon0L_lagrange !< counters of node-based number ofnon-zero items
257    integer(kind=kint)                   :: i, j, ierr
258    integer(kind=kint)                   :: numI_node, numI_lagrange
259    integer(kind=kint)                   :: ndof, nn
260    type(hecmwST_matrix), optional       :: conMAT
261
262    !     ----  For Parallel Contact with Multi-Partition Domains
263    if(paraContactFlag.and.present(conMAT)) then
264      conMAT%N  = hecMAT%N
265      conMAT%NP = hecMAT%NP
266      conMAT%ndof = hecMAT%ndof
267      if(associated(conMAT%indexL).and.associated(conMAT%indexU))deallocate(conMAT%indexL,conMAT%indexU)
268      allocate(conMAT%indexL(0:conMAT%NP), conMAT%indexU(0:conMAT%NP), stat=ierr)
269      if ( ierr /= 0) stop " Allocation error, conMAT%indexL-conMAT%indexU "
270      conMAT%indexL = 0 ; conMAT%indexU = 0
271      if(associated(conMAT%itemL).and.associated(conMAT%itemU))deallocate(conMAT%itemL,conMAT%itemU)
272      allocate(conMAT%itemL(numNon0_node), conMAT%itemU(numNon0_node), stat=ierr)
273      if ( ierr /= 0) stop " Allocation error, conMAT%itemL-conMAT%itemU "
274      conMAT%itemL = 0 ; conMAT%itemU = 0
275      !
276      conMAT%NPL = numNon0_node
277      conMAT%NPU = numNon0_node
278    endif
279
280    if(associated(hecMAT%indexL).and.associated(hecMAT%indexU))deallocate(hecMAT%indexL,hecMAT%indexU)
281    allocate(hecMAT%indexL(0:hecMAT%NP), hecMAT%indexU(0:hecMAT%NP), stat=ierr)
282    if ( ierr /= 0) stop " Allocation error, hecMAT%indexL-hecMAT%indexU "
283    hecMAT%indexL = 0 ; hecMAT%indexU = 0
284    if(associated(hecMAT%itemL).and.associated(hecMAT%itemU))deallocate(hecMAT%itemL,hecMAT%itemU)
285    allocate(hecMAT%itemL(numNon0_node), hecMAT%itemU(numNon0_node), stat=ierr)
286    if ( ierr /= 0) stop " Allocation error, hecMAT%itemL-hecMAT%itemU "
287    hecMAT%itemL = 0 ; hecMAT%itemU = 0
288
289    if(associated(fstrMAT%indexL_lagrange).and.associated(fstrMAT%indexU_lagrange)) &
290      deallocate(fstrMAT%indexL_lagrange,fstrMAT%indexU_lagrange)
291    if(associated(fstrMAT%itemL_lagrange).and.associated(fstrMAT%itemU_lagrange)) &
292      deallocate(fstrMAT%itemL_lagrange,fstrMAT%itemU_lagrange)
293    if( fstr_is_contact_active() ) then
294      allocate(fstrMAT%indexL_lagrange(0:fstrMAT%num_lagrange), fstrMAT%indexU_lagrange(0:hecMAT%NP), stat=ierr)
295      if ( ierr /= 0) stop " Allocation error, fstrMAT%indexL_lagrange-fstrMAT%indexU_lagrange "
296      fstrMAT%indexL_lagrange = 0 ; fstrMAT%indexU_lagrange = 0
297      allocate(fstrMAT%itemL_lagrange(numNon0_lagrange), fstrMAT%itemU_lagrange(numNon0_lagrange), stat=ierr)
298      if ( ierr /= 0) stop " Allocation error, fstrMAT%itemL_lagrange-fstrMAT%itemU_lagrange "
299      fstrMAT%itemL_lagrange = 0 ; fstrMAT%itemU_lagrange = 0
300    endif
301
302    hecMAT%NPL = numNon0_node
303    hecMAT%NPU = numNon0_node
304
305    fstrMAT%numL_lagrange = numNon0_lagrange
306    fstrMAT%numU_lagrange = numNon0_lagrange
307
308    countNon0L_node = 0
309    countNon0U_node = 0
310    countNon0U_lagrange = 0
311    do i = 1, hecMAT%NP
312
313      list_nodeRelated(i)%num_node = count(list_nodeRelated(i)%id_node /= 0)
314      numI_node = list_nodeRelated(i)%num_node
315      if( fstr_is_contact_active() ) &
316        numI_lagrange = list_nodeRelated(i)%num_lagrange
317
318      do j = 1, numI_node
319        if( list_nodeRelated(i)%id_node(j) < i ) then
320          countNon0L_node = countNon0L_node + 1
321          hecMAT%itemL(countNon0L_node) = list_nodeRelated(i)%id_node(j)
322        elseif( list_nodeRelated(i)%id_node(j) > i ) then
323          countNon0U_node = countNon0U_node + 1
324          hecMAT%itemU(countNon0U_node) = list_nodeRelated(i)%id_node(j)
325        endif
326      enddo
327      hecMAT%indexL(i) = countNon0L_node
328      hecMAT%indexU(i) = countNon0U_node
329
330      if( fstr_is_contact_active() ) then
331        do j = 1, numI_lagrange
332          countNon0U_lagrange = countNon0U_lagrange + 1
333          fstrMAT%itemU_lagrange(countNon0U_lagrange) = list_nodeRelated(i)%id_lagrange(j)
334        enddo
335        fstrMAT%indexU_lagrange(i) = countNon0U_lagrange
336      endif
337
338      deallocate(list_nodeRelated(i)%id_node)
339      if(associated(list_nodeRelated(i)%id_lagrange)) deallocate(list_nodeRelated(i)%id_lagrange)
340
341    end do
342
343    !     ----  For Parallel Contact with Multi-Partition Domains
344    if(paraContactFlag.and.present(conMAT)) then
345      conMAT%itemL(:)   = hecMAT%itemL(:)
346      conMAT%indexL(:)  = hecMAT%indexL(:)
347      conMAT%itemU(:)   = hecMAT%itemU(:)
348      conMAT%indexU(:)  = hecMAT%indexU(:)
349    endif
350
351    if( fstr_is_contact_active() ) then
352      countNon0L_lagrange = 0
353      do i = 1, fstrMAT%num_lagrange
354        numI_lagrange = list_nodeRelated(hecMAT%NP+i)%num_lagrange
355        do j = 1, numI_lagrange
356          countNon0L_lagrange = countNon0L_lagrange + 1
357          fstrMAT%itemL_lagrange(countNon0L_lagrange) = list_nodeRelated(hecMAT%NP+i)%id_lagrange(j)
358        enddo
359        fstrMAT%indexL_lagrange(i) = countNon0L_lagrange
360        deallocate(list_nodeRelated(hecMAT%NP+i)%id_lagrange)
361      enddo
362    endif
363
364    deallocate(list_nodeRelated)
365
366    ndof = hecMAT%NDOF
367    nn = ndof*ndof
368    if(associated(hecMAT%AL)) deallocate(hecMAT%AL)
369    allocate(hecMAT%AL(nn*hecMAT%NPL), stat=ierr)
370    if ( ierr /= 0 ) stop " Allocation error, hecMAT%AL "
371    hecMAT%AL = 0.0D0
372
373    if(associated(hecMAT%AU)) deallocate(hecMAT%AU)
374    allocate(hecMAT%AU(nn*hecMAT%NPU), stat=ierr)
375    if ( ierr /= 0 ) stop " Allocation error, hecMAT%AU "
376    hecMAT%AU = 0.0D0
377
378    if(associated(fstrMAT%AL_lagrange)) deallocate(fstrMAT%AL_lagrange)
379    if(associated(fstrMAT%AU_lagrange)) deallocate(fstrMAT%AU_lagrange)
380    if(associated(fstrMAT%Lagrange)) deallocate(fstrMAT%Lagrange)
381
382    if( fstr_is_contact_active() ) then
383      allocate(fstrMAT%AL_lagrange(ndof*fstrMAT%numL_lagrange), stat=ierr)
384      if ( ierr /= 0 ) stop " Allocation error, fstrMAT%AL_lagrange "
385      fstrMAT%AL_lagrange = 0.0D0
386      allocate(fstrMAT%AU_lagrange(ndof*fstrMAT%numU_lagrange), stat=ierr)
387      if ( ierr /= 0 ) stop " Allocation error, fstrMAT%AU_lagrange "
388      fstrMAT%AU_lagrange = 0.0D0
389      allocate(fstrMAT%Lagrange(fstrMAT%num_lagrange))
390      fstrMAT%Lagrange = 0.0D0
391    endif
392
393    if(associated(hecMAT%B)) deallocate(hecMAT%B)
394    allocate(hecMAT%B(hecMAT%NP*ndof+fstrMAT%num_lagrange))
395    hecMAT%B = 0.0D0
396
397    if(associated(hecMAT%X)) deallocate(hecMAT%X)
398    allocate(hecMAT%X(hecMAT%NP*ndof+fstrMAT%num_lagrange))
399    hecMAT%X = 0.0D0
400
401    if(associated(hecMAT%D)) deallocate(hecMAT%D)
402    allocate(hecMAT%D(hecMAT%NP*ndof**2+fstrMAT%num_lagrange))
403    hecMAT%D = 0.0D0
404    !
405    !     ----  For Parallel Contact with Multi-Partition Domains
406    if(paraContactFlag.and.present(conMAT)) then
407      if(associated(conMAT%AL)) deallocate(conMAT%AL)
408      allocate(conMAT%AL(nn*conMAT%NPL), stat=ierr)
409      if ( ierr /= 0 ) stop " Allocation error, conMAT%AL "
410      conMAT%AL = 0.0D0
411
412      if(associated(conMAT%AU)) deallocate(conMAT%AU)
413      allocate(conMAT%AU(nn*conMAT%NPU), stat=ierr)
414      if ( ierr /= 0 ) stop " Allocation error, conMAT%AU "
415      conMAT%AU = 0.0D0
416
417      if(associated(conMAT%B)) deallocate(conMAT%B)
418      allocate(conMAT%B(conMAT%NP*ndof+fstrMAT%num_lagrange))
419      conMAT%B = 0.0D0
420
421      if(associated(conMAT%X)) deallocate(conMAT%X)
422      allocate(conMAT%X(conMAT%NP*ndof+fstrMAT%num_lagrange))
423      conMAT%X = 0.0D0
424
425      if(associated(conMAT%D)) deallocate(conMAT%D)
426      allocate(conMAT%D(conMAT%NP*ndof**2+fstrMAT%num_lagrange))
427      conMAT%D = 0.0D0
428    endif
429
430  end subroutine ConstructNewMatrixStructure
431
432
433  !> Insert a Lagrange multiplier in list of related Lagrange multipliers
434  subroutine insert_lagrange(i,id_lagrange,list_node,countNon0_lagrange)
435
436    type(nodeRelated)  :: list_node !< type nodeRelated
437    integer(kind=kint) :: i, id_lagrange !< local number of node in current contact pair
438    !< Lagrange multiplier ID
439    integer(kind=kint) :: countNon0_lagrange !< counter of node-based number of non-zero items
440    !< in Lagrange multiplier-related matrix
441    integer(kind=kint) :: ierr, num_lagrange, location
442    integer(kind=kint) :: id_lagrange_save(1000)
443
444    character(len=1)   :: answer
445
446    ierr = 0
447
448    num_lagrange = count(list_node%id_lagrange /= 0 )
449
450    !      if( i == 1 .and. num_lagrange /= 0) return
451    if( i == 1 .and. num_lagrange /= 0 .and. .not. permission) then
452      1  write(*,*) '##Error: node is both slave and master node simultaneously !'
453      write(*,*) '         Please check contact surface definition !'
454      write(*,'(''          Do you want to continue(y/n)):'',$)')
455      read(*,'(A1)',err=1) answer
456      if(answer == 'Y' .OR. answer == 'y')then
457        permission = .true.
458      else
459        stop
460      endif
461    endif
462
463    if (num_lagrange == 0)then
464      list_node%num_lagrange = 1
465      list_node%id_lagrange(1) = id_lagrange
466      countNon0_lagrange = countNon0_lagrange + 1
467    else
468      id_lagrange_save(1:num_lagrange) = list_node%id_lagrange(1:num_lagrange)
469      location = find_locationINarray(id_lagrange,num_lagrange,list_node%id_lagrange)
470      if(location /= 0)then
471        num_lagrange = num_lagrange + 1
472        if( num_lagrange > size(list_node%id_lagrange)) then
473          deallocate(list_node%id_lagrange)
474          allocate(list_node%id_lagrange(num_lagrange),stat=ierr)
475          if( ierr /= 0 ) stop " Allocation error, list_nodeRelated%id_lagrange "
476        endif
477        list_node%num_lagrange = num_lagrange
478        list_node%id_lagrange(location) = id_lagrange
479        if(location /= 1) list_node%id_lagrange(1:location-1) = id_lagrange_save(1:location-1)
480        if(location /= num_lagrange) list_node%id_lagrange(location+1:num_lagrange) = id_lagrange_save(location:num_lagrange-1)
481        countNon0_lagrange = countNon0_lagrange + 1
482      endif
483    endif
484
485  end subroutine insert_lagrange
486
487  !> Insert a node in list of related nodes
488  subroutine insert_node(id_node,list_node,countNon0_node)
489
490    type(nodeRelated)  :: list_node !< type nodeRelated
491    integer(kind=kint) :: id_node !< local number of node in current contact pair
492    !< global number of node
493    integer(kind=kint) :: countNon0_node !< counter of node-based number of non-zero items in displacement-related matrix
494    integer(kind=kint) :: ierr, num_node, location
495    integer(kind=kint) :: id_node_save(1000)
496
497    ierr = 0
498
499    num_node = list_node%num_node
500
501    id_node_save(1:num_node) = list_node%id_node(1:num_node)
502    location = find_locationINarray(id_node,num_node,list_node%id_node)
503    if(location /= 0)then
504      num_node = num_node + 1
505      if( num_node > size(list_node%id_node)) then
506        deallocate(list_node%id_node)
507        allocate(list_node%id_node(num_node),stat=ierr)
508        if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_node "
509      endif
510      list_node%num_node = num_node
511      list_node%id_node(location) = id_node
512      if(location /= 1) list_node%id_node(1:location-1) = id_node_save(1:location-1)
513      if(location /= num_node) list_node%id_node(location+1:num_node) = id_node_save(location:num_node-1)
514      countNon0_node = countNon0_node + 1
515    endif
516
517  end subroutine insert_node
518
519
520  !> Find location of an item in an array by bisection method
521  integer function find_locationINarray(item,n,a)
522
523    integer(kind=kint)           :: item, n !< item to be found; length of array
524    integer(kind=kint), pointer  :: a(:) !< array
525    integer(kind=kint)           :: l, r, m
526
527    find_locationINarray = 0
528
529    l = 1 ; r = n ; m = (l+r)/2
530    if( item == a(l) .or. item == a(r) )then
531      return
532    elseif( item < a(l) )then
533      find_locationINarray = 1
534      return
535    elseif( item > a(r) )then
536      find_locationINarray = n + 1
537      return
538    endif
539
540    do while ( l <= r)
541      if( item > a(m) ) then
542        l = m + 1
543        m = (l + r)/2
544      elseif( item < a(m) ) then
545        r = m - 1
546        m = (l + r)/2
547      elseif( item == a(m) )then
548        return
549      endif
550    enddo
551
552    find_locationINarray = m + 1
553
554  end function find_locationINarray
555
556
557  !> Reallocate memory for list_relatedNodes
558  subroutine reallocate_memory(num,list_node)
559
560    type(nodeRelated)  :: list_node !< type nodeRelated
561    integer(kind=kint) :: num !< length to be added
562    integer(kind=kint) :: num_node_org !< original number of related nodes
563    !< before reallocation
564    integer(kind=kint) :: id_save(1000)
565    integer(kind=kint) :: ierr
566
567    num_node_org = size(list_node%id_node)
568    id_save(1:num_node_org) = list_node%id_node(1:num_node_org)
569    deallocate(list_node%id_node)
570    allocate(list_node%id_node(num_node_org+num),stat=ierr)
571    if( ierr /= 0) stop " reAllocation error, list_nodeRelated%id_node "
572    list_node%id_node = 0
573    list_node%id_node(1:num_node_org) = id_save(1:num_node_org)
574
575  end subroutine reallocate_memory
576
577
578  !> Copy Lagrange multipliers
579  subroutine fstr_copy_lagrange_contact(fstrSOLID,fstrMAT)
580
581    type(fstr_solid)                        :: fstrSOLID                !< type fstr_solid
582    type(fstrST_matrix_contact_lagrange)    :: fstrMAT                  !< fstrST_matrix_contact_lagrange
583    integer (kind=kint)                    :: id_lagrange, i, j
584
585    id_lagrange = 0
586
587    do i = 1, size(fstrSOLID%contacts)
588      do j = 1, size(fstrSOLID%contacts(i)%slave)
589        if( fstrSOLID%contacts(i)%states(j)%state == CONTACTFREE ) cycle
590        id_lagrange = id_lagrange + 1
591        fstrMAT%Lagrange(id_lagrange)=fstrSOLID%contacts(i)%states(j)%multiplier(1)
592      enddo
593    enddo
594
595  end subroutine fstr_copy_lagrange_contact
596
597
598  !> \brief this function judges whether sitiffness matrix is symmetric or not
599  logical function fstr_is_matrixStruct_symmetric(fstrSOLID,hecMESH)
600
601    type(fstr_solid )        :: fstrSOLID
602    type(hecmwST_local_mesh) :: hecMESH
603    integer (kind=kint)      :: is_in_contact
604
605    is_in_contact = 0
606    if( any(fstrSOLID%contacts(:)%fcoeff /= 0.0d0) )  &
607      is_in_contact = 1
608    call hecmw_allreduce_I1(hecMESH, is_in_contact, HECMW_MAX)
609    fstr_is_matrixStruct_symmetric = (is_in_contact == 0)
610
611  end function fstr_is_matrixStruct_symmetric
612
613
614end module fstr_matrix_con_contact
615