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