10 private :: update_abort
29 type (hecmwST_matrix) :: hecMAT
30 type (hecmwST_local_mesh) :: hecMESH
31 type (fstr_solid) :: fstrSOLID
32 real(kind=kreal),
intent(in) :: time
33 real(kind=kreal),
intent(in) :: tincr
34 integer,
intent(in) :: iter
36 integer(kind=kint) :: nodLOCAL(20)
37 real(kind=kreal) :: ecoord(3, 20), stiff(60, 60)
38 real(kind=kreal) :: thick
39 integer(kind=kint) :: ndof, itype, is, iE, ic_type, nn, icel, iiS, i, j
41 real(kind=kreal) :: total_disp(6, 20), du(6, 20), ddu(6, 20)
42 real(kind=kreal) :: tt(20), tt0(20), ttn(20), qf(20*6), coords(3, 3)
43 integer :: ig0, ig, ik, in, ierror, isect, ihead, cdsys_ID
44 integer :: ndim, initt
46 real(kind=kreal),
optional :: strainenergy
47 real(kind=kreal) :: tmp
48 real(kind=kreal) :: ddaux(3,3)
51 fstrsolid%QFORCE=0.0d0
78 do itype = 1, hecmesh%n_elem_type
79 is = hecmesh%elem_type_index(itype-1)+1
80 ie = hecmesh%elem_type_index(itype )
81 ic_type= hecmesh%elem_type_item(itype)
82 if (hecmw_is_etype_link(ic_type)) cycle
83 if (hecmw_is_etype_patch(ic_type)) cycle
84 nn = hecmw_get_max_node(ic_type)
85 if( nn>20 ) stop
"Elemental nodes>20!"
98 iis = hecmesh%elem_node_index(icel-1)
101 nodlocal(j) = hecmesh%elem_node_item (iis+j)
103 ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
105 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
106 if( iselastoplastic(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype) .or. &
107 fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton )
then
108 tt0(j)=fstrsolid%last_temp( nodlocal(j) )
111 if( hecmesh%hecmw_flag_initcon == 1 ) tt0(j) = hecmesh%node_init_val_item(nodlocal(j))
112 if( initt>0 ) tt0(j) =
g_initialcnd(initt)%realval(nodlocal(j))
114 ttn(j) = fstrsolid%last_temp( nodlocal(j) )
115 tt(j) = fstrsolid%temperature( nodlocal(j) )
121 nodlocal(j) = hecmesh%elem_node_item (iis+j)
123 ddu(i,j) = hecmat%X(ndof*nodlocal(j)+i-ndof)
124 du(i,j) = fstrsolid%dunode(ndof*nodlocal(j)+i-ndof)
125 total_disp(i,j) = fstrsolid%unode(ndof*nodlocal(j)+i-ndof)
129 isect = hecmesh%section_ID(icel)
130 cdsys_id = hecmesh%section%sect_orien_ID(isect)
131 if( cdsys_id > 0 )
call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
134 if( getspacedimension( ic_type ) == 2 ) thick = 1.0d0
136 if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 )
then
138 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
139 call update_c2( ic_type,nn,ecoord(1:3,1:nn),fstrsolid%elements(icel)%gausses(:), &
140 thick,fstrsolid%elements(icel)%iset, &
141 total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof), &
142 tt(1:nn), tt0(1:nn), ttn(1:nn) )
144 call update_c2( ic_type,nn,ecoord(1:3,1:nn),fstrsolid%elements(icel)%gausses(:), &
145 thick,fstrsolid%elements(icel)%iset, &
146 total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof) )
149 else if( ic_type == 301 )
then
150 isect= hecmesh%section_ID(icel)
151 ihead = hecmesh%section%sect_R_index(isect-1)
152 thick = hecmesh%section%sect_R_item(ihead+1)
153 call update_c1( ic_type,nn,ecoord(:,1:nn), thick, total_disp(1:3,1:nn), du(1:3,1:nn), &
154 qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:) )
156 else if( ic_type == 361 )
then
158 if( fstrsolid%sections(isect)%elemopt361 ==
kel361fi )
then
159 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
160 call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
161 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
163 call update_c3( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
164 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr )
166 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361bbar )
then
167 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
168 call update_c3d8bbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
169 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
171 call update_c3d8bbar( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
172 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr )
174 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361ic )
then
175 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
176 call update_c3d8ic( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), cdsys_id, coords,&
177 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, &
178 fstrsolid%elements(icel)%aux, ddaux(1:3,1:3), tt(1:nn), tt0(1:nn), ttn(1:nn) )
180 call update_c3d8ic( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), cdsys_id, coords,&
181 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, &
182 fstrsolid%elements(icel)%aux, ddaux(1:3,1:3) )
184 fstrsolid%elements(icel)%aux(1:3,1:3) = fstrsolid%elements(icel)%aux(1:3,1:3) + ddaux(1:3,1:3)
185 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361fbar )
then
186 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
187 call update_c3d8fbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
188 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
190 call update_c3d8fbar( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
191 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr )
195 else if (ic_type == 341 .or. ic_type == 351 .or. ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 )
then
196 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
197 call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
198 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
200 call update_c3( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
201 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr )
204 else if( ic_type == 611)
then
205 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
206 isect = hecmesh%section_ID(icel)
207 ihead = hecmesh%section%sect_R_index(isect-1)
208 CALL updatest_beam(ic_type, nn, ecoord, total_disp(1:6,1:nn), du(1:6,1:nn), &
209 & hecmesh%section%sect_R_item(ihead+1:), fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof))
211 else if( ic_type == 641 )
then
212 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
213 isect = hecmesh%section_ID(icel)
214 ihead = hecmesh%section%sect_R_index(isect-1)
215 call updatest_beam_641(ic_type, nn, ecoord, total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
216 & fstrsolid%elements(icel)%gausses(:), hecmesh%section%sect_R_item(ihead+1:), qf(1:nn*ndof))
218 else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) )
then
219 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
220 isect = hecmesh%section_ID(icel)
221 ihead = hecmesh%section%sect_R_index(isect-1)
222 thick = hecmesh%section%sect_R_item(ihead+1)
223 call updatest_shell_mitc(ic_type, nn, ndof, ecoord(1:3, 1:nn), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
224 & fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 0)
226 else if( ic_type == 761 )
then
227 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
228 isect = hecmesh%section_ID(icel)
229 ihead = hecmesh%section%sect_R_index(isect-1)
230 thick = hecmesh%section%sect_R_item(ihead+1)
231 call updatest_shell_mitc33(731, 3, 6, ecoord(1:3, 1:3), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
232 & fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 2)
234 else if( ic_type == 781 )
then
235 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
236 isect = hecmesh%section_ID(icel)
237 ihead = hecmesh%section%sect_R_index(isect-1)
238 thick = hecmesh%section%sect_R_item(ihead+1)
239 call updatest_shell_mitc33(741, 4, 6, ecoord(1:3, 1:4), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
240 & fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 1)
242 else if ( ic_type == 3414 )
then
243 if(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype /= incomp_newtonian)
then
244 write(*, *)
'###ERROR### : This element is not supported for this material'
245 write(*, *)
'ic_type = ', ic_type,
', mtype = ', fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype
246 call hecmw_abort(hecmw_comm_get_comm())
249 ( ic_type, nn, ecoord(:,1:nn), total_disp(1:4,1:nn), du(1:4,1:nn), &
250 fstrsolid%elements(icel)%gausses(:) )
261 write(*, *)
'###ERROR### : Element type not supported for nonlinear static analysis'
262 write(*, *)
' ic_type = ', ic_type
263 call hecmw_abort(hecmw_comm_get_comm())
271 fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i) = fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i)+qf(ndof*(j-1)+i)
276 if(
present(strainenergy))
then
277 ndim = getspacedimension( fstrsolid%elements(icel)%etype )
280 tmp = 0.5d0*( fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i)+qf(ndim*(j-1)+i) )*ddu(i,j)
282 strainenergy = strainenergy+tmp
283 fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i) = qf(ndim*(j-1)+i)
297 call hecmw_update_3_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node)
298 else if( ndof==2 )
then
299 call hecmw_update_2_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node)
300 else if( ndof==4 )
then
301 call hecmw_update_4_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node)
302 else if( ndof==6 )
then
303 call hecmw_update_m_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node,6)
316 type(hecmwst_local_mesh) :: hecmesh
318 real(kind=kreal) :: tincr
319 integer(kind=kint) :: itype, is, iE, ic_type, icel, ngauss, i
321 if(
associated( fstrsolid%temperature ) )
then
322 do i = 1, hecmesh%n_node
323 fstrsolid%last_temp(i) = fstrsolid%temperature(i)
327 do itype = 1, hecmesh%n_elem_type
328 is = hecmesh%elem_type_index(itype-1) + 1
329 ie = hecmesh%elem_type_index(itype )
330 ic_type= hecmesh%elem_type_item(itype)
331 if( ic_type == 301 ) ic_type = 111
332 if( hecmw_is_etype_link(ic_type) ) cycle
333 if( hecmw_is_etype_patch(ic_type) ) cycle
335 ngauss = numofquadpoints( ic_type )
337 if( iselastoplastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) )
then
341 elseif( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton )
then
342 if( tincr>0.d0 )
then
347 elseif( isviscoelastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) )
then
348 if( tincr > 0.d0 )
then
356 fstrsolid%elements(icel)%gausses(i)%strain_bak = fstrsolid%elements(icel)%gausses(i)%strain
357 fstrsolid%elements(icel)%gausses(i)%stress_bak = fstrsolid%elements(icel)%gausses(i)%stress
363 subroutine update_abort( ic_type, flag )
364 integer(kind=kint),
intent(in) :: ic_type
365 integer(kind=kint),
intent(in) :: flag
368 write(*,*)
'###ERROR### : Element type not supported for static analysis'
369 else if( flag == 2 )
then
370 write(*,*)
'###ERROR### : Element type not supported for nonlinear static analysis'
372 write(*,*)
' ic_type = ', ic_type
373 call hecmw_abort(hecmw_comm_get_comm())
This module provide functions for elastoplastic calculation.
subroutine updateepstate(gauss)
Clear elatoplastic state.
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 kel361bbar
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords)
This subroutine fetch coords defined by local coordinate system.
integer(kind=kint), parameter kel361fi
section control
integer(kind=kint), parameter kel361ic
type(fstr_param), target fstrpr
GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP.
integer(kind=kint), parameter kel361fbar
type(tinitialcondition), dimension(:), pointer, save g_initialcnd
This modules just summarizes all modules used in static analysis.
This module provides functions for creep calculation.
subroutine updateviscostate(gauss)
Update viscoplastic state.
This module provides functions for viscoelastic calculation.
subroutine updateviscoelasticstate(gauss)
Update viscoplastic state.