10 private :: update_abort
33 type (hecmwST_matrix) :: hecMAT
34 type (hecmwST_local_mesh) :: hecMESH
35 type (fstr_solid) :: fstrSOLID
36 real(kind=kreal),
intent(in) :: time
37 real(kind=kreal),
intent(in) :: tincr
38 integer,
intent(in) :: iter
40 integer(kind=kint) :: nodLOCAL(fstrSOLID%max_ncon)
41 real(kind=kreal) :: ecoord(3, fstrsolid%max_ncon)
42 real(kind=kreal) :: thick, thick0(6)
43 integer(kind=kint) :: ndof, itype, is, iE, ic_type, nn, icel, iiS, i, j, nbase
45 real(kind=kreal) :: total_disp(6, fstrsolid%max_ncon), du(6, fstrsolid%max_ncon), ddu(6, fstrsolid%max_ncon)
46 real(kind=kreal) :: triad_tri(9, fstrsolid%max_ncon)
47 real(kind=kreal) :: triad_cur(9, fstrsolid%max_ncon)
48 real(kind=kreal) :: triad_ref(9, fstrsolid%max_ncon)
49 real(kind=kreal) :: shell_drill(fstrsolid%max_ncon)
50 real(kind=kreal) :: tt(fstrsolid%max_ncon), tt0(fstrsolid%max_ncon), ttn(fstrsolid%max_ncon)
51 real(kind=kreal) :: qf(fstrsolid%max_ncon*6), coords(3, 3)
52 integer :: isect, ihead, cdsys_ID
53 integer :: ndim, initt
55 real(kind=kreal),
optional :: strainenergy
56 real(kind=kreal) :: tmp, lambda(1), ddlambda(1)
57 real(kind=kreal) :: ddaux(3,3)
60 fstrsolid%QFORCE=0.0d0
88 do itype = 1, hecmesh%n_elem_type
89 is = hecmesh%elem_type_index(itype-1)+1
90 ie = hecmesh%elem_type_index(itype )
91 ic_type= hecmesh%elem_type_item(itype)
92 if (hecmw_is_etype_link(ic_type)) cycle
93 if (hecmw_is_etype_patch(ic_type)) cycle
109 iis = hecmesh%elem_node_index(icel-1)
110 nn = hecmesh%elem_node_index(icel)-iis
114 do j = 1,
size(fstrsolid%elements(icel)%gausses)
115 thick0(1:6) = fstrsolid%elements(icel)%gausses(j)%stress(1:6)
116 thick = thick + dsqrt(dot_product(thick0(1:6),thick0(1:6)))
118 if( thick < 1.d-10 )
then
119 do j = 1,
size(fstrsolid%elements(icel)%gausses)
120 if(
associated(fstrsolid%elements(icel)%gausses(j)%fstatus) ) &
121 & fstrsolid%elements(icel)%gausses(j)%fstatus = 0.d0
126 nodlocal(j) = hecmesh%elem_node_item (iis+j)
128 ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
131 ddu(i,j) = hecmat%X(ndof*nodlocal(j)+i-ndof)
132 du(i,j) = fstrsolid%dunode(ndof*nodlocal(j)+i-ndof)
133 total_disp(i,j) = fstrsolid%unode(ndof*nodlocal(j)+i-ndof)
136 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
137 if( iselastoplastic(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype) .or. &
138 fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton )
then
139 tt0(j)=fstrsolid%last_temp( nodlocal(j) )
142 if( hecmesh%hecmw_flag_initcon == 1 ) tt0(j) = hecmesh%node_init_val_item(nodlocal(j))
143 if( initt>0 ) tt0(j) =
g_initialcnd(initt)%realval(nodlocal(j))
145 ttn(j) = fstrsolid%last_temp( nodlocal(j) )
146 tt(j) = fstrsolid%temperature( nodlocal(j) )
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 cdsys_id = hecmesh%section%sect_orien_ID(isect)
154 if( cdsys_id > 0 )
call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords, icel)
156 if (ic_type == 741 .or. ic_type == 743 .or. ic_type == 731)
then
158 nbase = 9*(nodlocal(j)-1)
159 triad_tri(1:9,j) = 0.0d0
160 triad_cur(1:9,j) = 0.0d0
161 triad_ref(1:9,j) = 0.0d0
162 shell_drill(j) = 0.0d0
163 if(
associated(fstrsolid%shell_dtriad) ) triad_tri(1:9,j) = fstrsolid%shell_dtriad(nbase+1:nbase+9)
164 if(
associated(fstrsolid%shell_triad) ) triad_cur(1:9,j) = fstrsolid%shell_triad(nbase+1:nbase+9)
165 if(
associated(fstrsolid%shell_ref_triad) ) triad_ref(1:9,j) = fstrsolid%shell_ref_triad(nbase+1:nbase+9)
166 if(
associated(fstrsolid%shell_ddrill) ) shell_drill(j) = fstrsolid%shell_ddrill(nodlocal(j))
171 if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 )
then
172 call update_c2( ic_type,nn,ecoord(1:3,1:nn),fstrsolid%elements(icel)%gausses(:), &
173 thick,fstrsolid%elements(icel)%iset, &
174 total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof), &
175 tt(1:nn), tt0(1:nn), ttn(1:nn) )
177 else if( ic_type == 301 )
then
178 call update_c1( ic_type,nn,ecoord(:,1:nn), thick, total_disp(1:3,1:nn), du(1:3,1:nn), &
179 qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:) )
181 else if( ic_type == 361 )
then
182 if( fstrsolid%sections(isect)%elemopt361 ==
kel361fi )
then
183 call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
184 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
185 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361bbar )
then
186 call update_c3d8bbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
187 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
188 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361ic )
then
189 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,&
190 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, &
191 fstrsolid%elements(icel)%aux, ddaux(1:3,1:3), tt(1:nn), tt0(1:nn), ttn(1:nn) )
192 fstrsolid%elements(icel)%aux(1:3,1:3) = fstrsolid%elements(icel)%aux(1:3,1:3) + ddaux(1:3,1:3)
193 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361fbar )
then
194 call update_c3d8fbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
195 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
196 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361up )
then
197 lambda(1) = -0.5d0*fstrsolid%elements(icel)%p(1)
199 ( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), &
200 cdsys_id, coords, qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), &
201 iter, time, tincr, 1, lambda, ddlambda, tt(1:nn), tt0(1:nn) )
202 lambda(1) = lambda(1) + ddlambda(1)
203 fstrsolid%elements(icel)%p(1) = -2.0d0*lambda(1)
205 else if (ic_type == 341 .or. ic_type == 351 .or. ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 )
then
206 if( ic_type==341 .and. fstrsolid%sections(isect)%elemopt341 ==
kel341sesns ) cycle
207 call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
208 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
210 else if( ic_type == 511)
then
211 call update_connector( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), &
212 qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:) )
214 else if( ic_type == 611)
then
215 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
216 CALL updatest_beam(ic_type, nn, ecoord, total_disp(1:6,1:nn), du(1:6,1:nn), &
217 & hecmesh%section%sect_R_item(ihead+1:), fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof))
219 else if( ic_type == 641 )
then
220 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
221 call updatest_beam_641(ic_type, nn, ecoord, total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
222 & fstrsolid%elements(icel)%gausses(:), hecmesh%section%sect_R_item(ihead+1:), qf(1:nn*ndof))
224 else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) )
then
225 call update_shell_mitc(ic_type, nn, ndof, ecoord(1:3,1:nn), total_disp(1:ndof,1:nn), &
226 du(1:ndof,1:nn), fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 0, &
227 element=fstrsolid%elements(icel), ndtriad=triad_tri(1:9,1:nn), &
228 ndreftriad=triad_ref(1:9,1:nn), &
229 ndcurtriad=triad_cur(1:9,1:nn), nddrill=shell_drill(1:nn))
231 else if( ic_type == 761 )
then
232 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
233 call update_shell_mitc33(731, 3, 6, ecoord(1:3, 1:3), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
234 & fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 2)
236 else if( ic_type == 781 )
then
237 if(
fstrpr%nlgeom )
call update_abort( ic_type, 2 )
238 call update_shell_mitc33(741, 4, 6, ecoord(1:3, 1:4), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
239 & fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 1)
241 else if ( ic_type == 3414 )
then
242 if(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype /= incomp_newtonian) &
243 &
call update_abort( ic_type, 3, fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype )
245 ( ic_type, nn, ecoord(:,1:nn), total_disp(1:4,1:nn), du(1:4,1:nn), &
246 fstrsolid%elements(icel)%gausses(:) )
249 else if ( ic_type == 881 .or. ic_type == 891 )
then
250 call update_c3_sesns( ic_type, nn, nodlocal, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
251 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
254 write(*, *)
'###ERROR### : Element type not supported for nonlinear static analysis'
255 write(*, *)
' ic_type = ', ic_type
256 call hecmw_abort(hecmw_comm_get_comm())
262 call update_dummy( ndof, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), &
263 & du(1:3,1:nn), qf(1:nn*ndof), fstrsolid%elements(icel) )
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)
296 call hecmw_update_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node, ndof)
307 type(hecmwst_local_mesh) :: hecmesh
309 real(kind=kreal) :: tincr
310 integer(kind=kint) :: itype, is, iE, ic_type, icel, ngauss, i
311 integer(kind=kint) :: ishell
313 if(
associated( fstrsolid%temperature ) )
then
314 do i = 1, hecmesh%n_node
315 fstrsolid%last_temp(i) = fstrsolid%temperature(i)
319 do itype = 1, hecmesh%n_elem_type
320 is = hecmesh%elem_type_index(itype-1) + 1
321 ie = hecmesh%elem_type_index(itype )
322 ic_type= hecmesh%elem_type_item(itype)
323 if( ic_type == 301 ) ic_type = 111
324 if( hecmw_is_etype_link(ic_type) ) cycle
325 if( hecmw_is_etype_patch(ic_type) ) cycle
327 ngauss = numofquadpoints( ic_type )
329 if( iselastoplastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) )
then
333 elseif( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton )
then
334 if( tincr>0.d0 )
then
339 elseif( isviscoelastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) )
then
340 if( tincr > 0.d0 )
then
348 fstrsolid%elements(icel)%gausses(i)%strain_bak = fstrsolid%elements(icel)%gausses(i)%strain
349 fstrsolid%elements(icel)%gausses(i)%stress_bak = fstrsolid%elements(icel)%gausses(i)%stress
350 fstrsolid%elements(icel)%gausses(i)%strain_energy_bak = fstrsolid%elements(icel)%gausses(i)%strain_energy
353 if(
associated( fstrsolid%elements(icel)%shell_layer_gausses ) )
then
354 do ishell = 1,
size( fstrsolid%elements(icel)%shell_layer_gausses )
355 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_bak = &
356 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain
357 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%stress_bak = &
358 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%stress
359 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_energy_bak = &
360 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_energy
366 do i = 1, hecmesh%n_node
367 fstrsolid%QFORCE_bak(i) = fstrsolid%QFORCE(i)
372 subroutine update_abort( ic_type, flag, mtype )
373 integer(kind=kint),
intent(in) :: ic_type
374 integer(kind=kint),
intent(in) :: flag
375 integer(kind=kint),
intent(in),
optional :: mtype
378 write(*,*)
'###ERROR### : Element type not supported for static analysis'
379 else if( flag == 2 )
then
380 write(*,*)
'###ERROR### : Element type not supported for nonlinear static analysis'
381 else if( flag == 3 )
then
382 write(*,*)
'###ERROR### : This element is not supported for this material'
384 write(*,*)
' ic_type = ', ic_type
385 if(
present(mtype) )
write(*,*)
' mtype = ', mtype
386 call hecmw_abort(hecmw_comm_get_comm())
This module provide functions for elastoplastic calculation.
subroutine, public updateepstate(gauss)
Clear elatoplastic state.
This module defined elemact data and function.
subroutine update_dummy(ndof, nn, ecoord, u, du, qf, element)
integer, parameter kelact_inactive
Finite-rotation nodal kinematics for NLGEOM.
subroutine, public fstr_ensure_finite_rotation_state(hecMESH, fstrSOLID, ndof)
Build the per-node reference frames once, by averaging element shell triads at shared nodes....
This module provides function to calculate to do updates.
subroutine fstr_updatestate(hecMESH, fstrSOLID, tincr)
Update elastiplastic status.
subroutine fstr_updatenewton(hecMESH, hecMAT, fstrSOLID, time, tincr, iter, strainEnergy)
Update displacement, stress, strain and internal forces.
This module defines common data and basic structures for analysis.
integer(kind=kint), parameter kel361bbar
integer(kind=kint), parameter kel341sesns
integer(kind=kint), parameter kel361up
integer(kind=kint), parameter kel361fi
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords, icel)
This subroutine fetch coords defined by local coordinate system.
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.