1!------------------------------------------------------------------------------- 2! Copyright (c) 2019 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5!> \brief This module contains subroutines for nonlinear explicit dynamic analysis 6 7module fstr_dynamic_nlexplicit 8 use m_fstr 9 use m_static_lib 10 use m_dynamic_output 11 use m_fstr_EIG_setMASS 12 use m_dynamic_mat_ass_bc_ac 13 use m_dynamic_mat_ass_bc 14 use m_dynamic_mat_ass_bc_vl 15 use m_dynamic_mat_ass_load 16 use m_fstr_Update 17 use m_fstr_Restart 18 use fstr_matrix_con_contact 19 use m_dynamic_mat_ass_couple 20 use m_fstr_rcap_io 21 22contains 23 24 !C================================================================C 25 !C-- subroutine fstr_solve_LINEAR_DYNAMIC 26 !C================================================================C 27 subroutine fstr_solve_dynamic_nlexplicit(hecMESH,hecMAT,fstrSOLID,fstrEIG & 28 ,fstrDYN,fstrRESULT,fstrPARAM,infoCTChange & 29 ,fstrCPL, restrt_step_num ) 30 implicit none 31 type(hecmwST_local_mesh) :: hecMESH 32 type(hecmwST_matrix) :: hecMAT 33 type(fstr_eigen) :: fstrEIG 34 type(fstr_solid) :: fstrSOLID 35 type(hecmwST_result_data) :: fstrRESULT 36 type(fstr_param) :: fstrPARAM 37 type(fstr_dynamic) :: fstrDYN 38 type(fstrST_matrix_contact_lagrange) :: fstrMAT !< type fstrST_matrix_contact_lagrange 39 type(fstr_info_contactChange) :: infoCTChange !< fstr_info_contactChange 40 type(fstr_couple) :: fstrCPL !for COUPLE 41 type(hecmwST_matrix), pointer :: hecMATmpc 42 integer(kind=kint), allocatable :: mark(:) 43 integer(kind=kint) :: nnod, ndof, nn, numnp 44 integer(kind=kint) :: i, j, ids, ide, kk 45 integer(kind=kint) :: kkk0, kkk1 46 integer(kind=kint) :: ierror 47 integer(kind=kint) :: iiii5, iexit 48 integer(kind=kint) :: revocap_flag 49 real(kind=kreal), allocatable :: prevB(:) 50 real(kind=kreal) :: a1, a2, a3, b1, b2, b3, c1, c2 51 real(kind=kreal) :: bsize, res 52 real(kind=kreal) :: time_1, time_2 53 integer(kind=kint) :: restrt_step_num 54 real(kind=kreal), parameter :: PI = 3.14159265358979323846D0 55 56 a1 = 0.0d0; a2 = 0.0d0; a3 = 0.0d0; b1 = 0.0d0; b2 = 0.0d0; b3 = 0.0d0 57 c1 = 0.0d0; c2 = 0.0d0 58 59 call hecmw_mpc_mat_init_explicit(hecMESH, hecMAT, hecMATmpc) 60 61 hecMAT%NDOF=hecMESH%n_dof 62 nnod=hecMESH%n_node 63 ndof=hecMAT%NDOF 64 nn=ndof*ndof 65 66 if( fstrPARAM%fg_couple == 1) then 67 if( fstrPARAM%fg_couple_type==5 .or. & 68 fstrPARAM%fg_couple_type==6 ) then 69 allocate( prevB(hecMAT%NP*ndof) ,stat=ierror ) 70 prevB = 0.0d0 71 if( ierror /= 0 ) then 72 write(idbg,*) 'stop due to allocation error <fstr_solve_NONLINEAR_DYNAMIC, prevB>' 73 write(idbg,*) ' rank = ', hecMESH%my_rank,' ierror = ',ierror 74 call flush(idbg) 75 call hecmw_abort( hecmw_comm_get_comm()) 76 endif 77 endif 78 endif 79 80 fstrSOLID%dunode(:) =0.d0 81 82 a1 = 1.d0/fstrDYN%t_delta**2 83 a2 = 1.d0/(2.d0*fstrDYN%t_delta) 84 85 call setMASS(fstrSOLID,hecMESH,hecMAT,fstrEIG) 86 call hecmw_mpc_trans_mass(hecMESH, hecMAT, fstrEIG%mass) 87 88 allocate(mark(hecMAT%NP * hecMAT%NDOF)) 89 call hecmw_mpc_mark_slave(hecMESH, hecMAT, mark) 90 91 do j = 1 ,ndof*nnod 92 fstrDYN%VEC1(j) = (a1 + a2 *fstrDYN%ray_m) * fstrEIG%mass(j) 93 if(mark(j) == 1) fstrDYN%VEC1(j) = 1.d0 94 if(dabs(fstrDYN%VEC1(j)) < 1.0e-20) then 95 if( hecMESH%my_rank == 0 ) then 96 write(*,*) 'stop due to fstrDYN%VEC(j) = 0 , j = ', j 97 write(imsg,*) 'stop due to fstrDYN%VEC(j) = 0 , j = ', j 98 end if 99 call hecmw_abort( hecmw_comm_get_comm()) 100 endif 101 end do 102 103 deallocate(mark) 104 105 !C-- output of initial state 106 if( restrt_step_num == 1 ) then 107 do j = 1 ,ndof*nnod 108 fstrDYN%DISP(j,3) = fstrDYN%DISP(j,1) - fstrDYN%VEL (j,1)/(2.d0*a2) + fstrDYN%ACC (j,1)/ (2.d0*a1) 109 fstrDYN%DISP(j,2) = fstrDYN%DISP(j,1) - fstrDYN%VEL (j,1)/ a2 + fstrDYN%ACC (j,1)/ (2.d0*a1) * 4.d0 110 end do 111 112 call fstr_dynamic_Output(hecMESH, fstrSOLID, fstrDYN, fstrPARAM) 113 call dynamic_output_monit(hecMESH, fstrPARAM, fstrDYN, fstrEIG, fstrSOLID) 114 end if 115 116 if( associated( fstrSOLID%contacts ) ) then 117 call initialize_contact_output_vectors(fstrSOLID,hecMAT) 118 call forward_increment_Lagrange(1,ndof,fstrDYN%VEC1,hecMESH,fstrSOLID,infoCTChange,& 119 & fstrDYN%DISP(:,2),fstrSOLID%ddunode) 120 endif 121 122 do i= restrt_step_num, fstrDYN%n_step 123 124 fstrDYN%i_step = i 125 fstrDYN%t_curr = fstrDYN%t_delta * i 126 127 !C-- mechanical boundary condition 128 call dynamic_mat_ass_load (hecMESH, hecMAT, fstrSOLID, fstrDYN, fstrPARAM) 129 do j=1, hecMESH%n_node* hecMESH%n_dof 130 hecMAT%B(j)=hecMAT%B(j)-fstrSOLID%QFORCE(j) 131 end do 132 133 !C ******************************************************************************** 134 !C for couple analysis 135 if( fstrPARAM%fg_couple == 1 ) then 136 if( fstrPARAM%fg_couple_type==5 .or. & 137 fstrPARAM%fg_couple_type==6 ) then 138 do j = 1, hecMAT%NP * ndof 139 prevB(j) = hecMAT%B(j) 140 enddo 141 endif 142 endif 143 do 144 if( fstrPARAM%fg_couple == 1 ) then 145 if( fstrPARAM%fg_couple_type==1 .or. & 146 fstrPARAM%fg_couple_type==3 .or. & 147 fstrPARAM%fg_couple_type==5 ) call fstr_rcap_get( fstrCPL ) 148 if( fstrPARAM%fg_couple_first /= 0 ) then 149 bsize = dfloat( i ) / dfloat( fstrPARAM%fg_couple_first ) 150 if( bsize > 1.0 ) bsize = 1.0 151 do kkk0 = 1, fstrCPL%coupled_node_n 152 kkk1 = 3 * kkk0 153 fstrCPL%trac(kkk1-2) = bsize * fstrCPL%trac(kkk1-2) 154 fstrCPL%trac(kkk1-1) = bsize * fstrCPL%trac(kkk1-1) 155 fstrCPL%trac(kkk1 ) = bsize * fstrCPL%trac(kkk1 ) 156 enddo 157 endif 158 if( fstrPARAM%fg_couple_window > 0 ) then 159 j = i - restrt_step_num + 1 160 kk = fstrDYN%n_step - restrt_step_num + 1 161 bsize = 0.5*(1.0-cos(2.0*PI*dfloat(j)/dfloat(kk))) 162 do kkk0 = 1, fstrCPL%coupled_node_n 163 kkk1 = 3 * kkk0 164 fstrCPL%trac(kkk1-2) = bsize * fstrCPL%trac(kkk1-2) 165 fstrCPL%trac(kkk1-1) = bsize * fstrCPL%trac(kkk1-1) 166 fstrCPL%trac(kkk1 ) = bsize * fstrCPL%trac(kkk1 ) 167 enddo 168 endif 169 call dynamic_mat_ass_couple( hecMESH, hecMAT, fstrSOLID, fstrCPL ) 170 endif 171 !C ******************************************************************************** 172 173 call hecmw_mpc_trans_rhs(hecMESH, hecMAT, hecMATmpc) 174 175 do j = 1 ,ndof*nnod 176 hecMATmpc%B(j) = hecMATmpc%B(j) + 2.d0*a1* fstrEIG%mass(j) * fstrDYN%DISP(j,1) & 177 + (- a1 + a2 * fstrDYN%ray_m) * fstrEIG%mass(j) * fstrDYN%DISP(j,3) 178 end do 179 180 !C 181 !C-- geometrical boundary condition 182 183 call dynamic_explicit_ass_bc(hecMESH, hecMATmpc, fstrSOLID, fstrDYN) 184 call dynamic_explicit_ass_vl(hecMESH, hecMATmpc, fstrSOLID, fstrDYN) 185 call dynamic_explicit_ass_ac(hecMESH, hecMATmpc, fstrSOLID, fstrDYN) 186 !call dynamic_mat_ass_bc (hecMESH, hecMATmpc, fstrSOLID, fstrDYN, fstrPARAM, fstrMAT) 187 !call dynamic_mat_ass_bc_vl(hecMESH, hecMATmpc, fstrSOLID, fstrDYN, fstrPARAM, fstrMAT) 188 !call dynamic_mat_ass_bc_ac(hecMESH, hecMATmpc, fstrSOLID, fstrDYN, fstrPARAM, fstrMAT) 189 190 ! Finish the calculation 191 do j = 1 ,ndof*nnod 192 hecMATmpc%X(j) = hecMATmpc%B(j) / fstrDYN%VEC1(j) 193 if(dabs(hecMATmpc%X(j)) > 1.0d+5) then 194 if( hecMESH%my_rank == 0 ) then 195 print *, 'Displacement increment too large, please adjust your step size!',i,hecMATmpc%X(j) 196 write(imsg,*) 'Displacement increment too large, please adjust your step size!',i,hecMATmpc%B(j),fstrDYN%VEC1(j) 197 end if 198 call hecmw_abort( hecmw_comm_get_comm()) 199 end if 200 end do 201 call hecmw_mpc_tback_sol(hecMESH, hecMAT, hecMATmpc) 202 203 !C ***************************************************** 204 !C for couple analysis 205 if( fstrPARAM%fg_couple == 1 ) then 206 if( fstrPARAM%fg_couple_type>1 ) then 207 do j=1, fstrCPL%coupled_node_n 208 if( fstrCPL%dof == 3 ) then 209 kkk0 = j*3 210 kkk1 = fstrCPL%coupled_node(j)*3 211 212 fstrCPL%disp (kkk0-2) = hecMAT%X(kkk1-2) 213 fstrCPL%disp (kkk0-1) = hecMAT%X(kkk1-1) 214 fstrCPL%disp (kkk0 ) = hecMAT%X(kkk1 ) 215 216 fstrCPL%velo (kkk0-2) = -b1*fstrDYN%ACC(kkk1-2,1) - b2*fstrDYN%VEL(kkk1-2,1) + & 217 b3*( hecMAT%X(kkk1-2) - fstrDYN%DISP(kkk1-2,1) ) 218 fstrCPL%velo (kkk0-1) = -b1*fstrDYN%ACC(kkk1-1,1) - b2*fstrDYN%VEL(kkk1-1,1) + & 219 b3*( hecMAT%X(kkk1-1) - fstrDYN%DISP(kkk1-1,1) ) 220 fstrCPL%velo (kkk0 ) = -b1*fstrDYN%ACC(kkk1,1) - b2*fstrDYN%VEL(kkk1,1) + & 221 b3*( hecMAT%X(kkk1) - fstrDYN%DISP(kkk1,1) ) 222 fstrCPL%accel(kkk0-2) = -a1*fstrDYN%ACC(kkk1-2,1) - a2*fstrDYN%VEL(kkk1-2,1) + & 223 a3*( hecMAT%X(kkk1-2) - fstrDYN%DISP(kkk1-2,1) ) 224 fstrCPL%accel(kkk0-1) = -a1*fstrDYN%ACC(kkk1-1,1) - a2*fstrDYN%VEL(kkk1-1,1) + & 225 a3*( hecMAT%X(kkk1-1) - fstrDYN%DISP(kkk1-1,1) ) 226 fstrCPL%accel(kkk0 ) = -a1*fstrDYN%ACC(kkk1,1) - a2*fstrDYN%VEL(kkk1,1) + & 227 a3*( hecMAT%X(kkk1) - fstrDYN%DISP(kkk1,1) ) 228 else 229 kkk0 = j*2 230 kkk1 = fstrCPL%coupled_node(j)*2 231 232 fstrCPL%disp (kkk0-1) = hecMAT%X(kkk1-1) 233 fstrCPL%disp (kkk0 ) = hecMAT%X(kkk1 ) 234 235 fstrCPL%velo (kkk0-1) = -b1*fstrDYN%ACC(kkk1-1,1) - b2*fstrDYN%VEL(kkk1-1,1) + & 236 b3*( hecMAT%X(kkk1-1) - fstrDYN%DISP(kkk1-1,1) ) 237 fstrCPL%velo (kkk0 ) = -b1*fstrDYN%ACC(kkk1,1) - b2*fstrDYN%VEL(kkk1,1) + & 238 b3*( hecMAT%X(kkk1) - fstrDYN%DISP(kkk1,1) ) 239 fstrCPL%accel(kkk0-1) = -a1*fstrDYN%ACC(kkk1-1,1) - a2*fstrDYN%VEL(kkk1-1,1) + & 240 a3*( hecMAT%X(kkk1-1) - fstrDYN%DISP(kkk1-1,1) ) 241 fstrCPL%accel(kkk0 ) = -a1*fstrDYN%ACC(kkk1,1) - a2*fstrDYN%VEL(kkk1,1) + & 242 a3*( hecMAT%X(kkk1) - fstrDYN%DISP(kkk1,1) ) 243 endif 244 end do 245 call fstr_rcap_send( fstrCPL ) 246 endif 247 248 select case ( fstrPARAM%fg_couple_type ) 249 case (4) 250 call fstr_rcap_get( fstrCPL ) 251 case (5) 252 call fstr_get_convergence( revocap_flag ) 253 if( revocap_flag==0 ) then 254 do j = 1, hecMAT%NP * ndof 255 hecMAT%B(j) = prevB(j) 256 enddo 257 cycle 258 endif 259 case (6) 260 call fstr_get_convergence( revocap_flag ) 261 if( revocap_flag==0 ) then 262 do j = 1, hecMAT%NP * ndof 263 hecMAT%B(j) = prevB(j) 264 enddo 265 call fstr_rcap_get( fstrCPL ) 266 cycle 267 else 268 if( i /= fstrDYN%n_step ) call fstr_rcap_get( fstrCPL ) 269 endif 270 end select 271 endif 272 exit 273 enddo 274 275 !C ***************************************************** 276 !C-- contact corrector 277 !C 278 do j = 1 ,ndof*nnod 279 fstrSOLID%unode(j) = fstrDYN%DISP(j,1) 280 fstrSOLID%dunode(j) = hecMAT%X(j)-fstrDYN%DISP(j,1) 281 enddo 282 if( associated( fstrSOLID%contacts ) ) then 283 !call fstr_scan_contact_state( 1, fstrDYN%t_delta, kcaSLAGRANGE, hecMESH, fstrSOLID, infoCTChange ) 284 call forward_increment_Lagrange(1,ndof,fstrDYN%VEC1,hecMESH,fstrSOLID,infoCTChange,& 285 & fstrDYN%DISP(:,2),fstrSOLID%ddunode) 286 do j = 1 ,ndof*nnod 287 hecMAT%X(j) = hecMAT%X(j) + fstrSOLID%ddunode(j) 288 enddo 289 endif 290 291 !C-- new displacement, velocity and accelaration 292 do j = 1 ,ndof*nnod 293 fstrDYN%ACC (j,1) = a1*(hecMAT%X(j) - 2.d0*fstrDYN%DISP(j,1) + fstrDYN%DISP(j,3)) 294 fstrDYN%VEL (j,1) = a2*(hecMAT%X(j) - fstrDYN%DISP(j,3)) 295 fstrSOLID%unode(j) = fstrDYN%DISP(j,1) 296 fstrSOLID%dunode(j) = hecMAT%X(j)-fstrDYN%DISP(j,1) 297 fstrDYN%DISP(j,3) = fstrDYN%DISP(j,1) 298 fstrDYN%DISP(j,1) = hecMAT%X(j) 299 hecMAT%X(j) = fstrSOLID%dunode(j) 300 end do 301 302 ! ----- update strain, stress, and internal force 303 call fstr_UpdateNewton( hecMESH, hecMAT, fstrSOLID, fstrDYN%t_curr, fstrDYN%t_delta, 1 ) 304 305 do j = 1 ,ndof*nnod 306 fstrSOLID%unode(j) = fstrSOLID%unode(j) + fstrSOLID%dunode(j) 307 end do 308 call fstr_UpdateState( hecMESH, fstrSOLID, fstrDYN%t_delta ) 309 310 if( fstrDYN%restart_nout > 0 .and. & 311 (mod(i,fstrDYN%restart_nout).eq.0 .or. i.eq.fstrDYN%n_step) ) then 312 call fstr_write_restart_dyna_nl(i,hecMESH,fstrSOLID,fstrDYN,fstrPARAM) 313 end if 314 ! 315 !C-- output new displacement, velocity and accelaration 316 call fstr_dynamic_Output(hecMESH, fstrSOLID, fstrDYN, fstrPARAM) 317 call dynamic_output_monit(hecMESH, fstrPARAM, fstrDYN, fstrEIG, fstrSOLID) 318 319 enddo 320 321 if( fstrPARAM%fg_couple == 1) then 322 if( fstrPARAM%fg_couple_type==5 .or. & 323 fstrPARAM%fg_couple_type==6 ) then 324 deallocate( prevB ,stat=ierror ) 325 if( ierror /= 0 ) then 326 write(idbg,*) 'stop due to deallocation error <fstr_solve_NONLINEAR_DYNAMIC, prevB>' 327 write(idbg,*) ' rank = ', hecMESH%my_rank,' ierror = ',ierror 328 call flush(idbg) 329 call hecmw_abort( hecmw_comm_get_comm()) 330 endif 331 endif 332 endif 333 334 call hecmw_mpc_mat_finalize_explicit(hecMESH, hecMAT, hecMATmpc) 335 336 end subroutine fstr_solve_dynamic_nlexplicit 337 338 !< This subroutine implements Forward increment Lagrange multiplier method( NJ Carpenter et al. Int.J.Num.Meth.Eng.,32(1991),103-128 ) 339 subroutine forward_increment_Lagrange(cstep,ndof,mmat,hecMESH,fstrSOLID,infoCTChange,wkarray,uc) 340 integer, intent(in) :: cstep 341 integer, intent(in) :: ndof 342 real(kind=kreal), intent(in) :: mmat(:) 343 type( hecmwST_local_mesh ), intent(in) :: hecMESH !< type mesh 344 type(fstr_solid), intent(inout) :: fstrSOLID 345 type(fstr_info_contactChange) :: infoCTChange 346 real(kind=kreal), intent(out) :: wkarray(:) 347 real(kind=kreal), intent(out) :: uc(:) 348 integer :: i, j, k, m, grpid, slave, nn, iSS, sid, etype, iter 349 real(kind=kreal) :: fdum, conv, dlambda, shapefunc(l_max_surface_node), lambda(3) 350 351 call fstr_scan_contact_state_exp( cstep, hecMESH, fstrSOLID, infoCTChange ) 352 if( .not. infoCTChange%active ) return 353 354 uc = 0.0d0 355 356 iter = 0 357 do 358 wkarray = 0.0d0 359 do i=1,size(fstrSOLID%contacts) 360 do j= 1, size(fstrSOLID%contacts(i)%slave) 361 if( fstrSOLID%contacts(i)%states(j)%state == CONTACTFREE ) cycle 362 if( fstrSOLID%contacts(i)%states(j)%distance>epsilon(1.d0) ) then 363 fstrSOLID%contacts(i)%states(j)%state = CONTACTFREE 364 cycle 365 endif 366 if( iter==0 ) then 367 fstrSOLID%contacts(i)%states(j)%multiplier(:) =0.d0 368 fstrSOLID%contacts(i)%states(j)%wkdist =0.d0 369 cycle 370 endif 371 slave = fstrSOLID%contacts(i)%slave(j) 372 373 sid = fstrSOLID%contacts(i)%states(j)%surface 374 nn = size( fstrSOLID%contacts(i)%master(sid)%nodes ) 375 etype = fstrSOLID%contacts(i)%master(sid)%etype 376 call getShapeFunc( etype, fstrSOLID%contacts(i)%states(j)%lpos(:), shapefunc ) 377 wkarray( slave ) = -fstrSOLID%contacts(i)%states(j)%multiplier(1) 378 do k=1,nn 379 iSS = fstrSOLID%contacts(i)%master(sid)%nodes(k) 380 wkarray( iSS ) = wkarray( iSS ) + shapefunc(k) * fstrSOLID%contacts(i)%states(j)%multiplier(1) 381 enddo 382 enddo 383 enddo 384 385 if(iter > 0)then 386 do i=1,size(fstrSOLID%contacts) 387 do j= 1, size(fstrSOLID%contacts(i)%slave) 388 if( fstrSOLID%contacts(i)%states(j)%state == CONTACTFREE ) cycle 389 slave = fstrSOLID%contacts(i)%slave(j) 390 sid = fstrSOLID%contacts(i)%states(j)%surface 391 nn = size( fstrSOLID%contacts(i)%master(sid)%nodes ) 392 etype = fstrSOLID%contacts(i)%master(sid)%etype 393 call getShapeFunc( etype, fstrSOLID%contacts(i)%states(j)%lpos(:), shapefunc ) 394 fstrSOLID%contacts(i)%states(j)%wkdist = -wkarray( slave )/mmat( (slave-1)*ndof+1 ) 395 do k=1,nn 396 iSS = fstrSOLID%contacts(i)%master(sid)%nodes(k) 397 fstrSOLID%contacts(i)%states(j)%wkdist = fstrSOLID%contacts(i)%states(j)%wkdist & 398 + shapefunc(k) * wkarray(iSS) / mmat( (iSS-1)*ndof+1 ) 399 enddo 400 enddo 401 enddo 402 endif 403 404 conv = 0.d0 405 wkarray = 0.d0 406 do i=1,size(fstrSOLID%contacts) 407 do j= 1, size(fstrSOLID%contacts(i)%slave) 408 if( fstrSOLID%contacts(i)%states(j)%state == CONTACTFREE ) cycle 409 slave = fstrSOLID%contacts(i)%slave(j) 410 sid = fstrSOLID%contacts(i)%states(j)%surface 411 nn = size( fstrSOLID%contacts(i)%master(sid)%nodes ) 412 etype = fstrSOLID%contacts(i)%master(sid)%etype 413 call getShapeFunc( etype, fstrSOLID%contacts(i)%states(j)%lpos(:), shapefunc ) 414 fdum = 1.d0/mmat( (slave-1)*ndof+1 ) 415 do k=1,nn 416 iSS = fstrSOLID%contacts(i)%master(sid)%nodes(k) 417 fdum = fdum + shapefunc(k)*shapefunc(k)/mmat( (iSS-1)*ndof+1 ) 418 enddo 419 dlambda= (fstrSOLID%contacts(i)%states(j)%distance-fstrSOLID%contacts(i)%states(j)%wkdist) /fdum 420 conv = conv + dlambda*dlambda; 421 fstrSOLID%contacts(i)%states(j)%multiplier(1) = fstrSOLID%contacts(i)%states(j)%multiplier(1) + dlambda 422 if( fstrSOLID%contacts(i)%fcoeff>0.d0 ) then 423 if( fstrSOLID%contacts(i)%states(j)%state == CONTACTSLIP ) then 424 fstrSOLID%contacts(i)%states(j)%multiplier(2) = & 425 fstrSOLID%contacts(i)%fcoeff * fstrSOLID%contacts(i)%states(j)%multiplier(1) 426 else ! stick 427 ! fstrSOLID%contacts(i)%states(j)%multiplier(2) = 428 endif 429 endif 430 lambda = fstrSOLID%contacts(i)%states(j)%multiplier(1)* fstrSOLID%contacts(i)%states(j)%direction 431 wkarray((slave-1)*ndof+1:(slave-1)*ndof+3) = lambda(:) 432 do k=1,nn 433 iSS = fstrSOLID%contacts(i)%master(sid)%nodes(k) 434 wkarray((iSS-1)*ndof+1:(iSS-1)*ndof+3) = wkarray((iSS-1)*ndof+1:(iSS-1)*ndof+3) -lambda(:)*shapefunc(k) 435 enddo 436 enddo 437 enddo 438 if( dsqrt(conv)<1.d-8 ) exit 439 iter = iter+1 440 enddo 441 442 do i=1,hecMESH%n_node*ndof 443 uc(i) = wkarray(i)/mmat(i) 444 enddo 445 end subroutine forward_increment_Lagrange 446 447end module fstr_dynamic_nlexplicit 448