FrontISTR  5.2.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 !-------------------------------------------------------------------------------
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(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
40 
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
45 
46  real(kind=kreal), optional :: strainenergy
47  real(kind=kreal) :: tmp
48  real(kind=kreal) :: ddaux(3,3)
49 
50  ndof = hecmat%NDOF
51  fstrsolid%QFORCE=0.0d0
52 
53  tt0 = 0.d0
54  ttn = 0.d0
55  tt = 0.d0
56 
57  ! if initial temperature exists
58  initt = 0
59  if( associated(g_initialcnd) ) then
60  do j=1,size(g_initialcnd)
61  if( g_initialcnd(j)%cond_name=="temperature" ) then
62  initt=j
63  exit
64  endif
65  end do
66  endif
67 
68  ! --------------------------------------------------------------------
69  ! updated
70  ! 1. stress and strain : ep^(k) = ep^(k-1)+dep^(k)
71  ! sgm^(k) = sgm^(k-1)+dsgm^(k)
72  ! 2. Internal Force : Q^(k-1) ( u^(k-1) )
73  ! --------------------------------------------------------------------
74  ! ----------------------------------------------------------------------------------
75  ! calculate the Strain and Stress and Internal Force ( Equivalent Nodal Force )
76  ! ----------------------------------------------------------------------------------
77 
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!"
86 
87  !element loop
88  !$omp parallel default(none), &
89  !$omp& private(icel,iiS,j,nodLOCAL,i,ecoord,ddu,du,total_disp, &
90  !$omp& cdsys_ID,coords,thick,qf,isect,ihead,tmp,ndim,ddaux), &
91  !$omp& shared(iS,iE,hecMESH,nn,fstrSOLID,ndof,hecMAT,ic_type,fstrPR, &
92  !$omp& strainEnergy,iter,time,tincr,initt,g_InitialCnd), &
93  !$omp& firstprivate(tt0,ttn,tt)
94  !$omp do
95  do icel = is, ie
96 
97  ! ----- nodal coordinate
98  iis = hecmesh%elem_node_index(icel-1)
99 
100  do j = 1, nn
101  nodlocal(j) = hecmesh%elem_node_item (iis+j)
102  do i = 1, 3
103  ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
104  enddo
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) )
109  else
110  tt0(j) = 0.d0
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))
113  endif
114  ttn(j) = fstrsolid%last_temp( nodlocal(j) )
115  tt(j) = fstrsolid%temperature( nodlocal(j) )
116  endif
117  enddo
118 
119  ! nodal displacement
120  do j = 1, nn
121  nodlocal(j) = hecmesh%elem_node_item (iis+j)
122  do i = 1, ndof
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)
126  enddo
127  enddo
128 
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)
132 
133  ! ===== calculate the Internal Force
134  if( getspacedimension( ic_type ) == 2 ) thick = 1.0d0
135 
136  if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then
137 
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) )
143  else
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) )
147  endif
148 
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(:) )
155 
156  else if( ic_type == 361 ) then
157 
158  if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then ! full integration element
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) )
162  else
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 )
165  endif
166  else if( fstrsolid%sections(isect)%elemopt361 == kel361bbar ) then ! B-bar element
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) )
170  else
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 )
173  endif
174  else if( fstrsolid%sections(isect)%elemopt361 == kel361ic ) then ! incompatible element
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) )
179  else
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) )
183  endif
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 ! F-bar element
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) )
189  else
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 )
192  endif
193  endif
194 
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) )
199  else
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 )
202  endif
203 
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))
210 
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))
217 
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)
225 
226  else if( ic_type == 761 ) then !for shell-solid mixed analysis
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)
233 
234  else if( ic_type == 781 ) then !for shell-solid mixed analysis
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)
241 
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())
247  endif
248  call update_c3_vp &
249  ( ic_type, nn, ecoord(:,1:nn), total_disp(1:4,1:nn), du(1:4,1:nn), &
250  fstrsolid%elements(icel)%gausses(:) )
251  qf = 0.0d0
252 
253  ! else if ( ic_type==731) then
254  ! call UPDATE_S3(xx,yy,zz,ee,pp,thick,local_stf)
255  ! call fstr_local_stf_restore_temp(local_stf, nn*ndof, stiffness)
256  ! else if ( ic_type==741) then
257  ! call UPDATE_S4(xx,yy,zz,ee,pp,thick,local_stf)
258  ! call fstr_local_stf_restore_temp(local_stf, nn*ndof, stiffness)
259 
260  else
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())
264 
265  endif
266 
267  ! ----- calculate the global internal force ( Q(u_{n+1}^{k-1}) )
268  do j = 1, nn
269  do i = 1, ndof
270  !$omp atomic
271  fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i) = fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i)+qf(ndof*(j-1)+i)
272  enddo
273  enddo
274 
275  ! ----- calculate strain energy
276  if(present(strainenergy))then
277  ndim = getspacedimension( fstrsolid%elements(icel)%etype )
278  do j = 1, nn
279  do i = 1, ndim
280  tmp = 0.5d0*( fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i)+qf(ndim*(j-1)+i) )*ddu(i,j)
281  !$omp atomic
282  strainenergy = strainenergy+tmp
283  fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i) = qf(ndim*(j-1)+i)
284  enddo
285  enddo
286  endif
287 
288  enddo ! icel
289  !$omp end do
290  !$omp end parallel
291  enddo ! itype
292 
293  !C
294  !C Update for fstrSOLID%QFORCE
295  !C
296  if( ndof==3 ) then
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)
304  endif
305 
306  end subroutine fstr_updatenewton
307 
308 
310  subroutine fstr_updatestate( hecMESH, fstrSOLID, tincr)
311  use m_fstr
312  use m_static_lib
313  use m_elastoplastic
314  use mcreep
315  use mviscoelastic
316  type(hecmwst_local_mesh) :: hecmesh
317  type(fstr_solid) :: fstrSOLID
318  real(kind=kreal) :: tincr
319  integer(kind=kint) :: itype, is, iE, ic_type, icel, ngauss, i
320 
321  if( associated( fstrsolid%temperature ) ) then
322  do i = 1, hecmesh%n_node
323  fstrsolid%last_temp(i) = fstrsolid%temperature(i)
324  end do
325  endif
326 
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
334 
335  ngauss = numofquadpoints( ic_type )
336  do icel = is, ie
337  if( iselastoplastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
338  do i = 1, ngauss
339  call updateepstate( fstrsolid%elements(icel)%gausses(i) )
340  enddo
341  elseif( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton ) then
342  if( tincr>0.d0 ) then
343  do i = 1, ngauss
344  call updateviscostate( fstrsolid%elements(icel)%gausses(i) )
345  enddo
346  endif
347  elseif( isviscoelastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
348  if( tincr > 0.d0 ) then
349  do i = 1, ngauss
350  call updateviscoelasticstate( fstrsolid%elements(icel)%gausses(i) )
351  enddo
352  endif
353  endif
354 
355  do i = 1, ngauss
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
358  enddo
359  enddo
360  enddo
361  end subroutine fstr_updatestate
362 
363  subroutine update_abort( ic_type, flag )
364  integer(kind=kint), intent(in) :: ic_type
365  integer(kind=kint), intent(in) :: flag
366 
367  if( flag == 1 ) then
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'
371  endif
372  write(*,*) ' ic_type = ', ic_type
373  call hecmw_abort(hecmw_comm_get_comm())
374  end subroutine
375 
376 end module m_fstr_update
This module provide functions for elastoplastic calculation.
subroutine updateepstate(gauss)
Clear elatoplastic state.
This module provides function to calcualte to do updates.
Definition: fstr_Update.f90:6
subroutine fstr_updatestate(hecMESH, fstrSOLID, tincr)
Update elastiplastic status.
subroutine fstr_updatenewton(hecMESH, hecMAT, fstrSOLID, time, tincr, iter, strainEnergy)
変位/応力・ひずみ/内力のアップデート
Definition: fstr_Update.f90:26
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint), parameter kel361bbar
Definition: m_fstr.f90:69
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords)
This subroutine fetch coords defined by local coordinate system.
Definition: m_fstr.f90:1032
integer(kind=kint), parameter kel361fi
section control
Definition: m_fstr.f90:68
integer(kind=kint), parameter kel361ic
Definition: m_fstr.f90:70
type(fstr_param), target fstrpr
GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP.
Definition: m_fstr.f90:190
integer(kind=kint), parameter kel361fbar
Definition: m_fstr.f90:71
type(tinitialcondition), dimension(:), pointer, save g_initialcnd
Definition: m_fstr.f90:135
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
This module provides functions for creep calculation.
Definition: creep.f90:6
subroutine updateviscostate(gauss)
Update viscoplastic state.
Definition: creep.f90:217
This module provides functions for viscoelastic calculation.
Definition: Viscoelastic.f90:6
subroutine updateviscoelasticstate(gauss)
Update viscoplastic state.