FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
fstr_Update.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 module m_fstr_update
7  use m_fstr
8  implicit none
9 
10  private :: update_abort
11 
12 contains
13 
14  !=====================================================================*
15  ! UPDATE_C3
25  subroutine fstr_updatenewton ( hecMESH, hecMAT, fstrSOLID, time, tincr,iter, strainEnergy)
26  !=====================================================================*
27  use m_static_lib
28 
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
35 
36  integer(kind=kint) :: nodLOCAL(fstrSOLID%max_ncon)
37  real(kind=kreal) :: ecoord(3, fstrsolid%max_ncon)
38  real(kind=kreal) :: thick
39  integer(kind=kint) :: ndof, itype, is, iE, ic_type, nn, icel, iiS, i, j
40 
41  real(kind=kreal) :: total_disp(6, fstrsolid%max_ncon), du(6, fstrsolid%max_ncon), ddu(6, fstrsolid%max_ncon)
42  real(kind=kreal) :: tt(fstrsolid%max_ncon), tt0(fstrsolid%max_ncon), ttn(fstrsolid%max_ncon)
43  real(kind=kreal) :: qf(fstrsolid%max_ncon*6), coords(3, 3)
44  integer :: isect, ihead, cdsys_ID
45  integer :: ndim, initt
46 
47  real(kind=kreal), optional :: strainenergy
48  real(kind=kreal) :: tmp
49  real(kind=kreal) :: ddaux(3,3)
50 
51  ndof = hecmat%NDOF
52  fstrsolid%QFORCE=0.0d0
53 
54  tt0 = 0.d0
55  ttn = 0.d0
56  tt = 0.d0
57 
58  ! if initial temperature exists
59  initt = 0
60  if( associated(g_initialcnd) ) then
61  do j=1,size(g_initialcnd)
62  if( g_initialcnd(j)%cond_name=="temperature" ) then
63  initt=j
64  exit
65  endif
66  end do
67  endif
68 
69  ! --------------------------------------------------------------------
70  ! updated
71  ! 1. stress and strain : ep^(k) = ep^(k-1)+dep^(k)
72  ! sgm^(k) = sgm^(k-1)+dsgm^(k)
73  ! 2. Internal Force : Q^(k-1) ( u^(k-1) )
74  ! --------------------------------------------------------------------
75  ! ----------------------------------------------------------------------------------
76  ! calculate the Strain and Stress and Internal Force ( Equivalent Nodal Force )
77  ! ----------------------------------------------------------------------------------
78 
79  do itype = 1, hecmesh%n_elem_type
80  is = hecmesh%elem_type_index(itype-1)+1
81  ie = hecmesh%elem_type_index(itype )
82  ic_type= hecmesh%elem_type_item(itype)
83  if (hecmw_is_etype_link(ic_type)) cycle
84  if (hecmw_is_etype_patch(ic_type)) cycle
85 
86  !element loop
87  !$omp parallel default(none), &
88  !$omp& private(icel,iiS,j,nn,nodLOCAL,i,ecoord,ddu,du,total_disp, &
89  !$omp& cdsys_ID,coords,thick,qf,isect,ihead,tmp,ndim,ddaux), &
90  !$omp& shared(iS,iE,hecMESH,fstrSOLID,ndof,hecMAT,ic_type,fstrPR, &
91  !$omp& strainEnergy,iter,time,tincr,initt,g_InitialCnd), &
92  !$omp& firstprivate(tt0,ttn,tt)
93  !$omp do
94  do icel = is, ie
95 
96  ! ----- nodal coordinate, displacement and temperature
97  iis = hecmesh%elem_node_index(icel-1)
98  nn = hecmesh%elem_node_index(icel)-iis
99  !if( nn>150 ) stop "elemental nodes > 150!"
100 
101  do j = 1, nn
102  nodlocal(j) = hecmesh%elem_node_item (iis+j)
103  do i = 1, 3
104  ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
105  enddo
106  do i = 1, ndof
107  ddu(i,j) = hecmat%X(ndof*nodlocal(j)+i-ndof)
108  du(i,j) = fstrsolid%dunode(ndof*nodlocal(j)+i-ndof)
109  total_disp(i,j) = fstrsolid%unode(ndof*nodlocal(j)+i-ndof)
110  enddo
111 
112  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
113  if( iselastoplastic(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype) .or. &
114  fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton ) then
115  tt0(j)=fstrsolid%last_temp( nodlocal(j) )
116  else
117  tt0(j) = 0.d0
118  if( hecmesh%hecmw_flag_initcon == 1 ) tt0(j) = hecmesh%node_init_val_item(nodlocal(j))
119  if( initt>0 ) tt0(j) = g_initialcnd(initt)%realval(nodlocal(j))
120  endif
121  ttn(j) = fstrsolid%last_temp( nodlocal(j) )
122  tt(j) = fstrsolid%temperature( nodlocal(j) )
123  endif
124  enddo
125 
126  isect = hecmesh%section_ID(icel)
127  ihead = hecmesh%section%sect_R_index(isect-1)
128  thick = hecmesh%section%sect_R_item(ihead+1)
129  cdsys_id = hecmesh%section%sect_orien_ID(isect)
130  if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
131 
132  ! ===== calculate the Internal Force
133  if( getspacedimension( ic_type ) == 2 ) thick = 1.0d0
134 
135  if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then
136  call update_c2( ic_type,nn,ecoord(1:3,1:nn),fstrsolid%elements(icel)%gausses(:), &
137  thick,fstrsolid%elements(icel)%iset, &
138  total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof), &
139  tt(1:nn), tt0(1:nn), ttn(1:nn) )
140 
141  else if( ic_type == 301 ) then
142  call update_c1( ic_type,nn,ecoord(:,1:nn), thick, total_disp(1:3,1:nn), du(1:3,1:nn), &
143  qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:) )
144 
145  else if( ic_type == 361 ) then
146  if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then ! full integration element
147  call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
148  qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
149  else if( fstrsolid%sections(isect)%elemopt361 == kel361bbar ) then ! B-bar element
150  call update_c3d8bbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
151  qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
152  else if( fstrsolid%sections(isect)%elemopt361 == kel361ic ) then ! incompatible element
153  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,&
154  qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, &
155  fstrsolid%elements(icel)%aux, ddaux(1:3,1:3), tt(1:nn), tt0(1:nn), ttn(1:nn) )
156  fstrsolid%elements(icel)%aux(1:3,1:3) = fstrsolid%elements(icel)%aux(1:3,1:3) + ddaux(1:3,1:3)
157  else if( fstrsolid%sections(isect)%elemopt361 == kel361fbar ) then ! F-bar element
158  call update_c3d8fbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
159  qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
160  endif
161 
162  else if (ic_type == 341 .or. ic_type == 351 .or. ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
163  if( ic_type==341 .and. fstrsolid%sections(isect)%elemopt341 == kel341sesns ) cycle ! skip smoothed fem
164  call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
165  qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
166 
167  else if( ic_type == 511) then
168  call update_connector( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), &
169  qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:), tincr )
170 
171  else if( ic_type == 611) then
172  if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
173  CALL updatest_beam(ic_type, nn, ecoord, total_disp(1:6,1:nn), du(1:6,1:nn), &
174  & hecmesh%section%sect_R_item(ihead+1:), fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof))
175 
176  else if( ic_type == 641 ) then
177  if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
178  call updatest_beam_641(ic_type, nn, ecoord, total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
179  & fstrsolid%elements(icel)%gausses(:), hecmesh%section%sect_R_item(ihead+1:), qf(1:nn*ndof))
180 
181  else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
182  if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
183  call updatest_shell_mitc(ic_type, nn, ndof, ecoord(1:3, 1:nn), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
184  & fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 0)
185 
186  else if( ic_type == 761 ) then !for shell-solid mixed analysis
187  if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
188  call updatest_shell_mitc33(731, 3, 6, ecoord(1:3, 1:3), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
189  & fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 2)
190 
191  else if( ic_type == 781 ) then !for shell-solid mixed analysis
192  if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
193  call updatest_shell_mitc33(741, 4, 6, ecoord(1:3, 1:4), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
194  & fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 1)
195 
196  else if ( ic_type == 3414 ) then
197  if(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype /= incomp_newtonian) &
198  & call update_abort( ic_type, 3, fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype )
199  call update_c3_vp &
200  ( ic_type, nn, ecoord(:,1:nn), total_disp(1:4,1:nn), du(1:4,1:nn), &
201  fstrsolid%elements(icel)%gausses(:) )
202  qf = 0.0d0
203 
204  else if ( ic_type == 881 .or. ic_type == 891 ) then !for selective es/ns smoothed fem
205  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, &
206  qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) )
207 
208  else
209  write(*, *) '###ERROR### : Element type not supported for nonlinear static analysis'
210  write(*, *) ' ic_type = ', ic_type
211  call hecmw_abort(hecmw_comm_get_comm())
212 
213  endif
214 
215  ! ----- calculate the global internal force ( Q(u_{n+1}^{k-1}) )
216  do j = 1, nn
217  do i = 1, ndof
218  !$omp atomic
219  fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i) = fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i)+qf(ndof*(j-1)+i)
220  enddo
221  enddo
222 
223  ! ----- calculate strain energy
224  if(present(strainenergy))then
225  ndim = getspacedimension( fstrsolid%elements(icel)%etype )
226  do j = 1, nn
227  do i = 1, ndim
228  tmp = 0.5d0*( fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i)+qf(ndim*(j-1)+i) )*ddu(i,j)
229  !$omp atomic
230  strainenergy = strainenergy+tmp
231  fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i) = qf(ndim*(j-1)+i)
232  enddo
233  enddo
234  endif
235 
236  enddo ! icel
237  !$omp end do
238  !$omp end parallel
239  enddo ! itype
240 
241  !C
242  !C Update for fstrSOLID%QFORCE
243  !C
244  call hecmw_update_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node, ndof)
245  end subroutine fstr_updatenewton
246 
247 
249  subroutine fstr_updatestate( hecMESH, fstrSOLID, tincr)
250  use m_fstr
251  use m_static_lib
252  use m_elastoplastic
253  use mcreep
254  use mviscoelastic
255  type(hecmwst_local_mesh) :: hecmesh
256  type(fstr_solid) :: fstrSOLID
257  real(kind=kreal) :: tincr
258  integer(kind=kint) :: itype, is, iE, ic_type, icel, ngauss, i
259 
260  if( associated( fstrsolid%temperature ) ) then
261  do i = 1, hecmesh%n_node
262  fstrsolid%last_temp(i) = fstrsolid%temperature(i)
263  end do
264  endif
265 
266  do itype = 1, hecmesh%n_elem_type
267  is = hecmesh%elem_type_index(itype-1) + 1
268  ie = hecmesh%elem_type_index(itype )
269  ic_type= hecmesh%elem_type_item(itype)
270  if( ic_type == 301 ) ic_type = 111
271  if( hecmw_is_etype_link(ic_type) ) cycle
272  if( hecmw_is_etype_patch(ic_type) ) cycle
273 
274  ngauss = numofquadpoints( ic_type )
275  do icel = is, ie
276  if( iselastoplastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
277  do i = 1, ngauss
278  call updateepstate( fstrsolid%elements(icel)%gausses(i) )
279  enddo
280  elseif( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton ) then
281  if( tincr>0.d0 ) then
282  do i = 1, ngauss
283  call updateviscostate( fstrsolid%elements(icel)%gausses(i) )
284  enddo
285  endif
286  elseif( isviscoelastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
287  if( tincr > 0.d0 ) then
288  do i = 1, ngauss
289  call updateviscoelasticstate( fstrsolid%elements(icel)%gausses(i) )
290  enddo
291  endif
292  endif
293 
294  do i = 1, ngauss
295  fstrsolid%elements(icel)%gausses(i)%strain_bak = fstrsolid%elements(icel)%gausses(i)%strain
296  fstrsolid%elements(icel)%gausses(i)%stress_bak = fstrsolid%elements(icel)%gausses(i)%stress
297  fstrsolid%elements(icel)%gausses(i)%strain_energy_bak = fstrsolid%elements(icel)%gausses(i)%strain_energy
298  enddo
299  enddo
300  enddo
301 
302  do i = 1, hecmesh%n_node
303  fstrsolid%QFORCE_bak(i) = fstrsolid%QFORCE(i)
304  end do
305 
306  end subroutine fstr_updatestate
307 
308  subroutine update_abort( ic_type, flag, mtype )
309  integer(kind=kint), intent(in) :: ic_type
310  integer(kind=kint), intent(in) :: flag
311  integer(kind=kint), intent(in), optional :: mtype
312 
313  if( flag == 1 ) then
314  write(*,*) '###ERROR### : Element type not supported for static analysis'
315  else if( flag == 2 ) then
316  write(*,*) '###ERROR### : Element type not supported for nonlinear static analysis'
317  else if( flag == 3 ) then
318  write(*,*) '###ERROR### : This element is not supported for this material'
319  endif
320  write(*,*) ' ic_type = ', ic_type
321  if( present(mtype) ) write(*,*) ' mtype = ', mtype
322  call hecmw_abort(hecmw_comm_get_comm())
323  end subroutine
324 
325 end module m_fstr_update
mcreep::updateviscostate
subroutine updateviscostate(gauss)
Update viscoplastic state.
Definition: creep.f90:214
mviscoelastic
This module provides functions for viscoelastic calculation.
Definition: Viscoelastic.f90:6
m_fstr::kel361fi
integer(kind=kint), parameter kel361fi
Definition: m_fstr.f90:77
m_fstr_update
This module provides function to calculate to do updates.
Definition: fstr_Update.f90:6
m_fstr::kel361bbar
integer(kind=kint), parameter kel361bbar
Definition: m_fstr.f90:78
m_fstr::fstr_solid
Definition: m_fstr.f90:238
mviscoelastic::updateviscoelasticstate
subroutine updateviscoelasticstate(gauss)
Update viscoplastic state.
Definition: Viscoelastic.f90:268
m_fstr::fstrpr
type(fstr_param), target fstrpr
GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP.
Definition: m_fstr.f90:208
m_elastoplastic
This module provide functions for elastoplastic calculation.
Definition: Elastoplastic.f90:6
m_fstr::kel341sesns
integer(kind=kint), parameter kel341sesns
Definition: m_fstr.f90:75
m_elastoplastic::updateepstate
subroutine, public updateepstate(gauss)
Clear elatoplastic state.
Definition: Elastoplastic.f90:876
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
m_fstr::g_initialcnd
type(tinitialcondition), dimension(:), pointer, save g_initialcnd
Definition: m_fstr.f90:151
m_fstr_update::fstr_updatenewton
subroutine fstr_updatenewton(hecMESH, hecMAT, fstrSOLID, time, tincr, iter, strainEnergy)
Update displacement, stress, strain and internal forces.
Definition: fstr_Update.f90:26
m_fstr::get_coordsys
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords)
This subroutine fetch coords defined by local coordinate system.
Definition: m_fstr.f90:1073
m_fstr::kel361fbar
integer(kind=kint), parameter kel361fbar
Definition: m_fstr.f90:80
m_fstr::kel361ic
integer(kind=kint), parameter kel361ic
Definition: m_fstr.f90:79
m_static_lib
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
m_fstr_update::fstr_updatestate
subroutine fstr_updatestate(hecMESH, fstrSOLID, tincr)
Update elastiplastic status.
Definition: fstr_Update.f90:250
mcreep
This module provides functions for creep calculation.
Definition: creep.f90:6