28 ,fstrDYN,fstrRESULT,fstrPARAM,infoCTChange &
29 ,fstrCPL, restrt_step_num )
31 type(hecmwst_local_mesh) :: hecMESH
32 type(hecmwst_matrix) :: hecMAT
35 type(hecmwst_result_data) :: fstrRESULT
39 type(fstr_info_contactchange) :: infoCTChange
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
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
59 call hecmw_mpc_mat_init_explicit(hecmesh, hecmat, hecmatmpc)
61 hecmat%NDOF=hecmesh%n_dof
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 )
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
75 call hecmw_abort( hecmw_comm_get_comm())
80 fstrsolid%dunode(:) =0.d0
82 a1 = 1.d0/fstrdyn%t_delta**2
83 a2 = 1.d0/(2.d0*fstrdyn%t_delta)
85 call setmass(fstrsolid,hecmesh,hecmat,fstreig)
86 call hecmw_mpc_trans_mass(hecmesh, hecmat, fstreig%mass)
88 allocate(mark(hecmat%NP * hecmat%NDOF))
89 call hecmw_mpc_mark_slave(hecmesh, hecmat, mark)
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
99 call hecmw_abort( hecmw_comm_get_comm())
106 if( restrt_step_num == 1 )
then
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
116 if(
associated( fstrsolid%contacts ) )
then
117 call initialize_contact_output_vectors(fstrsolid,hecmat)
119 & fstrdyn%DISP(:,2),fstrsolid%ddunode)
122 do i= restrt_step_num, fstrdyn%n_step
125 fstrdyn%t_curr = fstrdyn%t_delta * i
129 do j=1, hecmesh%n_node* hecmesh%n_dof
130 hecmat%B(j)=hecmat%B(j)-fstrsolid%QFORCE(j)
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)
144 if( fstrparam%fg_couple == 1 )
then
145 if( fstrparam%fg_couple_type==1 .or. &
146 fstrparam%fg_couple_type==3 .or. &
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
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 )
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
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 )
173 call hecmw_mpc_trans_rhs(hecmesh, hecmat, hecmatmpc)
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)
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)
198 call hecmw_abort( hecmw_comm_get_comm())
201 call hecmw_mpc_tback_sol(hecmesh, hecmat, hecmatmpc)
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
210 kkk1 = fstrcpl%coupled_node(j)*3
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 )
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) )
230 kkk1 = fstrcpl%coupled_node(j)*2
232 fstrcpl%disp (kkk0-1) = hecmat%X(kkk1-1)
233 fstrcpl%disp (kkk0 ) = hecmat%X(kkk1 )
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) )
248 select case ( fstrparam%fg_couple_type )
253 if( revocap_flag==0 )
then
254 do j = 1, hecmat%NP * ndof
255 hecmat%B(j) = prevb(j)
261 if( revocap_flag==0 )
then
262 do j = 1, hecmat%NP * ndof
263 hecmat%B(j) = prevb(j)
279 fstrsolid%unode(j) = fstrdyn%DISP(j,1)
280 fstrsolid%dunode(j) = hecmat%X(j)-fstrdyn%DISP(j,1)
282 if(
associated( fstrsolid%contacts ) )
then
285 & fstrdyn%DISP(:,2),fstrsolid%ddunode)
287 hecmat%X(j) = hecmat%X(j) + fstrsolid%ddunode(j)
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)
303 call fstr_updatenewton( hecmesh, hecmat, fstrsolid, fstrdyn%t_curr, fstrdyn%t_delta, 1 )
306 fstrsolid%unode(j) = fstrsolid%unode(j) + fstrsolid%dunode(j)
310 if( fstrdyn%restart_nout > 0 .and. &
311 (mod(i,fstrdyn%restart_nout).eq.0 .or. i.eq.fstrdyn%n_step) )
then
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
329 call hecmw_abort( hecmw_comm_get_comm())
334 call hecmw_mpc_mat_finalize_explicit(hecmesh, hecmat, hecmatmpc)
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
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)
351 call fstr_scan_contact_state_exp( cstep, hecmesh, fstrsolid, infoctchange )
352 if( .not. infoctchange%active )
return
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
367 fstrsolid%contacts(i)%states(j)%multiplier(:) =0.d0
368 fstrsolid%contacts(i)%states(j)%wkdist =0.d0
371 slave = fstrsolid%contacts(i)%slave(j)
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)
379 iss = fstrsolid%contacts(i)%master(sid)%nodes(k)
380 wkarray( iss ) = wkarray( iss ) + shapefunc(k) * fstrsolid%contacts(i)%states(j)%multiplier(1)
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 )
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 )
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 )
416 iss = fstrsolid%contacts(i)%master(sid)%nodes(k)
417 fdum = fdum + shapefunc(k)*shapefunc(k)/mmat( (iss-1)*ndof+1 )
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)
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(:)
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)
438 if( dsqrt(conv)<1.d-8 )
exit
442 do i=1,hecmesh%n_node*ndof
443 uc(i) = wkarray(i)/mmat(i)
This module contains subroutines for nonlinear explicit dynamic analysis.
subroutine fstr_solve_dynamic_nlexplicit(hecMESH, hecMAT, fstrSOLID, fstrEIG, fstrDYN, fstrRESULT, fstrPARAM, infoCTChange, fstrCPL, restrt_step_num)
subroutine forward_increment_lagrange(cstep, ndof, mmat, hecMESH, fstrSOLID, infoCTChange, wkarray, uc)
This module contains functions to set acceleration boundary condition in dynamic analysis.
subroutine dynamic_explicit_ass_ac(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, iter)
This module contains functions to set velocity boundary condition in dynamic analysis.
subroutine dynamic_explicit_ass_vl(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, iter)
This module contains functions to set displacement boundary condition in dynamic analysis.
subroutine dynamic_explicit_ass_bc(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, iter)
This subroutine setup disp boundary condition.
This module contains functions relates to coupling analysis.
subroutine dynamic_mat_ass_couple(hecMESH, hecMAT, fstrSOLID, fstrCPL)
This module contains function to set boundary condition of external load in dynamic analysis.
subroutine dynamic_mat_ass_load(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, fstrPARAM, iter)
This function sets boundary condition of external load.
This module provides functions to output result.
subroutine dynamic_output_monit(hecMESH, fstrPARAM, fstrDYNAMIC, fstrEIG, fstrSOLID)
subroutine fstr_dynamic_output(hecMESH, fstrSOLID, fstrDYNAMIC, fstrPARAM)
Output result.
Set up lumped mass matrix.
subroutine setmass(fstrSOLID, hecMESH, hecMAT, fstrEIG)
subroutine, public fstr_rcap_send(fstrCPL)
subroutine, public fstr_rcap_get(fstrCPL)
subroutine fstr_get_convergence(revocap_flag)
This module provides functions to read in and write out restart fiels.
subroutine fstr_write_restart_dyna_nl(cstep, hecMESH, fstrSOLID, fstrDYNAMIC, fstrPARAM, contactNode)
write out restart file for nonlinear dynamic analysis
This module provides function to calcualte to do updates.
subroutine fstr_updatestate(hecMESH, fstrSOLID, tincr)
Update elastiplastic status.
subroutine fstr_updatenewton(hecMESH, hecMAT, fstrSOLID, time, tincr, iter, strainEnergy)
変位/応力・ひずみ/内力のアップデート
This module defined coomon data and basic structures for analysis.
integer(kind=kint), parameter imsg
integer(kind=kint), parameter idbg
This modules just summarizes all modules used in static analysis.
Data for coupling analysis.
Data for DYNAMIC ANSLYSIS (fstrDYNAMIC)
Package of data used by Lanczos eigenvalue solver.
FSTR INNER CONTROL PARAMETERS (fstrPARAM)