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