FrontISTR  5.9.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
27  subroutine fstr_updatenewton ( hecMESH, hecMAT, fstrSOLID, time, tincr,iter, strainEnergy)
28  !=====================================================================*
29  use m_static_lib
30  use m_elemact
32 
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
39 
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
44 
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
54 
55  real(kind=kreal), optional :: strainenergy
56  real(kind=kreal) :: tmp, lambda(1), ddlambda(1)
57  real(kind=kreal) :: ddaux(3,3)
58 
59  ndof = hecmat%NDOF
60  fstrsolid%QFORCE=0.0d0
61  call fstr_ensure_finite_rotation_state( hecmesh, fstrsolid, ndof )
62 
63  tt0 = 0.d0
64  ttn = 0.d0
65  tt = 0.d0
66 
67  ! if initial temperature exists
68  initt = 0
69  if( associated(g_initialcnd) ) then
70  do j=1,size(g_initialcnd)
71  if( g_initialcnd(j)%cond_name=="temperature" ) then
72  initt=j
73  exit
74  endif
75  end do
76  endif
77 
78  ! --------------------------------------------------------------------
79  ! updated
80  ! 1. stress and strain : ep^(k) = ep^(k-1)+dep^(k)
81  ! sgm^(k) = sgm^(k-1)+dsgm^(k)
82  ! 2. Internal Force : Q^(k-1) ( u^(k-1) )
83  ! --------------------------------------------------------------------
84  ! ----------------------------------------------------------------------------------
85  ! calculate the Strain and Stress and Internal Force ( Equivalent Nodal Force )
86  ! ----------------------------------------------------------------------------------
87 
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
94 
95  !element loop
96  !$omp parallel default(none), &
97  !$omp& private(icel,iiS,j,nn,nodLOCAL,i,nbase,ecoord,ddu,du,total_disp,triad_tri,triad_cur, &
98  !$omp& triad_ref, &
99  !$omp& shell_drill, &
100  !$omp& cdsys_ID,coords,thick,qf,isect,ihead,tmp,ndim,ddaux,thick0, &
101  !$omp& lambda,ddlambda), &
102  !$omp& shared(iS,iE,hecMESH,fstrSOLID,ndof,hecMAT,ic_type,fstrPR, &
103  !$omp& strainEnergy,iter,time,tincr,initt,g_InitialCnd), &
104  !$omp& firstprivate(tt0,ttn,tt)
105  !$omp do
106  do icel = is, ie
107 
108  ! ----- nodal coordinate, displacement and temperature
109  iis = hecmesh%elem_node_index(icel-1)
110  nn = hecmesh%elem_node_index(icel)-iis
111  !if( nn>150 ) stop "elemental nodes > 150!"
112 
113  thick = 0.d0
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)))
117  enddo
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
122  enddo
123  end if
124 
125  do j = 1, nn
126  nodlocal(j) = hecmesh%elem_node_item (iis+j)
127  do i = 1, 3
128  ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
129  enddo
130  do i = 1, ndof
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)
134  enddo
135 
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) )
140  else
141  tt0(j) = 0.d0
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))
144  endif
145  ttn(j) = fstrsolid%last_temp( nodlocal(j) )
146  tt(j) = fstrsolid%temperature( nodlocal(j) )
147  endif
148  enddo
149 
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)
155 
156  if (ic_type == 741 .or. ic_type == 743 .or. ic_type == 731) then
157  do j = 1, nn
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))
167  enddo
168  endif
169 
170  ! ===== calculate the Internal Force
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) )
176 
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(:) )
180 
181  else if( ic_type == 361 ) then
182  if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then ! full integration element
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 ! B-bar element
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 ! incompatible element
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 ! F-bar element
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 ! UP element
197  lambda(1) = -0.5d0*fstrsolid%elements(icel)%p(1)
198  CALL update_c3_up &
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)
204  endif
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 ! skip smoothed fem
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) )
209 
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(:) )
213 
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))
218 
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))
223 
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))
230 
231  else if( ic_type == 761 ) then !for shell-solid mixed analysis
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)
235 
236  else if( ic_type == 781 ) then !for shell-solid mixed analysis
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)
240 
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 )
244  call update_c3_vp &
245  ( ic_type, nn, ecoord(:,1:nn), total_disp(1:4,1:nn), du(1:4,1:nn), &
246  fstrsolid%elements(icel)%gausses(:) )
247  qf = 0.0d0
248 
249  else if ( ic_type == 881 .or. ic_type == 891 ) then !for selective es/ns smoothed fem
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) )
252 
253  else
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())
257 
258  endif
259 
260  ! elemact element
261  if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) then
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) )
264  !qf(:) = fstrSOLID%elements(icel)%elemact_coeff*qf(:)
265  end if
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  call hecmw_update_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node, ndof)
297  end subroutine fstr_updatenewton
298 
299 
301  subroutine fstr_updatestate( hecMESH, fstrSOLID, tincr)
302  use m_fstr
303  use m_static_lib
304  use m_elastoplastic
305  use mcreep
306  use mviscoelastic
307  type(hecmwst_local_mesh) :: hecmesh
308  type(fstr_solid) :: fstrSOLID
309  real(kind=kreal) :: tincr
310  integer(kind=kint) :: itype, is, iE, ic_type, icel, ngauss, i
311  integer(kind=kint) :: ishell
312 
313  if( associated( fstrsolid%temperature ) ) then
314  do i = 1, hecmesh%n_node
315  fstrsolid%last_temp(i) = fstrsolid%temperature(i)
316  end do
317  endif
318 
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
326 
327  ngauss = numofquadpoints( ic_type )
328  do icel = is, ie
329  if( iselastoplastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
330  do i = 1, ngauss
331  call updateepstate( fstrsolid%elements(icel)%gausses(i) )
332  enddo
333  elseif( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton ) then
334  if( tincr>0.d0 ) then
335  do i = 1, ngauss
336  call updateviscostate( fstrsolid%elements(icel)%gausses(i) )
337  enddo
338  endif
339  elseif( isviscoelastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
340  if( tincr > 0.d0 ) then
341  do i = 1, ngauss
342  call updateviscoelasticstate( fstrsolid%elements(icel)%gausses(i) )
343  enddo
344  endif
345  endif
346 
347  do i = 1, ngauss
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
351  enddo
352 
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
361  enddo
362  endif
363  enddo
364  enddo
365 
366  do i = 1, hecmesh%n_node
367  fstrsolid%QFORCE_bak(i) = fstrsolid%QFORCE(i)
368  end do
369 
370  end subroutine fstr_updatestate
371 
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
376 
377  if( flag == 1 ) then
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'
383  endif
384  write(*,*) ' ic_type = ', ic_type
385  if( present(mtype) ) write(*,*) ' mtype = ', mtype
386  call hecmw_abort(hecmw_comm_get_comm())
387  end subroutine
388 
389 end module m_fstr_update
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.
Definition: fstr_Update.f90:6
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.
Definition: fstr_Update.f90:28
This module defines common data and basic structures for analysis.
Definition: m_fstr.F90:15
integer(kind=kint), parameter kel361bbar
Definition: m_fstr.F90:80
integer(kind=kint), parameter kel341sesns
Definition: m_fstr.F90:77
integer(kind=kint), parameter kel361up
Definition: m_fstr.F90:83
integer(kind=kint), parameter kel361fi
Definition: m_fstr.F90:79
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords, icel)
This subroutine fetch coords defined by local coordinate system.
Definition: m_fstr.F90:1116
integer(kind=kint), parameter kel361ic
Definition: m_fstr.F90:81
type(fstr_param), target fstrpr
GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP.
Definition: m_fstr.F90:213
integer(kind=kint), parameter kel361fbar
Definition: m_fstr.F90:82
type(tinitialcondition), dimension(:), pointer, save g_initialcnd
Definition: m_fstr.F90:154
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:214
This module provides functions for viscoelastic calculation.
Definition: Viscoelastic.f90:6
subroutine updateviscoelasticstate(gauss)
Update viscoplastic state.