10 private :: update_abort
27 subroutine fstr_updatenewton ( hecMESH, hecMAT, fstrSOLID, fstrPARAM, time, tincr,iter, strainEnergy)
33 type (hecmwST_matrix) :: hecMAT
34 type (hecmwST_local_mesh) :: hecMESH
35 type (fstr_solid) :: fstrSOLID
36 type (fstr_param) :: fstrPARAM
37 real(kind=kreal),
intent(in) :: time
38 real(kind=kreal),
intent(in) :: tincr
39 integer,
intent(in) :: iter
41 integer(kind=kint) :: nodLOCAL(fstrSOLID%max_ncon)
42 real(kind=kreal) :: ecoord(3, fstrsolid%max_ncon)
43 real(kind=kreal) :: thick, thick0(6)
44 integer(kind=kint) :: ndof, itype, is, iE, ic_type, nn, icel, iiS, i, j, nbase
46 real(kind=kreal) :: total_disp(6, fstrsolid%max_ncon), du(6, fstrsolid%max_ncon), ddu(6, fstrsolid%max_ncon)
47 real(kind=kreal) :: triad_tri(9, fstrsolid%max_ncon)
48 real(kind=kreal) :: triad_cur(9, fstrsolid%max_ncon)
49 real(kind=kreal) :: triad_ref(9, fstrsolid%max_ncon)
50 real(kind=kreal) :: shell_drill(fstrsolid%max_ncon)
51 real(kind=kreal) :: tt(fstrsolid%max_ncon), tt0(fstrsolid%max_ncon), ttn(fstrsolid%max_ncon)
52 real(kind=kreal) :: qf(fstrsolid%max_ncon*6), coords(3, 3)
53 integer :: isect, ihead, cdsys_ID
54 integer :: ndim, initt
56 real(kind=kreal),
optional :: strainenergy
57 real(kind=kreal) :: tmp, lambda(1), ddlambda(1)
58 real(kind=kreal) :: ddaux(3,3)
61 fstrsolid%QFORCE=0.0d0
89 do itype = 1, hecmesh%n_elem_type
90 is = hecmesh%elem_type_index(itype-1)+1
91 ie = hecmesh%elem_type_index(itype )
92 ic_type= hecmesh%elem_type_item(itype)
93 if (hecmw_is_etype_link(ic_type)) cycle
94 if (hecmw_is_etype_patch(ic_type)) cycle
110 iis = hecmesh%elem_node_index(icel-1)
111 nn = hecmesh%elem_node_index(icel)-iis
115 do j = 1,
size(fstrsolid%elements(icel)%gausses)
116 thick0(1:6) = fstrsolid%elements(icel)%gausses(j)%stress(1:6)
117 thick = thick + dsqrt(dot_product(thick0(1:6),thick0(1:6)))
119 if( thick < 1.d-10 )
then
120 do j = 1,
size(fstrsolid%elements(icel)%gausses)
121 if(
associated(fstrsolid%elements(icel)%gausses(j)%fstatus) ) &
122 & fstrsolid%elements(icel)%gausses(j)%fstatus = 0.d0
127 nodlocal(j) = hecmesh%elem_node_item (iis+j)
129 ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
132 ddu(i,j) = hecmat%X(ndof*nodlocal(j)+i-ndof)
133 du(i,j) = fstrsolid%dunode(ndof*nodlocal(j)+i-ndof)
134 total_disp(i,j) = fstrsolid%unode(ndof*nodlocal(j)+i-ndof)
137 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
138 if( iselastoplastic(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype) .or. &
139 fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton .or. &
140 fstrsolid%elements(icel)%gausses(1)%pMaterial%nlgeom_flag == updatelag )
then
141 tt0(j)=fstrsolid%last_temp( nodlocal(j) )
144 if( hecmesh%hecmw_flag_initcon == 1 ) tt0(j) = hecmesh%node_init_val_item(nodlocal(j))
145 if( initt>0 ) tt0(j) =
g_initialcnd(initt)%realval(nodlocal(j))
147 ttn(j) = fstrsolid%last_temp( nodlocal(j) )
148 tt(j) = fstrsolid%temperature( nodlocal(j) )
152 isect = hecmesh%section_ID(icel)
153 ihead = hecmesh%section%sect_R_index(isect-1)
154 thick = hecmesh%section%sect_R_item(ihead+1)
155 cdsys_id = hecmesh%section%sect_orien_ID(isect)
156 if( cdsys_id > 0 )
call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords, icel)
158 if (ic_type == 741 .or. ic_type == 743 .or. ic_type == 731)
then
160 nbase = 9*(nodlocal(j)-1)
161 triad_tri(1:9,j) = 0.0d0
162 triad_cur(1:9,j) = 0.0d0
163 triad_ref(1:9,j) = 0.0d0
164 shell_drill(j) = 0.0d0
165 if(
associated(fstrsolid%shell_dtriad) ) triad_tri(1:9,j) = fstrsolid%shell_dtriad(nbase+1:nbase+9)
166 if(
associated(fstrsolid%shell_triad) ) triad_cur(1:9,j) = fstrsolid%shell_triad(nbase+1:nbase+9)
167 if(
associated(fstrsolid%shell_ref_triad) ) triad_ref(1:9,j) = fstrsolid%shell_ref_triad(nbase+1:nbase+9)
168 if(
associated(fstrsolid%shell_ddrill) ) shell_drill(j) = fstrsolid%shell_ddrill(nodlocal(j))
172 if( ic_type == 761 .or. ic_type == 781 )
then
174 nbase = 9*(nodlocal(j+nn/2)-1)
175 triad_tri(1:9,j) = 0.0d0
176 triad_cur(1:9,j) = 0.0d0
177 triad_ref(1:9,j) = 0.0d0
178 shell_drill(j) = 0.0d0
179 if(
associated(fstrsolid%shell_dtriad) ) triad_tri(1:9,j) = fstrsolid%shell_dtriad(nbase+1:nbase+9)
180 if(
associated(fstrsolid%shell_triad) ) triad_cur(1:9,j) = fstrsolid%shell_triad(nbase+1:nbase+9)
181 if(
associated(fstrsolid%shell_ref_triad) ) triad_ref(1:9,j) = fstrsolid%shell_ref_triad(nbase+1:nbase+9)
182 if(
associated(fstrsolid%shell_ddrill) ) shell_drill(j) = fstrsolid%shell_ddrill(nodlocal(j+nn/2))
187 if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 )
then
188 call update_c2( ic_type,nn,ecoord(1:3,1:nn),fstrsolid%elements(icel)%gausses(:), &
189 thick,fstrsolid%elements(icel)%iset, &
190 total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof), &
191 tt(1:nn), tt0(1:nn), ttn(1:nn) )
193 else if( ic_type == 301 )
then
194 call update_c1( ic_type,nn,ecoord(:,1:nn), thick, total_disp(1:3,1:nn), du(1:3,1:nn), &
195 qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:) )
197 else if( ic_type == 361 )
then
198 if( fstrsolid%sections(isect)%elemopt361 ==
kel361fi )
then
199 call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
200 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
201 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361bbar )
then
202 call update_c3d8bbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
203 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
204 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361ic )
then
205 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,&
206 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, &
207 fstrsolid%elements(icel)%aux, ddaux(1:3,1:3), tt(1:nn), tt0(1:nn), ttn(1:nn) )
208 fstrsolid%elements(icel)%aux(1:3,1:3) = fstrsolid%elements(icel)%aux(1:3,1:3) + ddaux(1:3,1:3)
209 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361fbar )
then
210 call update_c3d8fbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
211 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
212 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361up )
then
213 lambda(1) = -0.5d0*fstrsolid%elements(icel)%p(1)
215 ( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), &
216 cdsys_id, coords, qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), &
217 iter, time, tincr, 1, lambda, ddlambda, tt(1:nn), tt0(1:nn) )
218 lambda(1) = lambda(1) + ddlambda(1)
219 fstrsolid%elements(icel)%p(1) = -2.0d0*lambda(1)
221 else if (ic_type == 341 .or. ic_type == 351 .or. ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 )
then
222 if( ic_type==341 .and. fstrsolid%sections(isect)%elemopt341 ==
kel341sesns ) cycle
223 call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
224 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
226 else if( ic_type == 511)
then
227 call update_connector( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), &
228 qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:) )
230 else if( ic_type == 611)
then
231 if( fstrparam%nlgeom )
call update_abort( ic_type, 2 )
232 CALL updatest_beam(ic_type, nn, ecoord, total_disp(1:6,1:nn), du(1:6,1:nn), &
233 hecmesh%section%sect_R_item(ihead+1:), fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), &
234 fstrsolid%sections(isect)%elemopt611)
236 else if( ic_type == 641 )
then
237 if( fstrparam%nlgeom )
call update_abort( ic_type, 2 )
238 call updatest_beam_641_from_611(ecoord(1:3,1:4), total_disp(1:ndof,1:4), du(1:ndof,1:4), &
239 fstrsolid%elements(icel)%gausses(:), hecmesh%section%sect_R_item(ihead+1:), &
240 qf(1:nn*ndof), fstrsolid%sections(isect)%elemopt611)
242 else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) )
then
243 call update_shell_mitc(ic_type, nn, ndof, ecoord(1:3,1:nn), total_disp(1:ndof,1:nn), &
244 du(1:ndof,1:nn), fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 0, &
245 element=fstrsolid%elements(icel), ndtriad=triad_tri(1:9,1:nn), &
246 ndreftriad=triad_ref(1:9,1:nn), &
247 ndcurtriad=triad_cur(1:9,1:nn), nddrill=shell_drill(1:nn), &
248 cdsys_id=cdsys_id, coords=coords)
250 else if( ic_type == 761 )
then
251 if( fstrparam%nlgeom )
then
252 call update_shell_mitc33(731, 3, 6, ecoord(1:3, 1:3), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
253 fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 2, &
254 element=fstrsolid%elements(icel), ndtriad=triad_tri(1:9,1:3), &
255 ndreftriad=triad_ref(1:9,1:3), ndcurtriad=triad_cur(1:9,1:3), nddrill=shell_drill(1:3), &
256 cdsys_id=cdsys_id, coords=coords)
258 call update_shell_mitc33(731, 3, 6, ecoord(1:3, 1:3), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
259 fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 2, &
260 cdsys_id=cdsys_id, coords=coords)
263 else if( ic_type == 781 )
then
264 call update_shell_mitc33(741, 4, 6, ecoord(1:3, 1:4), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
265 fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 1, &
266 element=fstrsolid%elements(icel), ndtriad=triad_tri(1:9,1:4), &
267 ndreftriad=triad_ref(1:9,1:4), ndcurtriad=triad_cur(1:9,1:4), nddrill=shell_drill(1:4), &
268 cdsys_id=cdsys_id, coords=coords)
270 else if ( ic_type == 3414 )
then
271 if(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype /= incomp_newtonian) &
272 &
call update_abort( ic_type, 3, fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype )
274 ( ic_type, nn, ecoord(:,1:nn), total_disp(1:4,1:nn), du(1:4,1:nn), &
275 fstrsolid%elements(icel)%gausses(:) )
278 else if ( ic_type == 881 .or. ic_type == 891 )
then
279 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, &
280 qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
284 write(*, *)
'###ERROR### : Element type not supported for nonlinear static analysis'
285 write(*, *)
' ic_type = ', ic_type
287 call hecmw_abort(hecmw_comm_get_comm())
293 call update_dummy( ndof, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), &
294 & du(1:3,1:nn), qf(1:nn*ndof), fstrsolid%elements(icel) )
302 fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i) = fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i)+qf(ndof*(j-1)+i)
307 if(
present(strainenergy))
then
308 ndim = getspacedimension( fstrsolid%elements(icel)%etype )
311 tmp = 0.5d0*( fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i)+qf(ndim*(j-1)+i) )*ddu(i,j)
313 strainenergy = strainenergy+tmp
314 fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i) = qf(ndim*(j-1)+i)
327 call hecmw_update_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node, ndof)
338 type(hecmwst_local_mesh) :: hecmesh
340 real(kind=kreal) :: tincr
341 integer(kind=kint) :: itype, is, iE, ic_type, icel, ngauss, i
342 integer(kind=kint) :: ishell
344 if(
associated( fstrsolid%temperature ) )
then
345 do i = 1, hecmesh%n_node
346 fstrsolid%last_temp(i) = fstrsolid%temperature(i)
350 do itype = 1, hecmesh%n_elem_type
351 is = hecmesh%elem_type_index(itype-1) + 1
352 ie = hecmesh%elem_type_index(itype )
353 ic_type= hecmesh%elem_type_item(itype)
354 if( ic_type == 301 ) ic_type = 111
355 if( hecmw_is_etype_link(ic_type) ) cycle
356 if( hecmw_is_etype_patch(ic_type) ) cycle
358 ngauss = numofquadpoints( ic_type )
360 if( iselastoplastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) )
then
364 elseif( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton )
then
365 if( tincr>0.d0 )
then
370 elseif( isviscoelastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) )
then
371 if( tincr > 0.d0 )
then
379 fstrsolid%elements(icel)%gausses(i)%strain_bak = fstrsolid%elements(icel)%gausses(i)%strain
380 fstrsolid%elements(icel)%gausses(i)%stress_bak = fstrsolid%elements(icel)%gausses(i)%stress
381 fstrsolid%elements(icel)%gausses(i)%strain_energy_bak = fstrsolid%elements(icel)%gausses(i)%strain_energy
384 if(
associated( fstrsolid%elements(icel)%shell_layer_gausses ) )
then
385 do ishell = 1,
size( fstrsolid%elements(icel)%shell_layer_gausses )
386 if( iselastoplastic( fstrsolid%elements(icel)%shell_layer_gausses(ishell)%pMaterial%mtype ) )
then
387 call updateepstate( fstrsolid%elements(icel)%shell_layer_gausses(ishell) )
389 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_bak = &
390 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain
391 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%stress_bak = &
392 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%stress
393 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_energy_bak = &
394 fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_energy
400 do i = 1, hecmesh%n_node
401 fstrsolid%QFORCE_bak(i) = fstrsolid%QFORCE(i)
406 subroutine update_abort( ic_type, flag, mtype )
407 integer(kind=kint),
intent(in) :: ic_type
408 integer(kind=kint),
intent(in) :: flag
409 integer(kind=kint),
intent(in),
optional :: mtype
413 write(*,*)
'###ERROR### : Element type not supported for static analysis'
414 else if( flag == 2 )
then
415 write(*,*)
'###ERROR### : Element type not supported for nonlinear static analysis'
416 else if( flag == 3 )
then
417 write(*,*)
'###ERROR### : This element is not supported for this material'
419 write(*,*)
' ic_type = ', ic_type
420 if(
present(mtype) )
write(*,*)
' mtype = ', mtype
422 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 the triads contributed by adjacent local shell...
This module provides function to calculate to do updates.
subroutine fstr_updatestate(hecMESH, fstrSOLID, tincr)
Update elastiplastic status.
subroutine fstr_updatenewton(hecMESH, hecMAT, fstrSOLID, fstrPARAM, 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
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.