FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
fstr_solve_NLGEOM.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
6
8 use m_fstr
17
18 implicit none
19
20contains
21
22 !======================================================================!
28 subroutine fstr_solve_nlgeom(hecMESH,hecMAT,fstrSOLID,fstrMAT,fstrPARAM,conMAT)
29 type (hecmwST_local_mesh) :: hecMESH
30 type (hecmwST_matrix ) :: hecMAT
31 type (fstr_param ) :: fstrPARAM
32 type (fstr_solid ) :: fstrSOLID
33 type (fstrST_matrix_contact_lagrange) :: fstrMAT
34 type (fstr_info_contactChange) :: infoCTChange, infoCTChange_bak
35 type (hecmwST_matrix ),optional :: conMAT
36
37 integer(kind=kint) :: ndof, nn
38 integer(kind=kint) :: j, i, tot_step, step_count, tot_step_print, CBbound
39 integer(kind=kint) :: sub_step
40 integer(kind=kint) :: restart_step_num, restart_substep_num
41 real(kind=kreal) :: ctime, dtime, endtime, factor
42 real(kind=kreal) :: time_1, time_2
43 logical :: ctchanged, is_OutPoint
44
45 if(hecmesh%my_rank==0) call fstr_timeinc_printstatus_init
46
47 hecmat%NDOF = hecmesh%n_dof
48
49 ndof = hecmat%NDOF
50 nn = ndof*ndof
51
52 if( fstrsolid%TEMP_ngrp_tot>0 .and. hecmesh%hecmw_flag_initcon==1 ) then
53 fstrsolid%last_temp = 0.0d0
54 fstrsolid%temperature = 0.0d0
55 do j=1, size(hecmesh%node_init_val_item)
56 i = hecmesh%node_init_val_index(j)
57 fstrsolid%last_temp(j) = hecmesh%node_init_val_item(i)
58 fstrsolid%temperature(j) = hecmesh%node_init_val_item(i)
59 end do
60 endif
61 if( fstrsolid%TEMP_ngrp_tot>0 .and. associated(g_initialcnd) ) then
62 fstrsolid%last_temp = 0.0d0
63 fstrsolid%temperature = 0.0d0
64 do j=1,size(g_initialcnd)
65 if( g_initialcnd(j)%cond_name=="temperature" ) then
66 if( .not. associated(fstrsolid%temperature) ) then
67 allocate( fstrsolid%temperature( hecmesh%n_node ) )
68 allocate( fstrsolid%temp_bak( hecmesh%n_node ) )
69 allocate( fstrsolid%last_temp( hecmesh%n_node ) )
70 endif
71 do i= 1, hecmesh%n_node
72 fstrsolid%last_temp(i) = g_initialcnd(j)%realval(i)
73 fstrsolid%temperature(i) = fstrsolid%last_temp(i)
74 enddo
75 endif
76 end do
77 endif
78
79 if( associated( fstrsolid%contacts ) ) then
80 call initialize_contact_output_vectors(fstrsolid,hecmat)
81 call setup_contact_elesurf_for_area( 1, hecmesh, fstrsolid )
82 endif
83
84 restart_step_num = 1
85 restart_substep_num = 1
86 fstrsolid%unode = 0.0d0
87 step_count = 0 !**
88 infoctchange%contactNode_previous = 0
89 infoctchange%contactNode_current = 0
90 if( fstrsolid%restart_nout < 0 ) then
91 call fstr_read_restart(restart_step_num,restart_substep_num,step_count,ctime,dtime,hecmesh,fstrsolid, &
92 fstrparam,infoctchange%contactNode_previous)
93 hecmat%Iarray(98) = 1
94 call fstr_set_time( ctime )
95 call fstr_set_timeinc_base( dtime )
96 fstrsolid%restart_nout = - fstrsolid%restart_nout
97 else
98 call fstr_static_output( 1, 0, 0.d0, hecmesh, fstrsolid, fstrparam, fstrpr%solution_type, .true. )
99 endif
100
101 fstrsolid%FACTOR = 0.0d0
102 call fstr_cutback_init( hecmesh, fstrsolid, fstrparam )
103 call fstr_cutback_save( fstrsolid, infoctchange, infoctchange_bak )
104
105 do tot_step=1, fstrsolid%nstep_tot
106 tot_step_print = tot_step+restart_step_num-1
107 if(hecmesh%my_rank==0) write(*,*) ''
108 if(hecmesh%my_rank==0) write(*,'(a,i5)') ' loading step=',tot_step_print
109
110 if( fstrsolid%TEMP_ngrp_tot>0 ) then
111 do j=1, hecmesh%n_node
112 fstrsolid%temp_bak(j) = fstrsolid%temperature(j)
113 end do
114 endif
115 call fstr_updatestate( hecmesh, fstrsolid, 0.0d0 )
116
117 ! -------------------------------------------------------------------------
118 ! STEP LOOP
119 ! -------------------------------------------------------------------------
120 sub_step = restart_substep_num
121 do while(.true.)
122
123 ! ----- time history of factor
124 call fstr_timeinc_settimeincrement( fstrsolid%step_ctrl(tot_step), fstrparam, sub_step, &
125 & fstrsolid%NRstat_i, fstrsolid%NRstat_r, fstrsolid%AutoINC_stat, fstrsolid%CutBack_stat )
126 if( fstrsolid%TEMP_irres > 0 ) then
127 fstrsolid%FACTOR(1) = 0.d0
128 fstrsolid%FACTOR(2) = 1.d0
129 call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time()+fstr_get_timeinc(), factor)
130 fstrsolid%TEMP_FACTOR = factor
131 else
132 call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time(),factor)
133 fstrsolid%FACTOR(1) = factor
134 call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time()+fstr_get_timeinc(), factor)
135 fstrsolid%FACTOR(2) = factor
136 endif
137
138 if(hecmesh%my_rank==0) then
139 write(*,'(A,I0,2(A,E12.4))') ' sub_step= ',sub_step,', &
140 & current_time=',fstr_get_time(), ', time_inc=',fstr_get_timeinc()
141 write(*,'(A,2f12.7)') ' loading_factor= ', fstrsolid%FACTOR
142 if( fstrsolid%TEMP_irres > 0 ) write(*,'(A,2f12.7)') ' readtemp_factor= ', fstrsolid%TEMP_FACTOR
143 endif
144
145 time_1 = hecmw_wtime()
146
147 ! analysis algorithm ( Newton-Rapshon Method )
148 if( .not. associated( fstrsolid%contacts ) ) then
149 call fstr_newton( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, &
150 restart_step_num, sub_step, fstr_get_time(), fstr_get_timeinc() )
151 else
152 if( fstrparam%contact_algo == kcaslagrange ) then
153 ! For Parallel Contact with Multi-Partition Domains
154 if(paracontactflag.and.present(conmat)) then
155 call fstr_newton_contactslag( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, fstrmat, &
156 restart_step_num, restart_substep_num, sub_step, fstr_get_time(), fstr_get_timeinc(), infoctchange, conmat )
157 else
158 call fstr_newton_contactslag( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, fstrmat, &
159 restart_step_num, restart_substep_num, sub_step, fstr_get_time(), fstr_get_timeinc(), infoctchange )
160 endif
161 else if( fstrparam%contact_algo == kcaalagrange ) then
162 call fstr_newton_contactalag( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, &
163 restart_step_num, restart_substep_num, sub_step, fstr_get_time(), fstr_get_timeinc(), infoctchange )
164 endif
165 endif
166
167 ! Time Increment
168 if( hecmesh%my_rank == 0 ) call fstr_timeinc_printstatus( fstrsolid%step_ctrl(tot_step), fstrparam, &
169 & tot_step_print, sub_step, fstrsolid%NRstat_i, fstrsolid%NRstat_r, &
170 & fstrsolid%AutoINC_stat, fstrsolid%CutBack_stat )
171 if( fstr_cutback_active() ) then
172
173 if( fstrsolid%CutBack_stat == 0 ) then ! converged
174 call fstr_cutback_save( fstrsolid, infoctchange, infoctchange_bak ) ! save analysis state
175 call fstr_proceed_time() ! current time += time increment
176
177 else ! not converged
178 cbbound = fstrparam%ainc(fstrsolid%step_ctrl(tot_step)%AincParam_id)%CBbound
179 if( fstrsolid%CutBack_stat == cbbound ) then
180 if( hecmesh%my_rank == 0 ) then
181 write(*,*) 'Number of successive cutback reached max number: ',cbbound
183 endif
184 call hecmw_abort( hecmw_comm_get_comm() )
185 endif
186 call fstr_cutback_load( fstrsolid, infoctchange, infoctchange_bak ) ! load analysis state
187 call fstr_set_contact_active( infoctchange%contactNode_current > 0 )
188
189 ! restore matrix structure for slagrange contact analysis
190 if( associated( fstrsolid%contacts ) .and. fstrparam%contact_algo == kcaslagrange ) then
191 if(paracontactflag.and.present(conmat)) then
192 call fstr_mat_con_contact( tot_step, hecmat, fstrsolid, fstrmat, infoctchange, conmat)
193 conmat%B(:) = 0.0d0
194 else
195 call fstr_mat_con_contact( tot_step, hecmat, fstrsolid, fstrmat, infoctchange)
196 endif
197 call solve_lineq_contact_init(hecmesh, hecmat, fstrmat, fstr_is_matrixstruct_symmetric(fstrsolid, hecmesh))
198 endif
199 if( hecmesh%my_rank == 0 ) write(*,*) '### State has been restored at time =',fstr_get_time()
200
201 !stop if # of substeps reached upper bound.
202 if( sub_step == fstrsolid%step_ctrl(tot_step)%num_substep ) then
203 if( hecmesh%my_rank == 0 ) then
204 write(*,'(a,i5,a,f6.3)') '### Number of substeps reached max number: at total_step=', &
205 & tot_step_print, ' time=', fstr_get_time()
206 endif
207 call hecmw_abort( hecmw_comm_get_comm())
208 endif
209
210 ! output time
211 time_2 = hecmw_wtime()
212 if( hecmesh%my_rank==0) write(imsg,'(a,",",2(I8,","),f10.2)') &
213 & 'step, substep, solve (sec) :', tot_step_print, sub_step, time_2 - time_1
214 cycle
215 endif
216 else
217 if( fstrsolid%CutBack_stat > 0 ) stop
218 call fstr_proceed_time() ! current time += time increment
219 endif
220
221 step_count = step_count + 1
222
223 ! ----- Restart
224 if( fstrsolid%restart_nout > 0 .and. mod(step_count,fstrsolid%restart_nout) == 0 ) then
225 call fstr_write_restart(tot_step,tot_step_print,sub_step,step_count,fstr_get_time(), &
226 & fstr_get_timeinc_base(), hecmesh,fstrsolid,fstrparam,.false.,infoctchange%contactNode_current)
227 endif
228
229 ! ----- Result output (include visualize output)
230 is_outpoint = fstr_timeinc_istimepoint( fstrsolid%step_ctrl(tot_step), fstrparam ) &
231 & .or. fstr_timeinc_isstepfinished( fstrsolid%step_ctrl(tot_step) )
232 call fstr_static_output( tot_step, step_count, fstr_get_time(), hecmesh, fstrsolid, fstrparam, &
233 & fstrpr%solution_type, is_outpoint )
234
235 time_2 = hecmw_wtime()
236 if( hecmesh%my_rank==0 ) then
237 write(imsg,'(A,",",2(I8,","),f10.2)') 'step, substep, solve (sec) :', tot_step_print, sub_step, time_2 - time_1
238 write(imsg,'(A,I0,",",1pE15.8)') '### stepcount (for output), time :', step_count, fstr_get_time()
239 endif
240
241 !if time reached the end time of step, exit loop.
242 if( fstr_timeinc_isstepfinished( fstrsolid%step_ctrl(tot_step) ) ) exit
243
244 if( sub_step == fstrsolid%step_ctrl(tot_step)%num_substep ) then
245 if( hecmesh%my_rank == 0 ) then
246 write(*,'(a,i5,a,f6.3)') '### Number of substeps reached max number: at total_step=', &
247 & tot_step_print, ' time=', fstr_get_time()
248 endif
249 if( hecmesh%my_rank == 0 ) call fstr_timeinc_printstatus_final(.false.)
250 stop !stop if # of substeps reached upper bound.
251 endif
252
253 sub_step = sub_step + 1
254 enddo !--- end of substep loop
255
256 ! ----- Restart at the end of step
257 if( fstrsolid%restart_nout > 0 ) then
258 call fstr_write_restart(tot_step,tot_step_print,sub_step,step_count,fstr_get_time(),fstr_get_timeinc_base(), &
259 & hecmesh,fstrsolid,fstrparam,.true.,infoctchange%contactNode_current)
260 endif
261 restart_substep_num = 1
262 if( fstrsolid%TEMP_irres > 0 ) exit
263 enddo !--- end of tot_step loop
264
265 call fstr_cutback_finalize( fstrsolid )
266
267 ! message
268 if(myrank == 0)then
270 write(imsg,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
271 write(*,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
272 endif
273
274 end subroutine fstr_solve_nlgeom
275
276 !C================================================================C
279 !C================================================================C
280 subroutine table_nlsta(hecMESH, fstrSOLID, cstep, time, f_t)
281 type ( hecmwST_local_mesh ), intent(in) :: hecMESH
282 type ( fstr_solid ), intent(in) :: fstrSOLID
283 integer(kind=kint), intent(in) :: cstep
284 real(kind=kreal), intent(in) :: time
285 real(kind=kreal), intent(out) :: f_t
286
287 integer(kind=kint) :: i
288 integer(kind=kint) :: jj_n_amp, jj1, jj2
289 integer(kind=kint) :: s1, s2, flag
290 real(kind=kreal) :: t_1, t_2, t_t, f_1, f_2, tincre
291
292 s1 = 0; s2 = 0
293 jj_n_amp = fstrsolid%step_ctrl( cstep )%amp_id
294
295 if( jj_n_amp <= 0 ) then ! Amplitude not defined
296 f_t = (time-fstrsolid%step_ctrl(cstep)%starttime)/fstrsolid%step_ctrl(cstep)%elapsetime
297 if( f_t>1.d0 ) f_t=1.d0
298 else
299 tincre = fstrsolid%step_ctrl( cstep )%initdt
300 jj1 = hecmesh%amp%amp_index(jj_n_amp - 1)
301 jj2 = hecmesh%amp%amp_index(jj_n_amp)
302
303 jj1 = jj1 + 2
304 t_t = time-fstrsolid%step_ctrl(cstep)%starttime
305
306 ! if(jj2 .eq. 0) then
307 ! f_t = 1.0
308 if(t_t .gt. hecmesh%amp%amp_table(jj2)) then
309 f_t = hecmesh%amp%amp_val(jj2)
310 else if(t_t .le. hecmesh%amp%amp_table(jj2)) then
311 flag=0
312 do i = jj1, jj2
313 if(t_t .le. hecmesh%amp%amp_table(i)) then
314 s2 = i
315 s1 = i - 1
316 flag = 1
317 endif
318 if( flag == 1 ) exit
319 end do
320
321 t_2 = hecmesh%amp%amp_table(s2)
322 t_1 = hecmesh%amp%amp_table(s1)
323 f_2 = hecmesh%amp%amp_val(s2)
324 f_1 = hecmesh%amp%amp_val(s1)
325 if( t_2-t_1 .lt. 1.0e-20) then
326 if(myrank == 0) then
327 write(imsg,*) 'stop due to t_2-t_1 <= 0'
328 endif
329 call hecmw_abort( hecmw_comm_get_comm())
330 endif
331 f_t = ((t_2*f_1 - t_1*f_2) + (f_2 - f_1)*t_t) / (t_2 - t_1)
332 endif
333
334 endif
335
336 end subroutine table_nlsta
337
338end module m_fstr_solve_nlgeom
This module provides functions of reconstructing.
logical function fstr_is_matrixstruct_symmetric(fstrsolid, hecmesh)
this function judges whether sitiffness matrix is symmetric or not
subroutine fstr_mat_con_contact(cstep, hecmat, fstrsolid, fstrmat, infoctchange, conmat)
this subroutine reconstructs node-based (stiffness) matrix structure \corresponding to contact state
This module provides functions to deal with cutback.
Definition: fstr_Cutback.f90:7
subroutine fstr_cutback_finalize(fstrsolid)
Finalizer of cutback variables.
subroutine fstr_cutback_init(hecmesh, fstrsolid, fstrparam)
Initializer of cutback variables.
subroutine fstr_cutback_load(fstrsolid, infoctchange, infoctchange_bak)
Load analysis status.
subroutine fstr_cutback_save(fstrsolid, infoctchange, infoctchange_bak)
Save analysis status.
logical function fstr_cutback_active()
This module provides functions on nonlinear analysis.
subroutine fstr_newton_contactslag(cstep, hecmesh, hecmat, fstrsolid, fstrparam, fstrmat, restart_step_num, restart_substep_num, sub_step, ctime, dtime, infoctchange, conmat)
This subroutine solve nonlinear solid mechanics problems by Newton-Raphson method....
subroutine fstr_newton(cstep, hecmesh, hecmat, fstrsolid, fstrparam, restrt_step_num, sub_step, ctime, dtime)
This subroutine solve nonlinear solid mechanics problems by Newton-Raphson method.
subroutine fstr_newton_contactalag(cstep, hecmesh, hecmat, fstrsolid, fstrparam, restart_step_num, restart_substep_num, sub_step, ctime, dtime, infoctchange)
This subroutine solve nonlinear solid mechanics problems by Newton-Raphson method combined with Neste...
This module provides functions to read in and write out restart fiels.
Definition: fstr_Restart.f90:8
subroutine fstr_write_restart(cstep, cstep_ext, substep, step_count, ctime, dtime, hecmesh, fstrsolid, fstrparam, is_stepfinished, contactnode)
write out restart file
subroutine fstr_read_restart(cstep, substep, step_count, ctime, dtime, hecmesh, fstrsolid, fstrparam, contactnode)
Read in restart file.
This module provides main suboruitne for nonliear calculation.
subroutine fstr_solve_nlgeom(hecmesh, hecmat, fstrsolid, fstrmat, fstrparam, conmat)
This module provides main suborutine for nonlinear calculation.
subroutine table_nlsta(hecmesh, fstrsolid, cstep, time, f_t)
This subroutine decide the loading increment considering the amplitude definition.
This module provides functions to deal with time and increment of stress analysis.
logical function fstr_timeinc_istimepoint(stepinfo, fstrparam)
real(kind=kreal) function fstr_get_timeinc()
real(kind=kreal) function fstr_get_timeinc_base()
subroutine fstr_set_timeinc_base(dtime_base)
real(kind=kreal) function fstr_get_time()
subroutine fstr_timeinc_printstatus(stepinfo, fstrparam, totstep, substep, nrstati, nrstatr, autoinc_stat, cutback_stat)
subroutine fstr_timeinc_settimeincrement(stepinfo, fstrparam, substep, nrstati, nrstatr, autoinc_stat, cutback_stat)
subroutine fstr_proceed_time()
subroutine fstr_timeinc_printstatus_final(success_flag)
subroutine fstr_timeinc_printstatus_init
subroutine fstr_set_time(time)
logical function fstr_timeinc_isstepfinished(stepinfo)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:80
integer(kind=kint), parameter imsg
Definition: m_fstr.f90:94
integer(kind=kint), parameter kcaslagrange
contact analysis algorithm
Definition: m_fstr.f90:53
integer(kind=kint), parameter kcaalagrange
Definition: m_fstr.f90:54
type(fstr_param), target fstrpr
GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP.
Definition: m_fstr.f90:190
type(tinitialcondition), dimension(:), pointer, save g_initialcnd
Definition: m_fstr.f90:135
logical paracontactflag
PARALLEL CONTACT FLAG.
Definition: m_fstr.f90:84
This module provides functions to solve sparse system of \linear equitions in the case of contact ana...
subroutine, public solve_lineq_contact_init(hecmesh, hecmat, fstrmat, is_sym)
This subroutine.
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
This module provides functions to output result.
subroutine fstr_static_output(cstep, istep, time, hecmesh, fstrsolid, fstrparam, flag, outflag)
Output result.