FrontISTR  6.0.0
Large-scale structural analysis program with finit element method
fstr_NodalStress.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_nodalstress
7  use m_fstr
8  use mmechgauss, only: fstr_element_average_plstrain, fstr_shell_output_point_count
9 
10  implicit none
11  private :: nodalstress_inv3, nodalstress_inv2, inverse_func
12 contains
13 
15  subroutine fstr_get_shell_gauss_output(hecMESH, fstrSOLID, icel, k, strain, stress)
16  use m_static_lib, only: elementstress_shell_mitc
17  use elementinfo, only: numofshellthicknessquadpoints, getshellthicknessquadpoint
18  type(hecmwst_local_mesh), intent(in) :: hecMESH
19  type(fstr_solid), intent(inout) :: fstrSOLID
20  integer(kind=kint), intent(in) :: icel, k
21  real(kind=kreal), intent(out) :: strain(6), stress(6)
22  integer(kind=kint) :: etype, nn, j, inode, isect, ihead, nlayer, nthick, ig, ilayer, ithick, cdsys_ID
23  real(kind=kreal) :: ecoord(3,9), edisp(6,9), point_strain(1,6), point_stress(1,6), thick, zeta, coords(3,3)
24 
25  strain = 0.d0
26  stress = 0.d0
27  if( k>fstr_shell_output_point_count(fstrsolid%elements(icel)) ) return
28  if( associated(fstrsolid%elements(icel)%shell_layer_gausses) ) then
29  strain = fstrsolid%elements(icel)%shell_layer_gausses(k)%strain_out
30  stress = fstrsolid%elements(icel)%shell_layer_gausses(k)%stress_out
31  return
32  endif
33  etype = hecmesh%elem_type(icel)
34  nn = hecmesh%elem_node_index(icel)-hecmesh%elem_node_index(icel-1)
35  nlayer = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
36  nthick = numofshellthicknessquadpoints(etype)
37  ig = (k-1)/(nlayer*nthick)+1
38  ilayer = mod((k-1)/nthick, nlayer)+1
39  ithick = mod(k-1, nthick)+1
40  call getshellthicknessquadpoint(etype, ithick, zeta)
41  do j=1,nn
42  inode = hecmesh%elem_node_item(hecmesh%elem_node_index(icel-1)+j)
43  ecoord(:,j) = hecmesh%node(3*inode-2:3*inode)
44  edisp(:,j) = fstrsolid%unode(6*inode-5:6*inode)
45  enddo
46  isect = hecmesh%section_ID(icel)
47  ihead = hecmesh%section%sect_R_index(isect-1)
48  thick = hecmesh%section%sect_R_item(ihead+1)
49  cdsys_id = hecmesh%section%sect_orien_ID(isect)
50  coords = 0.0d0
51  if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords, icel)
52  call elementstress_shell_mitc(etype, nn, 6, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses, &
53  edisp(:,1:nn), point_strain, point_stress, thick, zeta, ilayer, surface_gauss_index=ig, cdsys_id=cdsys_id, coords=coords)
54  strain = point_strain(1,:)
55  stress = point_stress(1,:)
56  end subroutine fstr_get_shell_gauss_output
57 
59  !----------------------------------------------------------------------*
60  subroutine fstr_nodalstress3d( hecMESH, fstrSOLID )
61  !----------------------------------------------------------------------*
62  use m_static_lib
63  use m_fstr_nodalkinematics, only: fstr_ensure_finite_rotation_state
64  type(hecmwst_local_mesh) :: hecmesh
65  type(fstr_solid) :: fstrSOLID
66  real(kind=kreal), pointer :: tnstrain(:), testrain(:), yield_ratio(:)
67  integer(kind=kint), pointer :: is_rot(:)
68  !C** local variables
69  integer(kind=kint) :: itype, icel, ic, is, iE, jS, i, j, ic_type, nn, ni, ID_area, nbase
70  integer(kind=kint) :: nodlocal(20), ntemp
71  integer(kind=kint), allocatable :: nnumber(:)
72  real(kind=kreal) :: estrain(6), estress(6), naturalcoord(3)
73  real(kind=kreal) :: enqm(12)
74  real(kind=kreal) :: ndstrain(20,6), ndstress(20,6), tdstrain(20,6)
75  real(kind=kreal) :: ecoord(3, 20), edisp(60), tt(20), t0(20)
76  real(kind=kreal) :: beam_disp(6,2)
77  real(kind=kreal) :: triad_cur(9,4), triad_ref(9,4)
78  real(kind=kreal), allocatable :: func(:,:), inv_func(:,:)
79 
80  !C** Shell33 variables
81  integer(kind=kint) :: isect, ihead, ntot_lyr, nlyr, flag33, truss, cdsys_ID
82  real(kind=kreal) :: mid_strain(9,6), mid_stress(9,6)
83  real(kind=kreal) :: thick, coords(3,3)
84  call fstr_solid_phys_clear(fstrsolid)
85 
86  allocate( nnumber(hecmesh%n_node) )
87  if( .not. associated(fstrsolid%is_rot) ) allocate( fstrsolid%is_rot(hecmesh%n_node) )
88  !allocate( fstrSOLID%yield_ratio(hecMESH%n_elem) )
89  nnumber = 0
90  fstrsolid%is_rot = 0
91  call fstr_ensure_finite_rotation_state(hecmesh, fstrsolid, 3)
92  !fstrSOLID%yield_ratio = 0.0d0
93 
94  tnstrain => fstrsolid%tnstrain
95  testrain => fstrsolid%testrain
96  is_rot => fstrsolid%is_rot
97  yield_ratio => fstrsolid%yield_ratio
98 
99  if( associated(tnstrain) ) tnstrain = 0.0d0
100 
101  !C** setting
102  ntot_lyr = fstrsolid%max_lyr
103  flag33 = fstrsolid%is_33shell
104  truss = fstrsolid%is_33beam
105 
106  !C +-------------------------------+
107  !C | according to ELEMENT TYPE |
108  !C +-------------------------------+
109  do itype = 1, hecmesh%n_elem_type
110  is = hecmesh%elem_type_index(itype-1) + 1
111  ie = hecmesh%elem_type_index(itype )
112  ic_type = hecmesh%elem_type_item(itype)
113  if( ic_type == fe_tet10nc ) ic_type = fe_tet10n
114  if( .not. (hecmw_is_etype_solid(ic_type) .or. ic_type == 781 &
115  & .or. ic_type == 761 .or. ic_type == fe_beam341 ) ) cycle
116  !C** set number of nodes and shape function
117  nn = hecmw_get_max_node( ic_type )
118  ni = numofquadpoints( ic_type )
119  allocate( func(ni,nn), inv_func(nn,ni) )
120  tdstrain = 0.0d0
121  if( ic_type == fe_tet10n ) then
122  ic = hecmw_get_max_node( fe_tet4n )
123  do i = 1, ni
124  call getquadpoint( ic_type, i, naturalcoord )
125  call getshapefunc( fe_tet4n, naturalcoord, func(i,1:ic) )
126  enddo
127  call inverse_func( ic, func, inv_func )
128  else if( ic_type == fe_hex8n ) then
129  do i = 1, ni
130  call getquadpoint( ic_type, i, naturalcoord )
131  call getshapefunc( ic_type, naturalcoord, func(i,1:nn) )
132  enddo
133  call inverse_func( ni, func, inv_func )
134  else if( ic_type == fe_prism15n ) then
135  ic = 0
136  do i = 1, ni
137  if( i==1 .or. i==2 .or. i==3 .or. i==7 .or. i==8 .or. i==9 ) then
138  ic = ic + 1
139  call getquadpoint( ic_type, i, naturalcoord )
140  call getshapefunc( fe_prism6n, naturalcoord, func(ic,1:6) )
141  endif
142  enddo
143  call inverse_func( ic, func, inv_func )
144  ni = ic
145  else if( ic_type == fe_hex20n ) then
146  ic = 0
147  do i = 1, ni
148  if( i==1 .or. i==3 .or. i==7 .or. i==9 .or. &
149  i==19 .or. i==21 .or. i==25 .or. i==27 ) then
150  ic = ic + 1
151  call getquadpoint( ic_type, i, naturalcoord )
152  call getshapefunc( fe_hex8n, naturalcoord, func(ic,1:8) )
153  endif
154  enddo
155  call inverse_func( ic, func, inv_func )
156  ni = ic
157  endif
158  !C** element loop
159  do icel = is, ie
160  js = hecmesh%elem_node_index(icel-1)
161  id_area = hecmesh%elem_ID(icel*2)
162  isect= hecmesh%section_ID(icel)
163  ihead = hecmesh%section%sect_R_index(isect-1)
164  thick = hecmesh%section%sect_R_item(ihead+1)
165  cdsys_id = hecmesh%section%sect_orien_ID(isect)
166  coords = 0.0d0
167  if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords, icel)
168  !initialize
169  enqm = 0.0d0
170  estrain = 0.0d0
171  estress = 0.0d0
172  ndstrain = 0.0d0
173  ndstress = 0.0d0
174  !if( ID_area == hecMESH%my_rank ) then
175 
176  !--- calculate nodal and elemental value
177  if( ic_type == 641 ) then
178  do j = 1, 4
179  nodlocal(j) = hecmesh%elem_node_item(js+j)
180  ecoord(1:3,j) = hecmesh%node(3*nodlocal(j)-2:3*nodlocal(j))
181  edisp(3*j-2:3*j) = fstrsolid%unode(3*nodlocal(j)-2:3*nodlocal(j))
182  end do
183  beam_disp(1:3,1:2) = reshape(edisp(1:6), (/ 3, 2 /))
184  beam_disp(4:6,1:2) = reshape(edisp(7:12), (/ 3, 2 /))
185  ntemp = 0
186  if( associated( fstrsolid%temperature ) ) then
187  ntemp = 1
188  do j = 1, 4
189  nodlocal(j) = hecmesh%elem_node_item(js+j)
190  t0(j) = fstrsolid%last_temp( nodlocal(j) )
191  tt(j) = fstrsolid%temperature( nodlocal(j) )
192  end do
193  end if
194  call nodalstress_beam( 2, ecoord(1:3,1:2), fstrsolid%elements(icel)%gausses, &
195  hecmesh%section%sect_R_item(ihead+1:), beam_disp, &
196  ndstrain(1:2,1:6), ndstress(1:2,1:6), fstrsolid%sections(isect)%elemopt611 )
197  call elementalstress_beam( fstrsolid%elements(icel)%gausses, estrain, estress, enqm )
198  fstrsolid%ENQM(icel*12-11:icel*12) = enqm(1:12)
199 
200 
201  elseif( ic_type == 781) then
202  triad_cur = 0.0d0
203  triad_ref = 0.0d0
204  do j = 1, 4
205  nodlocal(j ) = hecmesh%elem_node_item(js+j )
206  nodlocal(j+4) = hecmesh%elem_node_item(js+j+4)
207  is_rot(nodlocal(j+4)) = 1
208  ecoord(1:3,j ) = hecmesh%node(3*nodlocal(j )-2:3*nodlocal(j ))
209  ecoord(1:3,j+4) = hecmesh%node(3*nodlocal(j+4)-2:3*nodlocal(j+4))
210  edisp(6*j-5:6*j-3) = fstrsolid%unode(3*nodlocal(j )-2:3*nodlocal(j ))
211  edisp(6*j-2:6*j ) = fstrsolid%unode(3*nodlocal(j+4)-2:3*nodlocal(j+4))
212  nbase = 9*(nodlocal(j+4)-1)
213  if( associated(fstrsolid%shell_triad) ) triad_cur(1:9,j) = fstrsolid%shell_triad(nbase+1:nbase+9)
214  if( associated(fstrsolid%shell_ref_triad) ) triad_ref(1:9,j) = fstrsolid%shell_ref_triad(nbase+1:nbase+9)
215  enddo
216  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
217  mid_strain = 0.0d0
218  mid_stress = 0.0d0
219  do nlyr=1,ntot_lyr
220  if( associated(fstrsolid%shell_triad) ) then
221  call elementstress_shell_mitc( 741, 4, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
222  ndstrain(1:4,1:6), ndstress(1:4,1:6), thick, 1.0d0, nlyr, &
223  ndtriad=triad_cur, ndreftriad=triad_ref, cdsys_id=cdsys_id, coords=coords)
224  else
225  call elementstress_shell_mitc( 741, 4, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
226  ndstrain(1:4,1:6), ndstress(1:4,1:6), thick, 1.0d0, nlyr, cdsys_id=cdsys_id, coords=coords)
227  endif
228  call fstr_stress_add_shelllyr(4,fstrsolid,icel,nodlocal,nlyr,ndstrain(1:4,1:6),ndstress(1:4,1:6),1)
229  call fstr_add_shell_mid(fstrsolid%elements(icel),nlyr,4,ndstrain,ndstress,mid_strain,mid_stress)
230  !minus section
231  if( associated(fstrsolid%shell_triad) ) then
232  call elementstress_shell_mitc( 741, 4, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
233  ndstrain(1:4,1:6), ndstress(1:4,1:6), thick,-1.0d0, nlyr, &
234  ndtriad=triad_cur, ndreftriad=triad_ref, cdsys_id=cdsys_id, coords=coords)
235  else
236  call elementstress_shell_mitc( 741, 4, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
237  ndstrain(1:4,1:6), ndstress(1:4,1:6), thick,-1.0d0, nlyr, cdsys_id=cdsys_id, coords=coords)
238  endif
239  call fstr_stress_add_shelllyr(4,fstrsolid,icel,nodlocal,nlyr,ndstrain(1:4,1:6),ndstress(1:4,1:6),-1)
240  call fstr_add_shell_mid(fstrsolid%elements(icel),nlyr,4,ndstrain,ndstress,mid_strain,mid_stress)
241  enddo
242  ndstrain(1:4,1:6) = mid_strain(1:4,1:6)
243  ndstress(1:4,1:6) = mid_stress(1:4,1:6)
244  call fstr_getavg_shell(fstrsolid,icel,estrain,estress)
245 
246  elseif( ic_type == 761) then
247  triad_cur = 0.0d0
248  triad_ref = 0.0d0
249  do j = 1, 3
250  nodlocal(j ) = hecmesh%elem_node_item(js+j )
251  nodlocal(j+3) = hecmesh%elem_node_item(js+j+3)
252  is_rot(nodlocal(j+3)) = 1
253  ecoord(1:3,j ) = hecmesh%node(3*nodlocal(j )-2:3*nodlocal(j ))
254  ecoord(1:3,j+3) = hecmesh%node(3*nodlocal(j+3)-2:3*nodlocal(j+3))
255  edisp(6*j-5:6*j-3) = fstrsolid%unode(3*nodlocal(j )-2:3*nodlocal(j ))
256  edisp(6*j-2:6*j ) = fstrsolid%unode(3*nodlocal(j+3)-2:3*nodlocal(j+3))
257  nbase = 9*(nodlocal(j+3)-1)
258  if( associated(fstrsolid%shell_triad) ) triad_cur(1:9,j) = fstrsolid%shell_triad(nbase+1:nbase+9)
259  if( associated(fstrsolid%shell_ref_triad) ) triad_ref(1:9,j) = fstrsolid%shell_ref_triad(nbase+1:nbase+9)
260  enddo
261  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
262  mid_strain = 0.0d0
263  mid_stress = 0.0d0
264  do nlyr=1,ntot_lyr
265  if( associated(fstrsolid%shell_triad) ) then
266  call elementstress_shell_mitc( 731, 3, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
267  ndstrain(1:3,1:6), ndstress(1:3,1:6), thick, 1.0d0, nlyr, &
268  ndtriad=triad_cur(1:9,1:3), ndreftriad=triad_ref(1:9,1:3), cdsys_id=cdsys_id, coords=coords)
269  else
270  call elementstress_shell_mitc( 731, 3, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
271  ndstrain(1:3,1:6), ndstress(1:3,1:6), thick, 1.0d0, nlyr, cdsys_id=cdsys_id, coords=coords)
272  endif
273  call fstr_stress_add_shelllyr(3,fstrsolid,icel,nodlocal,nlyr,ndstrain(1:3,1:6),ndstress(1:3,1:6),1)
274  call fstr_add_shell_mid(fstrsolid%elements(icel),nlyr,3,ndstrain,ndstress,mid_strain,mid_stress)
275  !minus section
276  if( associated(fstrsolid%shell_triad) ) then
277  call elementstress_shell_mitc( 731, 3, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
278  ndstrain(1:3,1:6), ndstress(1:3,1:6), thick,-1.0d0, nlyr, &
279  ndtriad=triad_cur(1:9,1:3), ndreftriad=triad_ref(1:9,1:3), cdsys_id=cdsys_id, coords=coords)
280  else
281  call elementstress_shell_mitc( 731, 3, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
282  ndstrain(1:3,1:6), ndstress(1:3,1:6), thick,-1.0d0, nlyr, cdsys_id=cdsys_id, coords=coords)
283  endif
284  call fstr_stress_add_shelllyr(3,fstrsolid,icel,nodlocal,nlyr,ndstrain(1:3,1:6),ndstress(1:3,1:6),-1)
285  call fstr_add_shell_mid(fstrsolid%elements(icel),nlyr,3,ndstrain,ndstress,mid_strain,mid_stress)
286  enddo
287  ndstrain(1:3,1:6) = mid_strain(1:3,1:6)
288  ndstress(1:3,1:6) = mid_stress(1:3,1:6)
289  call fstr_getavg_shell(fstrsolid,icel,estrain,estress)
290 
291  else if( ic_type == 301 ) then
292  call nodalstress_c1( nn, fstrsolid%elements(icel)%gausses, &
293  ndstrain(1:nn,1:6), ndstress(1:nn,1:6) )
294  call elementstress_c1( fstrsolid%elements(icel)%gausses, estrain, estress )
295 
296  else if( ic_type == fe_tet10n .or. ic_type == fe_hex8n .or. &
297  ic_type == fe_prism15n .or. ic_type == fe_hex20n ) then
298  call nodalstress_inv3( ic_type, ni, fstrsolid%elements(icel)%gausses, &
299  inv_func, ndstrain(1:nn,1:6), ndstress(1:nn,1:6), &
300  tdstrain(1:nn,1:6) )
301  call elementstress_c3( ic_type, fstrsolid%elements(icel)%gausses, estrain, estress )
302 
303  else if ( ic_type == 881 .or. ic_type == 891 ) then !for selective es/ns smoothed fem
304  cycle
305  else
306  if( ic_type == 341 .and. fstrsolid%sections(isect)%elemopt341 == kel341sesns ) cycle
307 
308  call nodalstress_c3( ic_type, nn, fstrsolid%elements(icel)%gausses, &
309  ndstrain(1:nn,1:6), ndstress(1:nn,1:6) )
310  !call NodalStress_C3( ic_type, nn, fstrSOLID%elements(icel)%gausses, &
311  ! ndstrain(1:nn,1:6), ndstress(1:nn,1:6), tdstrain(1:nn,1:6) )
312  call elementstress_c3( ic_type, fstrsolid%elements(icel)%gausses, estrain, estress )
313 
314  endif
315 
316  !ADD VALUE and Count node
317  do j = 1, nn
318  ic = hecmesh%elem_node_item(js+j)
319  fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) + ndstrain(j,1:6)
320  fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) + ndstress(j,1:6)
321  if( associated(tnstrain) )then
322  tnstrain(6*(ic-1)+1:6*(ic-1)+6) = tnstrain(6*(ic-1)+1:6*(ic-1)+6) + tdstrain(j,1:6)
323  endif
324  nnumber(ic) = nnumber(ic) + 1
325  enddo
326 
327  fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) + estrain(1:6)
328  fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) + estress(1:6)
329 
330  !endif
331  enddo
332  deallocate( func, inv_func )
333  enddo
334 
335  !C** calculate nodal stress and strain
336  do i = 1, hecmesh%n_node
337  if( nnumber(i) == 0 ) cycle
338  fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
339  fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
340  if( associated(tnstrain) )then
341  tnstrain(6*(i-1)+1:6*(i-1)+6) = tnstrain(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
342  endif
343  enddo
344 
345  if( fstrsolid%is_smoothing_active ) call fstr_nodalstress3d_c3d4_sesns( &
346  & hecmesh, fstrsolid, nnumber, fstrsolid%STRAIN, fstrsolid%STRESS, fstrsolid%ESTRAIN, fstrsolid%ESTRESS )
347 
348  ! copy owner nodal stress/strain to external(ghost) slots so DD boundary MISES/principal match serial
349  call hecmw_update_r( hecmesh, fstrsolid%STRESS, hecmesh%n_node, 6 )
350  call hecmw_update_r( hecmesh, fstrsolid%STRAIN, hecmesh%n_node, 6 )
351  if( associated(tnstrain) ) call hecmw_update_r( hecmesh, tnstrain, hecmesh%n_node, 6 )
352 
353  if( flag33 == 1 )then
354  do nlyr = 1, ntot_lyr
355  do i = 1, hecmesh%n_node
356  if( nnumber(i) == 0 ) cycle
357  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
358  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
359  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
360  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
361  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
362  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
363  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
364  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
365  fstrsolid%SHELL%LAYER(nlyr)%PLUS%MISES(i) = &
366  & get_mises(fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6))
367  fstrsolid%SHELL%LAYER(nlyr)%MINUS%MISES(i) = &
368  & get_mises(fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6))
369  enddo
370  enddo
371  endif
372 
373  !C** calculate von MISES stress
374  do i = 1, hecmesh%n_node
375  fstrsolid%MISES(i) = get_mises(fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6))
376  enddo
377  do i = 1, hecmesh%n_elem
378  fstrsolid%EMISES(i) = get_mises(fstrsolid%ESTRESS(6*(i-1)+1:6*(i-1)+6))
379  enddo
380 
381  !C** calculate Elemental Plastic Strain
382  do i = 1, hecmesh%n_elem
383  fstrsolid%EPLSTRAIN(i) = fstr_element_average_plstrain(fstrsolid%elements(i))
384  enddo
385 
386  if( flag33 == 1 )then
387  if( fstrsolid%output_ctrl(3)%outinfo%on(27) .or. fstrsolid%output_ctrl(4)%outinfo%on(27) ) then
388  do nlyr = 1, ntot_lyr
389  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL%LAYER(nlyr)%PLUS)
390  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL%LAYER(nlyr)%MINUS)
391  enddo
392  endif
393  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL)
394  else
395  call make_principal(fstrsolid, hecmesh, fstrsolid%SOLID)
396  endif
397 
398  deallocate( nnumber )
399 
400  end subroutine fstr_nodalstress3d
401 
402  integer(kind=kint) function search_idx_senes( irow, asect, nid, sid )
403  integer(kind=kint), allocatable, intent(in) :: irow(:)
404  integer(kind=kint), allocatable, intent(in) :: asect(:)
405  integer(kind=kint), intent(in) :: nid
406  integer(kind=kint), intent(in) :: sid
407 
408  integer(kind=kint) :: i
409 
410  search_idx_senes = -1
411  do i=irow(nid-1)+1,irow(nid)
412  if( asect(i) == sid ) then
413  search_idx_senes = i
414  return
415  end if
416  end do
417 
418  end function
419 
420  subroutine fstr_nodalstress3d_c3d4_sesns( hecMESH, fstrSOLID, nnumber, &
421  Nodal_STRAIN, Nodal_STRESS, Elemental_STRAIN, Elemental_STRESS )
422  type(hecmwst_local_mesh),intent(in) :: hecMESH
423  type(fstr_solid),intent(inout) :: fstrSOLID
424  integer(kind=kint), allocatable, intent(inout) :: nnumber(:)
425  real(kind=kreal), pointer, intent(inout) :: nodal_strain(:)
426  real(kind=kreal), pointer, intent(inout) :: nodal_stress(:)
427  real(kind=kreal), pointer, intent(inout) :: elemental_strain(:)
428  real(kind=kreal), pointer, intent(inout) :: elemental_stress(:)
429 
430  integer(kind=kint) :: itype, iS, iE, jS, ic_type, icel, i, j, isect
431  integer(kind=kint) :: nsize, nid(2), idx(2), nd
432  integer(kind=kint) :: nnode, nlen
433  type(hecmwst_varray_int), allocatable :: nodal_sections(:)
434  real(kind=kreal) :: tmpval(6), nsecdup
435  integer(kind=kint), allocatable :: irow(:), asect(:)
436  real(kind=kreal), allocatable :: stress_hyd(:), strain_hyd(:)
437  real(kind=kreal), allocatable :: stress_dev(:)
438  real(kind=kreal), allocatable :: plstrain_dev(:)
439  real(kind=kreal) :: stress_hyd_ndave(6), strain_hyd_ndave(6)
440  real(kind=kreal) :: stress_dev_ndave(6)
441  real(kind=kreal), allocatable :: n_dup_dev(:), n_dup_hyd(:)
442  real(kind=kreal) :: edstrain(6), edstress(6)
443  real(kind=kreal) :: edplstrain
444 
445  nnode = hecmesh%n_node
446  nsize = size(nodal_strain)
447 
448  ! create section info at node
449  call hecmw_varray_int_initialize_all( nodal_sections, nnode, 2 )
450  do itype = 1, hecmesh%n_elem_type
451  ic_type = hecmesh%elem_type_item(itype)
452  if( ic_type /= 341 ) cycle
453 
454  is = hecmesh%elem_type_index(itype-1) + 1
455  ie = hecmesh%elem_type_index(itype )
456 
457  do icel=is,ie
458  isect= hecmesh%section_ID(icel)
459  if( fstrsolid%sections(isect)%elemopt341 /= kel341sesns ) cycle
460  js = hecmesh%elem_node_index(icel-1)
461  do i=1,4
462  nd = hecmesh%elem_node_item(js+i)
463  call hecmw_varray_int_add_if_not_exits( nodal_sections(nd), isect )
464  end do
465  end do
466  enddo
467 
468  ! create CRS arrays of nodal stress/strain with different sections
469  allocate(irow(0:nnode))
470  irow(0) = 0
471  do i=1,nnode
472  irow(i) = irow(i-1)+hecmw_varray_int_get_nitem(nodal_sections(i))
473  end do
474  nlen = irow(nnode)
475 
476  allocate(asect(nlen))
477  do i=1,nnode
478  if( irow(i-1) == irow(i) ) cycle
479  call hecmw_varray_int_get_item_all( nodal_sections(i), asect(irow(i-1)+1:irow(i)) )
480  end do
481 
482  ! add stress/strain from smoothed elements
483  allocate(stress_hyd(6*nlen), strain_hyd(6*nlen))
484  allocate(stress_dev(6*nlen))
485  allocate(plstrain_dev(nlen))
486  allocate(n_dup_dev(nlen),n_dup_hyd(nlen))
487 
488  stress_hyd(:) = 0.d0
489  strain_hyd(:) = 0.d0
490  stress_dev(:) = 0.d0
491  plstrain_dev(:) = 0.d0
492  n_dup_hyd(:) = 0.d0
493  n_dup_dev(:) = 0.d0
494  do itype = 1, hecmesh%n_elem_type
495  ic_type = hecmesh%elem_type_item(itype)
496  if( ic_type /= 881 .and. ic_type /= 891 ) cycle
497 
498  is = hecmesh%elem_type_index(itype-1) + 1
499  ie = hecmesh%elem_type_index(itype )
500 
501  do icel=is,ie
502  js = hecmesh%elem_node_index(icel-1)
503  isect= hecmesh%section_ID(icel)
504  if( ic_type == 881 ) then
505  nid(1) = hecmesh%elem_node_item(js+1)
506  idx(1) = search_idx_senes( irow, asect, nid(1), isect )
507 
508  !strain
509  strain_hyd(6*idx(1)-5:6*idx(1)) = fstrsolid%elements(icel)%gausses(1)%strain_out(1:6)
510  !stress
511  stress_hyd(6*idx(1)-5:6*idx(1)) = fstrsolid%elements(icel)%gausses(1)%stress_out(1:6)
512  !number of duplication
513  n_dup_hyd(idx(1)) = n_dup_hyd(idx(1)) + 1.d0
514  else if( ic_type == 891 ) then
515  nid(1:2) = hecmesh%elem_node_item(js+1:js+2)
516  idx(1) = search_idx_senes( irow, asect, nid(1), isect )
517  idx(2) = search_idx_senes( irow, asect, nid(2), isect )
518 
519  !stress
520  tmpval(1:6) = fstrsolid%elements(icel)%gausses(1)%stress_out(1:6)
521  stress_dev(6*idx(1)-5:6*idx(1)) = stress_dev(6*idx(1)-5:6*idx(1)) + tmpval(1:6)
522  stress_dev(6*idx(2)-5:6*idx(2)) = stress_dev(6*idx(2)-5:6*idx(2)) + tmpval(1:6)
523  !plastic strain
524  plstrain_dev(idx(1)) = plstrain_dev(idx(1)) + fstrsolid%elements(icel)%gausses(1)%plstrain
525  plstrain_dev(idx(2)) = plstrain_dev(idx(2)) + fstrsolid%elements(icel)%gausses(1)%plstrain
526  !number of duplication
527  n_dup_dev(idx(1)) = n_dup_dev(idx(1)) + 1.d0
528  n_dup_dev(idx(2)) = n_dup_dev(idx(2)) + 1.d0
529  end if
530  end do
531  enddo
532 
533  do i=1,nnode
534  if( irow(i-1) == irow(i) ) cycle
535  do j=irow(i-1)+1,irow(i)
536  if( n_dup_dev(j) < 1.0d-8 ) cycle
537  stress_dev(6*j-5:6*j) = stress_dev(6*j-5:6*j)/n_dup_dev(j)
538  plstrain_dev(j) = plstrain_dev(j)/n_dup_dev(j)
539  end do
540  end do
541 
542  ! average at node for nodal output
543  do i=1,nnode
544  if( irow(i-1) == irow(i) ) cycle
545  strain_hyd_ndave(:) = 0.d0
546  stress_hyd_ndave(:) = 0.d0
547  stress_dev_ndave(:) = 0.d0
548  do j=irow(i-1)+1,irow(i)
549  strain_hyd_ndave(1:6) = strain_hyd_ndave(1:6) + strain_hyd(6*j-5:6*j)
550  stress_hyd_ndave(1:6) = stress_hyd_ndave(1:6) + stress_hyd(6*j-5:6*j)
551  stress_dev_ndave(1:6) = stress_dev_ndave(1:6) + stress_dev(6*j-5:6*j)
552  end do
553  nsecdup = dble(irow(i)-irow(i-1))
554  strain_hyd_ndave(1:6) = strain_hyd_ndave(1:6)/nsecdup
555  stress_hyd_ndave(1:6) = stress_hyd_ndave(1:6)/nsecdup
556  stress_dev_ndave(1:6) = stress_dev_ndave(1:6)/nsecdup
557 
558  if( nnumber(i) == 0 ) then
559  nodal_strain(6*i-5:6*i) = strain_hyd_ndave(1:6)
560  nodal_stress(6*i-5:6*i) = stress_hyd_ndave(1:6)+stress_dev_ndave(1:6)
561  else
562  nodal_strain(6*i-5:6*i) = 0.5d0*(nodal_strain(6*i-5:6*i)+strain_hyd_ndave(1:6))
563  nodal_stress(6*i-5:6*i) = 0.5d0*(nodal_stress(6*i-5:6*i)+stress_hyd_ndave(1:6)+stress_dev_ndave(1:6))
564  endif
565  end do
566 
567  ! ELEMENTAL STRAIN and STRESS
568  do itype = 1, hecmesh%n_elem_type
569  ic_type = hecmesh%elem_type_item(itype)
570  if( ic_type /= 341 ) cycle
571 
572  is = hecmesh%elem_type_index(itype-1) + 1
573  ie = hecmesh%elem_type_index(itype )
574 
575  do icel=is,ie
576  isect= hecmesh%section_ID(icel)
577  if( fstrsolid%sections(isect)%elemopt341 /= kel341sesns ) cycle
578  js = hecmesh%elem_node_index(icel-1)
579  edstrain(1:6) = 0.d0
580  edstress(1:6) = 0.d0
581  edplstrain = 0.d0
582  do i=1,4
583  nd = hecmesh%elem_node_item(js+i)
584  idx(1) = search_idx_senes( irow, asect, hecmesh%elem_node_item(js+i), isect )
585  edstrain(1:6) = edstrain(1:6) + strain_hyd(6*idx(1)-5:6*idx(1))
586  edstress(1:6) = edstress(1:6) + stress_hyd(6*idx(1)-5:6*idx(1)) + stress_dev(6*idx(1)-5:6*idx(1))
587  edplstrain = edplstrain + plstrain_dev(idx(1))
588  end do
589  edstrain(1:6) = 0.25d0*edstrain(1:6)
590  edstress(1:6) = 0.25d0*edstress(1:6)
591  edplstrain = 0.25d0*edplstrain
592 
593  elemental_strain(6*(icel-1)+1:6*(icel-1)+6) = elemental_strain(6*(icel-1)+1:6*(icel-1)+6) + edstrain(1:6)
594  elemental_stress(6*(icel-1)+1:6*(icel-1)+6) = elemental_stress(6*(icel-1)+1:6*(icel-1)+6) + edstress(1:6)
595 
596  fstrsolid%elements(icel)%gausses(1)%strain_out(1:6) = elemental_strain(6*(icel-1)+1:6*(icel-1)+6)
597  fstrsolid%elements(icel)%gausses(1)%stress_out(1:6) = elemental_stress(6*(icel-1)+1:6*(icel-1)+6)
598  fstrsolid%elements(icel)%gausses(1)%plstrain = edplstrain
599  end do
600  enddo
601 
602  deallocate(stress_hyd, strain_hyd)
603  deallocate(stress_dev, plstrain_dev)
604  deallocate(n_dup_dev, n_dup_hyd)
605 
606  end subroutine
607 
608  subroutine fstr_stress_add_shelllyr(nn,fstrSOLID,icel,nodLOCAL,nlyr,strain,stress,flag)
609  implicit none
610  type(fstr_solid) :: fstrsolid
611  integer(kind=kint) :: nodlocal(20)
612  integer(kind=kint) :: nn, i, j, k, m, nlyr, icel, flag
613  real(kind=kreal) :: strain(nn, 6), stress(nn, 6)
614  type(fstr_solid_physic_val), pointer :: layer => null()
615 
616  do j = 1, nn
617  i = nodlocal(j)
618  m = nodlocal(j+nn)
619  if(flag == 1)then
620  layer => fstrsolid%SHELL%LAYER(nlyr)%PLUS
621  elseif(flag == -1)then
622  layer => fstrsolid%SHELL%LAYER(nlyr)%MINUS
623  endif
624  do k = 1, 6
625  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + strain(j,k)
626  layer%STRAIN(6*(m-1)+k) = layer%STRAIN(6*(m-1)+k) + strain(j,k)
627  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + stress(j,k)
628  layer%STRESS(6*(m-1)+k) = layer%STRESS(6*(m-1)+k) + stress(j,k)
629  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + strain(j,k)/nn
630  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + stress(j,k)/nn
631  enddo
632  enddo
633  end subroutine fstr_stress_add_shelllyr
634 
637  subroutine fstr_add_shell_mid(element,nlyr,nn,strain,stress,mid_strain,mid_stress)
638  implicit none
639  type(telement), intent(in) :: element
640  integer(kind=kint), intent(in) :: nlyr, nn
641  real(kind=kreal), intent(in) :: strain(:,:), stress(:,:)
642  real(kind=kreal), intent(inout) :: mid_strain(:,:), mid_stress(:,:)
643  real(kind=kreal) :: weight
644 
645  weight = 0.5d0*element%gausses(1)%pMaterial%shell_var(nlyr)%weight
646  mid_strain(1:nn,1:6) = mid_strain(1:nn,1:6) + weight*strain(1:nn,1:6)
647  mid_stress(1:nn,1:6) = mid_stress(1:nn,1:6) + weight*stress(1:nn,1:6)
648  end subroutine fstr_add_shell_mid
649 
652  subroutine fstr_getavg_shell(fstrSOLID,icel,estrain,estress)
653  implicit none
654  type (fstr_solid) :: fstrsolid
655  integer(kind=kint) :: j, nlyr, icel, ntot_lyr
656  real(kind=kreal) :: estrain(6), estress(6), weight
657  type(fstr_solid_physic_val), pointer :: layer => null()
658 
659  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
660  estrain = 0.0d0
661  estress = 0.0d0
662 
663  do nlyr = 1, ntot_lyr
664  layer => fstrsolid%SHELL%LAYER(nlyr)
665  weight = fstrsolid%elements(icel)%gausses(1)%pMaterial%shell_var(nlyr)%weight
666  do j = 1, 6
667  estrain(j) = estrain(j) &
668  & + weight*(0.5d0*layer%PLUS%ESTRAIN(6*(icel-1)+j) + 0.5d0*layer%MINUS%ESTRAIN(6*(icel-1)+j))
669  estress(j) = estress(j) &
670  & + weight*(0.5d0*layer%PLUS%ESTRESS(6*(icel-1)+j) + 0.5d0*layer%MINUS%ESTRESS(6*(icel-1)+j))
671  enddo
672  enddo
673  call set_shell_layer_surface_results(fstrsolid%elements(icel), fstrsolid, icel)
674  call get_shell_layer_gauss_average(fstrsolid%elements(icel), estrain, estress)
675  end subroutine fstr_getavg_shell
676 
677  subroutine set_shell_layer_surface_results(element, fstrSOLID, icel)
678  implicit none
679  type(telement), intent(in) :: element
680  type(fstr_solid), intent(inout) :: fstrSOLID
681  integer(kind=kint), intent(in) :: icel
682  integer(kind=kint) :: ilayer
683  real(kind=kreal) :: estrain(6), estress(6), eplstrain
684 
685  if( .not. associated(element%shell_layer_gausses) ) return
686  if( element%shell_nlayer <= 0 ) return
687 
688  do ilayer = 1, element%shell_nlayer
689  call get_shell_layer_surface_average(element, ilayer, 1, estrain, estress, eplstrain)
690  fstrsolid%SHELL%LAYER(ilayer)%PLUS%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) = estrain(1:6)
691  fstrsolid%SHELL%LAYER(ilayer)%PLUS%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) = estress(1:6)
692  fstrsolid%SHELL%LAYER(ilayer)%PLUS%EPLSTRAIN(icel) = eplstrain
693  call get_shell_layer_surface_average(element, ilayer, -1, estrain, estress, eplstrain)
694  fstrsolid%SHELL%LAYER(ilayer)%MINUS%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) = estrain(1:6)
695  fstrsolid%SHELL%LAYER(ilayer)%MINUS%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) = estress(1:6)
696  fstrsolid%SHELL%LAYER(ilayer)%MINUS%EPLSTRAIN(icel) = eplstrain
697  enddo
698  end subroutine set_shell_layer_surface_results
699 
700  subroutine get_shell_layer_surface_average(element, ilayer, flag, estrain, estress, eplstrain)
701  implicit none
702  type(telement), intent(in) :: element
703  integer(kind=kint), intent(in) :: ilayer, flag
704  real(kind=kreal), intent(out) :: estrain(6), estress(6), eplstrain
705  integer(kind=kint) :: ig, ithick, ishell, ierr, surface_ithick, npoint
706  real(kind=kreal) :: zeta_layer, weight, surface_zeta
707 
708  estrain(1:6) = 0.0d0
709  estress(1:6) = 0.0d0
710  eplstrain = 0.0d0
711  if( .not. associated(element%shell_layer_gausses) ) return
712  if( .not. associated(element%gausses) ) return
713  if( element%shell_nlayer <= 0 .or. element%shell_nthick <= 0 ) return
714 
715  surface_ithick = 0
716  if( flag > 0 ) then
717  surface_zeta = -huge(1.0d0)
718  else
719  surface_zeta = huge(1.0d0)
720  endif
721 
722  do ithick = 1, element%shell_nthick
723  call fstr_shell_layer_quadrature(element, ilayer, ithick, zeta_layer, weight, ierr)
724  if( ierr /= 0 ) cycle
725  if( flag > 0 ) then
726  if( zeta_layer > surface_zeta ) then
727  surface_zeta = zeta_layer
728  surface_ithick = ithick
729  endif
730  else
731  if( zeta_layer < surface_zeta ) then
732  surface_zeta = zeta_layer
733  surface_ithick = ithick
734  endif
735  endif
736  enddo
737  if( surface_ithick <= 0 ) return
738 
739  npoint = 0
740  do ig = 1, size(element%gausses)
741  ishell = fstr_shell_layer_gauss_index(element, ig, ilayer, surface_ithick)
742  if( ishell <= 0 ) cycle
743  estrain(1:6) = estrain(1:6) + element%shell_layer_gausses(ishell)%strain_out(1:6)
744  estress(1:6) = estress(1:6) + element%shell_layer_gausses(ishell)%stress_out(1:6)
745  eplstrain = eplstrain + element%shell_layer_gausses(ishell)%plstrain
746  npoint = npoint + 1
747  enddo
748  if( npoint > 0 ) then
749  estrain(1:6) = estrain(1:6) / npoint
750  estress(1:6) = estress(1:6) / npoint
751  eplstrain = eplstrain / npoint
752  endif
753  end subroutine get_shell_layer_surface_average
754 
755  subroutine get_shell_layer_gauss_average(element, estrain, estress)
756  implicit none
757  type(telement), intent(in) :: element
758  real(kind=kreal), intent(inout) :: estrain(6), estress(6)
759  integer(kind=kint) :: ig, ilayer, ithick, ishell
760  real(kind=kreal) :: weight, total_weight
761  real(kind=kreal) :: avg_strain(6), avg_stress(6)
762 
763  if( .not. associated(element%shell_layer_gausses) ) return
764  if( .not. associated(element%gausses) ) return
765  if( element%shell_nlayer <= 0 .or. element%shell_nthick <= 0 ) return
766  if( size(element%gausses) <= 0 ) return
767  if( .not. associated(element%gausses(1)%pMaterial) ) return
768  if( .not. associated(element%gausses(1)%pMaterial%shell_var) ) return
769  if( element%shell_nlayer > size(element%gausses(1)%pMaterial%shell_var) ) return
770 
771  avg_strain(1:6) = 0.0d0
772  avg_stress(1:6) = 0.0d0
773  total_weight = 0.0d0
774 
775  ! Average history points equally within each layer and weight layers by thickness.
776  do ig = 1, size(element%gausses)
777  do ilayer = 1, element%shell_nlayer
778  weight = element%gausses(1)%pMaterial%shell_var(ilayer)%weight
779  do ithick = 1, element%shell_nthick
780  ishell = fstr_shell_layer_gauss_index(element, ig, ilayer, ithick)
781  if( ishell <= 0 ) cycle
782  avg_strain(1:6) = avg_strain(1:6) &
783  + element%shell_layer_gausses(ishell)%strain_out(1:6) * weight
784  avg_stress(1:6) = avg_stress(1:6) &
785  + element%shell_layer_gausses(ishell)%stress_out(1:6) * weight
786  total_weight = total_weight + weight
787  enddo
788  enddo
789  enddo
790 
791  if( total_weight > 0.0d0 ) then
792  estrain(1:6) = avg_strain(1:6) / total_weight
793  estress(1:6) = avg_stress(1:6) / total_weight
794  endif
795  end subroutine get_shell_layer_gauss_average
796 
798  subroutine nodalstress_shellj2(element, nn, fstrSOLID, icel, nodLOCAL, strain, stress, estrain, estress)
799  use elementinfo, only: getquadpoint, getshapefunc
800  implicit none
801  type(telement), intent(in) :: element
802  integer(kind=kint), intent(in) :: nn, icel, nodlocal(:)
803  type(fstr_solid), intent(inout) :: fstrSOLID
804  real(kind=kreal), intent(out) :: strain(:,:), stress(:,:), estrain(6), estress(6)
805  real(kind=kreal) :: func(nn,nn), inv_func(nn,nn), coord(2), shape(nn)
806  real(kind=kreal) :: gp_strain(nn,6), gp_stress(nn,6), nod_strain(nn,6), nod_stress(nn,6)
807  real(kind=kreal) :: avg_strain(nn,6), avg_stress(nn,6), zeta, weight, total_weight
808  integer(kind=kint) :: ig, ilayer, ithick, ishell, ierr, j, i
809  type(fstr_solid_physic_val), pointer :: layer
810 
811  do ig = 1, nn
812  call getquadpoint(element%etype, ig, coord)
813  call getshapefunc(element%etype, coord, shape)
814  func(ig,:) = shape
815  enddo
816  call inverse_func(nn, func, inv_func)
817  avg_strain = 0.0d0
818  avg_stress = 0.0d0
819  estrain = 0.0d0
820  estress = 0.0d0
821  total_weight = 0.0d0
822 
823  do ilayer = 1, element%shell_nlayer
824  do ithick = 1, element%shell_nthick
825  call fstr_shell_layer_quadrature(element, ilayer, ithick, zeta, weight, ierr)
826  if( ierr /= 0 ) stop 'Invalid shell thickness quadrature'
827  do ig = 1, nn
828  ishell = fstr_shell_layer_gauss_index(element, ig, ilayer, ithick)
829  gp_strain(ig,:) = element%shell_layer_gausses(ishell)%strain_out(1:6)
830  gp_stress(ig,:) = element%shell_layer_gausses(ishell)%stress_out(1:6)
831  enddo
832  avg_strain = avg_strain+weight*gp_strain
833  avg_stress = avg_stress+weight*gp_stress
834  total_weight = total_weight+weight
835 
836  ! Layer +/- results use the outermost thickness integration points.
837  if( ithick /= 1 .and. ithick /= element%shell_nthick ) cycle
838  nod_strain = matmul(inv_func, gp_strain)
839  nod_stress = matmul(inv_func, gp_stress)
840  if( ithick == 1 ) then
841  layer => fstrsolid%SHELL%LAYER(ilayer)%MINUS
842  else
843  layer => fstrsolid%SHELL%LAYER(ilayer)%PLUS
844  endif
845  do j = 1, nn
846  i = nodlocal(j)
847  layer%STRAIN(6*i-5:6*i) = layer%STRAIN(6*i-5:6*i)+nod_strain(j,:)
848  layer%STRESS(6*i-5:6*i) = layer%STRESS(6*i-5:6*i)+nod_stress(j,:)
849  enddo
850  enddo
851  enddo
852  strain = matmul(inv_func, avg_strain)/total_weight
853  stress = matmul(inv_func, avg_stress)/total_weight
854  call set_shell_layer_surface_results(element, fstrsolid, icel)
855  call get_shell_layer_gauss_average(element, estrain, estress)
856  end subroutine nodalstress_shellj2
857 
858  !----------------------------------------------------------------------*
859  subroutine nodalstress_inv3( etype, ni, gausses, func, edstrain, edstress, tdstrain )
860  !----------------------------------------------------------------------*
861  use mmechgauss
862  integer(kind=kint) :: etype, ni
863  type(tgaussstatus) :: gausses(:)
864  real(kind=kreal) :: func(:, :), edstrain(:, :), edstress(:, :), tdstrain(:, :)
865  integer :: i, j, k, ic
866 
867  edstrain = 0.0d0
868  edstress = 0.0d0
869  tdstrain = 0.0d0
870 
871  if( etype == fe_hex8n ) then
872  do i = 1, ni
873  do j = 1, ni
874  do k = 1, 6
875  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
876  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
877  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
878  enddo
879  enddo
880  enddo
881  else if( etype == fe_tet10n ) then
882  do i = 1, ni
883  do j = 1, ni
884  do k = 1, 6
885  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
886  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
887  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
888  enddo
889  enddo
890  enddo
891  edstrain(5,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
892  edstress(5,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
893  tdstrain(5,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
894  edstrain(6,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
895  edstress(6,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
896  tdstrain(6,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
897  edstrain(7,1:6) = ( edstrain(3,1:6) + edstrain(1,1:6) ) / 2.0
898  edstress(7,1:6) = ( edstress(3,1:6) + edstress(1,1:6) ) / 2.0
899  tdstrain(7,1:6) = ( tdstrain(3,1:6) + tdstrain(1,1:6) ) / 2.0
900  edstrain(8,1:6) = ( edstrain(1,1:6) + edstrain(4,1:6) ) / 2.0
901  edstress(8,1:6) = ( edstress(1,1:6) + edstress(4,1:6) ) / 2.0
902  tdstrain(8,1:6) = ( tdstrain(1,1:6) + tdstrain(4,1:6) ) / 2.0
903  edstrain(9,1:6) = ( edstrain(2,1:6) + edstrain(4,1:6) ) / 2.0
904  edstress(9,1:6) = ( edstress(2,1:6) + edstress(4,1:6) ) / 2.0
905  tdstrain(9,1:6) = ( tdstrain(2,1:6) + tdstrain(4,1:6) ) / 2.0
906  edstrain(10,1:6) = ( edstrain(3,1:6) + edstrain(4,1:6) ) / 2.0
907  edstress(10,1:6) = ( edstress(3,1:6) + edstress(4,1:6) ) / 2.0
908  tdstrain(10,1:6) = ( tdstrain(3,1:6) + tdstrain(4,1:6) ) / 2.0
909  else if( etype == fe_prism15n ) then
910  do i = 1, ni
911  ic = 0
912  do j = 1, numofquadpoints(etype)
913  if( j==1 .or. j==2 .or. j==3 .or. j==7 .or. j==8 .or. j==9 ) then
914  ic = ic + 1
915  do k = 1, 6
916  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
917  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
918  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
919  enddo
920  endif
921  enddo
922  enddo
923  edstrain(7,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
924  edstress(7,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
925  tdstrain(7,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
926  edstrain(8,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
927  edstress(8,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
928  tdstrain(8,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
929  edstrain(9,1:6) = ( edstrain(3,1:6) + edstrain(1,1:6) ) / 2.0
930  edstress(9,1:6) = ( edstress(3,1:6) + edstress(1,1:6) ) / 2.0
931  tdstrain(9,1:6) = ( tdstrain(3,1:6) + tdstrain(1,1:6) ) / 2.0
932  edstrain(10,1:6) = ( edstrain(4,1:6) + edstrain(5,1:6) ) / 2.0
933  edstress(10,1:6) = ( edstress(4,1:6) + edstress(5,1:6) ) / 2.0
934  tdstrain(10,1:6) = ( tdstrain(4,1:6) + tdstrain(5,1:6) ) / 2.0
935  edstrain(11,1:6) = ( edstrain(5,1:6) + edstrain(6,1:6) ) / 2.0
936  edstress(11,1:6) = ( edstress(5,1:6) + edstress(6,1:6) ) / 2.0
937  tdstrain(11,1:6) = ( tdstrain(5,1:6) + tdstrain(6,1:6) ) / 2.0
938  edstrain(12,1:6) = ( edstrain(6,1:6) + edstrain(4,1:6) ) / 2.0
939  edstress(12,1:6) = ( edstress(6,1:6) + edstress(4,1:6) ) / 2.0
940  tdstrain(12,1:6) = ( tdstrain(6,1:6) + tdstrain(4,1:6) ) / 2.0
941  edstrain(13,1:6) = ( edstrain(1,1:6) + edstrain(4,1:6) ) / 2.0
942  edstress(13,1:6) = ( edstress(1,1:6) + edstress(4,1:6) ) / 2.0
943  tdstrain(13,1:6) = ( tdstrain(1,1:6) + tdstrain(4,1:6) ) / 2.0
944  edstrain(14,1:6) = ( edstrain(2,1:6) + edstrain(5,1:6) ) / 2.0
945  edstress(14,1:6) = ( edstress(2,1:6) + edstress(5,1:6) ) / 2.0
946  tdstrain(14,1:6) = ( tdstrain(2,1:6) + tdstrain(5,1:6) ) / 2.0
947  edstrain(15,1:6) = ( edstrain(3,1:6) + edstrain(6,1:6) ) / 2.0
948  edstress(15,1:6) = ( edstress(3,1:6) + edstress(6,1:6) ) / 2.0
949  tdstrain(15,1:6) = ( tdstrain(3,1:6) + tdstrain(6,1:6) ) / 2.0
950  else if( etype == fe_hex20n ) then
951  do i = 1, ni
952  ic = 0
953  do j = 1, numofquadpoints(etype)
954  if( j==1 .or. j==3 .or. j==7 .or. j==9 .or. &
955  j==19 .or. j==21 .or. j==25 .or. j==27 ) then
956  ic = ic + 1
957  do k = 1, 6
958  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
959  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
960  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
961  enddo
962  endif
963  enddo
964  enddo
965  edstrain(9,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
966  edstress(9,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
967  tdstrain(9,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
968  edstrain(10,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
969  edstress(10,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
970  tdstrain(10,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
971  edstrain(11,1:6) = ( edstrain(3,1:6) + edstrain(4,1:6) ) / 2.0
972  edstress(11,1:6) = ( edstress(3,1:6) + edstress(4,1:6) ) / 2.0
973  tdstrain(11,1:6) = ( tdstrain(3,1:6) + tdstrain(4,1:6) ) / 2.0
974  edstrain(12,1:6) = ( edstrain(4,1:6) + edstrain(1,1:6) ) / 2.0
975  edstress(12,1:6) = ( edstress(4,1:6) + edstress(1,1:6) ) / 2.0
976  tdstrain(12,1:6) = ( tdstrain(4,1:6) + tdstrain(1,1:6) ) / 2.0
977  edstrain(13,1:6) = ( edstrain(5,1:6) + edstrain(6,1:6) ) / 2.0
978  edstress(13,1:6) = ( edstress(5,1:6) + edstress(6,1:6) ) / 2.0
979  tdstrain(13,1:6) = ( tdstrain(5,1:6) + tdstrain(6,1:6) ) / 2.0
980  edstrain(14,1:6) = ( edstrain(6,1:6) + edstrain(7,1:6) ) / 2.0
981  edstress(14,1:6) = ( edstress(6,1:6) + edstress(7,1:6) ) / 2.0
982  tdstrain(14,1:6) = ( tdstrain(6,1:6) + tdstrain(7,1:6) ) / 2.0
983  edstrain(15,1:6) = ( edstrain(7,1:6) + edstrain(8,1:6) ) / 2.0
984  edstress(15,1:6) = ( edstress(7,1:6) + edstress(8,1:6) ) / 2.0
985  tdstrain(15,1:6) = ( tdstrain(7,1:6) + tdstrain(8,1:6) ) / 2.0
986  edstrain(16,1:6) = ( edstrain(8,1:6) + edstrain(5,1:6) ) / 2.0
987  edstress(16,1:6) = ( edstress(8,1:6) + edstress(5,1:6) ) / 2.0
988  tdstrain(16,1:6) = ( tdstrain(8,1:6) + tdstrain(5,1:6) ) / 2.0
989  edstrain(17,1:6) = ( edstrain(1,1:6) + edstrain(5,1:6) ) / 2.0
990  edstress(17,1:6) = ( edstress(1,1:6) + edstress(5,1:6) ) / 2.0
991  tdstrain(17,1:6) = ( tdstrain(1,1:6) + tdstrain(5,1:6) ) / 2.0
992  edstrain(18,1:6) = ( edstrain(2,1:6) + edstrain(6,1:6) ) / 2.0
993  edstress(18,1:6) = ( edstress(2,1:6) + edstress(6,1:6) ) / 2.0
994  tdstrain(18,1:6) = ( tdstrain(2,1:6) + tdstrain(6,1:6) ) / 2.0
995  edstrain(19,1:6) = ( edstrain(3,1:6) + edstrain(7,1:6) ) / 2.0
996  edstress(19,1:6) = ( edstress(3,1:6) + edstress(7,1:6) ) / 2.0
997  tdstrain(19,1:6) = ( tdstrain(3,1:6) + tdstrain(7,1:6) ) / 2.0
998  edstrain(20,1:6) = ( edstrain(4,1:6) + edstrain(8,1:6) ) / 2.0
999  edstress(20,1:6) = ( edstress(4,1:6) + edstress(8,1:6) ) / 2.0
1000  tdstrain(20,1:6) = ( tdstrain(4,1:6) + tdstrain(8,1:6) ) / 2.0
1001  endif
1002  end subroutine nodalstress_inv3
1003 
1004  function get_mises(s)
1005  implicit none
1006  real(kind=kreal) :: get_mises, s(1:6)
1007  real(kind=kreal) :: s11, s22, s33, s12, s23, s13, ps, smises
1008 
1009  s11 = s(1)
1010  s22 = s(2)
1011  s33 = s(3)
1012  s12 = s(4)
1013  s23 = s(5)
1014  s13 = s(6)
1015  ps = ( s11 + s22 + s33 ) / 3.0d0
1016  smises = 0.5d0 * ( (s11-ps)**2 + (s22-ps)**2 + (s33-ps)**2 ) + s12**2 + s23**2 + s13**2
1017  get_mises = dsqrt( 3.0d0 * smises )
1018 
1019  end function get_mises
1020 
1022  !----------------------------------------------------------------------*
1023  subroutine fstr_nodalstress2d( hecMESH, fstrSOLID )
1024  !----------------------------------------------------------------------*
1025  use m_static_lib
1026  type (hecmwST_local_mesh) :: hecMESH
1027  type (fstr_solid) :: fstrSOLID
1028  real(kind=kreal), pointer :: tnstrain(:), testrain(:)
1029  !C** local variables
1030  integer(kind=kint) :: itype, icel, ic, is, iE, jS, i, j, ic_type, nn, ni, ID_area
1031  real(kind=kreal) :: estrain(4), estress(4), naturalcoord(4)
1032  real(kind=kreal) :: edstrain(8,4), edstress(8,4), tdstrain(8,4)
1033  real(kind=kreal) :: s11, s22, s12, smises
1034  real(kind=kreal), allocatable :: func(:,:), inv_func(:,:)
1035  integer(kind=kint), allocatable :: nnumber(:)
1036 
1037  tnstrain => fstrsolid%tnstrain
1038  testrain => fstrsolid%testrain
1039  call fstr_solid_phys_clear(fstrsolid)
1040 
1041  allocate( nnumber(hecmesh%n_node) )
1042  if( .not. associated(fstrsolid%is_rot) ) allocate( fstrsolid%is_rot(hecmesh%n_node) )
1043  nnumber = 0
1044  fstrsolid%is_rot = 0
1045 
1046  if( associated(tnstrain) ) tnstrain = 0.0d0
1047 
1048  !C +-------------------------------+
1049  !C | according to ELEMENT TYPE |
1050  !C +-------------------------------+
1051  do itype = 1, hecmesh%n_elem_type
1052  is = hecmesh%elem_type_index(itype-1) + 1
1053  ie = hecmesh%elem_type_index(itype )
1054  ic_type = hecmesh%elem_type_item(itype)
1055  if( .not. hecmw_is_etype_surface(ic_type) ) cycle
1056  !C** set number of nodes and shape function
1057  nn = hecmw_get_max_node( ic_type )
1058  ni = numofquadpoints( ic_type )
1059  allocate( func(ni,nn), inv_func(nn,ni) )
1060  tdstrain = 0.0d0
1061  if( ic_type == fe_tri6n ) then
1062  ic = hecmw_get_max_node( fe_tri3n )
1063  do i = 1, ni
1064  call getquadpoint( ic_type, i, naturalcoord )
1065  call getshapefunc( fe_tri3n, naturalcoord, func(i,1:ic) )
1066  enddo
1067  call inverse_func( ic, func, inv_func )
1068  else if( ic_type == fe_quad4n ) then
1069  do i = 1, ni
1070  call getquadpoint( ic_type, i, naturalcoord )
1071  call getshapefunc( ic_type, naturalcoord, func(i,1:nn) )
1072  enddo
1073  call inverse_func( ni, func, inv_func )
1074  else if( ic_type == fe_quad8n ) then
1075  ic = 0
1076  do i = 1, ni
1077  if( i==1 .or. i==3 .or. i==7 .or. i==9 ) then
1078  ic = ic + 1
1079  call getquadpoint( ic_type, i, naturalcoord )
1080  call getshapefunc( fe_quad4n, naturalcoord, func(ic,1:4) )
1081  endif
1082  enddo
1083  call inverse_func( ic, func, inv_func )
1084  ni = ic
1085  endif
1086  !C** element loop
1087  do icel = is, ie
1088  js = hecmesh%elem_node_index(icel-1)
1089  id_area = hecmesh%elem_ID(icel*2)
1090  !--- calculate nodal stress and strain
1091  if( ic_type == fe_tri6n .or. ic_type == fe_quad4n .or. ic_type == fe_quad8n ) then
1092  call nodalstress_inv2( ic_type, ni, fstrsolid%elements(icel)%gausses, &
1093  inv_func, edstrain(1:nn,1:4), edstress(1:nn,1:4), &
1094  tdstrain(1:nn,1:4) )
1095  else
1096  call nodalstress_c2( ic_type, nn, fstrsolid%elements(icel)%gausses, &
1097  edstrain(1:nn,1:4), edstress(1:nn,1:4) )
1098  ! call NodalStress_C2( ic_type, nn, fstrSOLID%elements(icel)%gausses, &
1099  ! edstrain(1:nn,1:4), edstress(1:nn,1:4), tdstrain(1:nn,1:4) )
1100  endif
1101  do j = 1, nn
1102  ic = hecmesh%elem_node_item(js+j)
1103  fstrsolid%STRAIN(3*ic-2) = fstrsolid%STRAIN(3*ic-2) + edstrain(j,1)
1104  fstrsolid%STRAIN(3*ic-1) = fstrsolid%STRAIN(3*ic-1) + edstrain(j,2)
1105  fstrsolid%STRAIN(3*ic-0) = fstrsolid%STRAIN(3*ic-0) + edstrain(j,3)
1106  fstrsolid%STRESS(3*ic-2) = fstrsolid%STRESS(3*ic-2) + edstress(j,1)
1107  fstrsolid%STRESS(3*ic-1) = fstrsolid%STRESS(3*ic-1) + edstress(j,2)
1108  fstrsolid%STRESS(3*ic-0) = fstrsolid%STRESS(3*ic-0) + edstress(j,3)
1109 
1110  if( associated(tnstrain) ) then
1111  tnstrain(3*ic-2) = tnstrain(3*ic-2) + tdstrain(j,1)
1112  tnstrain(3*ic-1) = tnstrain(3*ic-1) + tdstrain(j,2)
1113  tnstrain(3*ic ) = tnstrain(3*ic ) + tdstrain(j,3)
1114  endif
1115  nnumber(ic) = nnumber(ic) + 1
1116  enddo
1117  !--- calculate elemental stress and strain
1118  ! if( ID_area == hecMESH%my_rank ) then
1119  call elementstress_c2( ic_type, fstrsolid%elements(icel)%gausses, estrain, estress )
1120 
1121  fstrsolid%ESTRAIN(3*icel-2) = estrain(1)
1122  fstrsolid%ESTRAIN(3*icel-1) = estrain(2)
1123  fstrsolid%ESTRAIN(3*icel-0) = estrain(3)
1124  fstrsolid%ESTRESS(3*icel-2) = estress(1)
1125  fstrsolid%ESTRESS(3*icel-1) = estress(2)
1126  fstrsolid%ESTRESS(3*icel-0) = estress(3)
1127 
1128  s11 = estress(1)
1129  s22 = estress(2)
1130  s12 = estress(3)
1131  smises = 0.5d0 * ((s11-s22)**2+(s11)**2+(s22)**2) + 3*s12**2
1132  fstrsolid%EMISES(icel) = sqrt( smises )
1133  ! endif
1134  enddo
1135  deallocate( func, inv_func )
1136  enddo
1137 
1138  !C** average over nodes
1139  do i = 1, hecmesh%n_node
1140  if( nnumber(i) == 0 ) cycle
1141  fstrsolid%STRAIN(3*i-2:3*i-0) = fstrsolid%STRAIN(3*i-2:3*i-0) / nnumber(i)
1142  fstrsolid%STRESS(3*i-2:3*i-0) = fstrsolid%STRESS(3*i-2:3*i-0) / nnumber(i)
1143  if( associated(tnstrain) ) tnstrain(3*i-2:3*i) = tnstrain(3*i-2:3*i) / nnumber(i)
1144  enddo
1145  !C** calculate von MISES stress
1146  do i = 1, hecmesh%n_node
1147  s11 = fstrsolid%STRESS(3*i-2)
1148  s22 = fstrsolid%STRESS(3*i-1)
1149  s12 = fstrsolid%STRESS(3*i-0)
1150  smises = 0.5d0 * ((s11-s22)**2+(s11)**2+(s22)**2) + 3*s12**2
1151  fstrsolid%MISES(i) = sqrt( smises )
1152  enddo
1153 
1154  deallocate( nnumber )
1155  end subroutine fstr_nodalstress2d
1156 
1157  !----------------------------------------------------------------------*
1158  subroutine nodalstress_inv2( etype, ni, gausses, func, edstrain, edstress, tdstrain )
1159  !----------------------------------------------------------------------*
1160  use mmechgauss
1161  integer(kind=kint) :: etype, ni
1162  type(tgaussstatus) :: gausses(:)
1163  real(kind=kreal) :: func(:,:), edstrain(:,:), edstress(:,:), tdstrain(:,:)
1164  integer :: i, j, k, ic
1165 
1166  edstrain = 0.0d0
1167  edstress = 0.0d0
1168  tdstrain = 0.0d0
1169 
1170  if( etype == fe_quad4n ) then
1171  do i = 1, ni
1172  do j = 1, ni
1173  do k = 1, 4
1174  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
1175  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
1176  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
1177  enddo
1178  enddo
1179  enddo
1180  else if( etype == fe_tri6n ) then
1181  do i = 1, ni
1182  do j = 1, ni
1183  do k = 1, 4
1184  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
1185  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
1186  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
1187  enddo
1188  enddo
1189  enddo
1190  edstrain(4,1:4) = ( edstrain(1,1:4) + edstrain(2,1:4) ) / 2.0
1191  edstress(4,1:4) = ( edstress(1,1:4) + edstress(2,1:4) ) / 2.0
1192  tdstrain(4,1:4) = ( tdstrain(1,1:4) + tdstrain(2,1:4) ) / 2.0
1193  edstrain(5,1:4) = ( edstrain(2,1:4) + edstrain(3,1:4) ) / 2.0
1194  edstress(5,1:4) = ( edstress(2,1:4) + edstress(3,1:4) ) / 2.0
1195  tdstrain(5,1:4) = ( tdstrain(2,1:4) + tdstrain(3,1:4) ) / 2.0
1196  edstrain(6,1:4) = ( edstrain(3,1:4) + edstrain(1,1:4) ) / 2.0
1197  edstress(6,1:4) = ( edstress(3,1:4) + edstress(1,1:4) ) / 2.0
1198  tdstrain(6,1:4) = ( tdstrain(3,1:4) + tdstrain(1,1:4) ) / 2.0
1199  else if( etype == fe_quad8n ) then
1200  do i = 1, ni
1201  ic = 0
1202  do j = 1, numofquadpoints(etype)
1203  if( j==1 .or. j==3 .or. j==7 .or. j==9 ) then
1204  ic = ic + 1
1205  do k = 1, 4
1206  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
1207  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
1208  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
1209  enddo
1210  endif
1211  enddo
1212  enddo
1213  edstrain(5,1:4) = ( edstrain(1,1:4) + edstrain(2,1:4) ) / 2.0
1214  edstress(5,1:4) = ( edstress(1,1:4) + edstress(2,1:4) ) / 2.0
1215  tdstrain(5,1:4) = ( tdstrain(1,1:4) + tdstrain(2,1:4) ) / 2.0
1216  edstrain(6,1:4) = ( edstrain(2,1:4) + edstrain(3,1:4) ) / 2.0
1217  edstress(6,1:4) = ( edstress(2,1:4) + edstress(3,1:4) ) / 2.0
1218  tdstrain(6,1:4) = ( tdstrain(2,1:4) + tdstrain(3,1:4) ) / 2.0
1219  edstrain(7,1:4) = ( edstrain(3,1:4) + edstrain(4,1:4) ) / 2.0
1220  edstress(7,1:4) = ( edstress(3,1:4) + edstress(4,1:4) ) / 2.0
1221  tdstrain(7,1:4) = ( tdstrain(3,1:4) + tdstrain(4,1:4) ) / 2.0
1222  edstrain(8,1:4) = ( edstrain(4,1:4) + edstrain(1,1:4) ) / 2.0
1223  edstress(8,1:4) = ( edstress(4,1:4) + edstress(1,1:4) ) / 2.0
1224  tdstrain(8,1:4) = ( tdstrain(4,1:4) + tdstrain(1,1:4) ) / 2.0
1225  endif
1226  end subroutine nodalstress_inv2
1227 
1228  !----------------------------------------------------------------------*
1229  subroutine inverse_func( n, a, inv_a )
1230  !----------------------------------------------------------------------*
1231  integer(kind=kint) :: n
1232  real(kind=kreal) :: a(:,:), inv_a(:,:)
1233  integer(kind=kint) :: i, j, k
1234  real(kind=kreal) :: buf
1235 
1236  do i = 1, n
1237  do j = 1, n
1238  if( i == j ) then
1239  inv_a(i,j) = 1.0
1240  else
1241  inv_a(i,j) = 0.0
1242  endif
1243  enddo
1244  enddo
1245 
1246  do i = 1, n
1247  buf = 1.0 / a(i,i)
1248  do j = 1, n
1249  a(i,j) = a(i,j) * buf
1250  inv_a(i,j) = inv_a(i,j) *buf
1251  enddo
1252  do j = 1, n
1253  if( i /= j ) then
1254  buf = a(j,i)
1255  do k = 1, n
1256  a(j,k) = a(j,k) - a(i,k) * buf
1257  inv_a(j,k) = inv_a(j,k) - inv_a(i,k) * buf
1258  enddo
1259  endif
1260  enddo
1261  enddo
1262  end subroutine inverse_func
1263 
1265  !----------------------------------------------------------------------*
1266  subroutine fstr_nodalstress6d( hecMESH, fstrSOLID )
1267  !----------------------------------------------------------------------*
1268  use m_static_lib
1269  use m_fstr_nodalkinematics, only: fstr_ensure_finite_rotation_state
1270  type (hecmwST_local_mesh) :: hecMESH
1271  type (fstr_solid) :: fstrSOLID
1272  !C** local variables
1273  integer(kind=kint) :: itype, icel, is, iE, jS, i, j, k, ic, ic_type, nn, isect, ihead, ID_area, nbase
1274  integer(kind=kint) :: nodLOCAL(20), ntot_lyr, nlyr, n_totlyr
1275  real(kind=kreal) :: ecoord(3,9), edisp(6,9), estrain(6), estress(6), ndstrain(9,6), ndstress(9,6)
1276  real(kind=kreal) :: enqm(12)
1277  real(kind=kreal) :: triad_cur(9,9), triad_ref(9,9)
1278  real(kind=kreal) :: thick, coords(3,3)
1279  integer(kind=kint) :: cdsys_ID
1280  real(kind=kreal) :: mid_strain(9,6), mid_stress(9,6)
1281  integer(kind=kint), allocatable :: nnumber(:)
1282  type(fstr_solid_physic_val), pointer :: layer => null()
1283 
1284  call fstr_solid_phys_clear(fstrsolid)
1285 
1286  n_totlyr = fstrsolid%max_lyr
1287 
1288  allocate( nnumber(hecmesh%n_node) )
1289  if( .not. associated(fstrsolid%is_rot) ) allocate( fstrsolid%is_rot(hecmesh%n_node) )
1290  nnumber = 0
1291  fstrsolid%is_rot = 0
1292  call fstr_ensure_finite_rotation_state(hecmesh, fstrsolid, 6)
1293 
1294  !C +-------------------------------+
1295  !C | according to ELEMENT TYPE |
1296  !C +-------------------------------+
1297  do itype = 1, hecmesh%n_elem_type
1298  is = hecmesh%elem_type_index(itype-1) + 1
1299  ie = hecmesh%elem_type_index(itype )
1300  ic_type = hecmesh%elem_type_item(itype)
1301  if( .not. (hecmw_is_etype_shell(ic_type) .or. ic_type == 611) ) then
1302  ntot_lyr = 0
1303  cycle
1304  end if
1305  if( ic_type == 611 ) ntot_lyr = 0
1306  nn = hecmw_get_max_node( ic_type )
1307  !C** element loop
1308  do icel = is, ie
1309  js = hecmesh%elem_node_index(icel-1)
1310  id_area = hecmesh%elem_ID(icel*2)
1311  do j = 1, nn
1312  nodlocal(j) = hecmesh%elem_node_item(js+j)
1313  ecoord(1,j) = hecmesh%node(3*nodlocal(j)-2)
1314  ecoord(2,j) = hecmesh%node(3*nodlocal(j)-1)
1315  ecoord(3,j) = hecmesh%node(3*nodlocal(j) )
1316  edisp(1,j) = fstrsolid%unode(6*nodlocal(j)-5)
1317  edisp(2,j) = fstrsolid%unode(6*nodlocal(j)-4)
1318  edisp(3,j) = fstrsolid%unode(6*nodlocal(j)-3)
1319  edisp(4,j) = fstrsolid%unode(6*nodlocal(j)-2)
1320  edisp(5,j) = fstrsolid%unode(6*nodlocal(j)-1)
1321  edisp(6,j) = fstrsolid%unode(6*nodlocal(j) )
1322  enddo
1323  isect = hecmesh%section_ID(icel)
1324  ihead = hecmesh%section%sect_R_index(isect-1)
1325  thick = hecmesh%section%sect_R_item(ihead+1)
1326  cdsys_id = hecmesh%section%sect_orien_ID(isect)
1327  coords = 0.0d0
1328  if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords, icel)
1329  !--- calculate elemental stress and strain
1330  if( ic_type == 611 ) then
1331  estrain = 0.0d0
1332  estress = 0.0d0
1333  enqm = 0.0d0
1334  ndstrain = 0.0d0
1335  ndstress = 0.0d0
1336  call nodalstress_beam( nn, ecoord, fstrsolid%elements(icel)%gausses, &
1337  & hecmesh%section%sect_R_item(ihead+1:), edisp(1:6,1:nn), &
1338  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), fstrsolid%sections(isect)%elemopt611 )
1339  call elementalstress_beam( fstrsolid%elements(icel)%gausses, estrain, estress, enqm )
1340  fstrsolid%ENQM(icel*12-11:icel*12) = enqm(1:12)
1341  else if( fstr_uses_plane_stress_j2_shell(ic_type, nn, fstrsolid%elements(icel)%gausses(1)%pMaterial) ) then
1342  ntot_lyr = fstrsolid%elements(icel)%shell_nlayer
1343  call nodalstress_shellj2(fstrsolid%elements(icel), nn, fstrsolid, icel, nodlocal, &
1344  ndstrain(1:nn,1:6), ndstress(1:nn,1:6), estrain, estress)
1345  else if( ic_type == 731 .or. ic_type == 741 .or. ic_type == 743 ) then
1346  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
1347  if( ic_type == 731 .or. ic_type == 741 ) then
1348  do j = 1, nn
1349  nbase = 9*(nodlocal(j)-1)
1350  triad_cur(1:9,j) = 0.0d0
1351  triad_ref(1:9,j) = 0.0d0
1352  if( associated(fstrsolid%shell_triad) ) triad_cur(1:9,j) = fstrsolid%shell_triad(nbase+1:nbase+9)
1353  if( associated(fstrsolid%shell_ref_triad) ) triad_ref(1:9,j) = fstrsolid%shell_ref_triad(nbase+1:nbase+9)
1354  enddo
1355  endif
1356  mid_strain = 0.0d0
1357  mid_stress = 0.0d0
1358  do nlyr=1,ntot_lyr
1359  if( ic_type == 731 .or. ic_type == 741 ) then
1360  call elementstress_shell_mitc( ic_type, nn, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
1361  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), thick, 1.0d0, nlyr, &
1362  & ndtriad=triad_cur(1:9,1:nn), ndreftriad=triad_ref(1:9,1:nn), cdsys_id=cdsys_id, coords=coords)
1363  else
1364  call elementstress_shell_mitc( ic_type, nn, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
1365  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), thick, 1.0d0, nlyr, cdsys_id=cdsys_id, coords=coords)
1366  endif
1367  do j = 1, nn
1368  i = nodlocal(j)
1369  layer => fstrsolid%SHELL%LAYER(nlyr)%PLUS
1370  do k = 1, 6
1371  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + ndstrain(j,k)
1372  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + ndstress(j,k)
1373  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + ndstrain(j,k)/nn
1374  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + ndstress(j,k)/nn
1375  enddo
1376  enddo
1377  call fstr_add_shell_mid(fstrsolid%elements(icel),nlyr,nn,ndstrain,ndstress,mid_strain,mid_stress)
1378  !minus section
1379  if( ic_type == 731 .or. ic_type == 741 ) then
1380  call elementstress_shell_mitc( ic_type, nn, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
1381  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), thick,-1.0d0, nlyr, &
1382  & ndtriad=triad_cur(1:9,1:nn), ndreftriad=triad_ref(1:9,1:nn), cdsys_id=cdsys_id, coords=coords)
1383  else
1384  call elementstress_shell_mitc( ic_type, nn, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
1385  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), thick,-1.0d0, nlyr, cdsys_id=cdsys_id, coords=coords)
1386  endif
1387  do j = 1, nn
1388  i = nodlocal(j)
1389  layer => fstrsolid%SHELL%LAYER(nlyr)%MINUS
1390  do k = 1, 6
1391  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + ndstrain(j,k)
1392  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + ndstress(j,k)
1393  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + ndstrain(j,k)/nn
1394  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + ndstress(j,k)/nn
1395  enddo
1396  enddo
1397  call fstr_add_shell_mid(fstrsolid%elements(icel),nlyr,nn,ndstrain,ndstress,mid_strain,mid_stress)
1398  enddo
1399  ndstrain(1:nn,1:6) = mid_strain(1:nn,1:6)
1400  ndstress(1:nn,1:6) = mid_stress(1:nn,1:6)
1401  call fstr_getavg_shell(fstrsolid,icel,estrain,estress)
1402  endif
1403 
1404  !if( ID_area == hecMESH%my_rank ) then
1405  !ADD VALUE and Count node
1406  do j = 1, nn
1407  ic = hecmesh%elem_node_item(js+j)
1408  fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) + ndstrain(j,1:6)
1409  fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) + ndstress(j,1:6)
1410  !if( associated(tnstrain) )then
1411  ! tnstrain(6*(ic-1)+1:6*(ic-1)+6) = tnstrain(6*(ic-1)+1:6*(ic-1)+6) + tdstrain(j,1:6)
1412  !endif
1413  nnumber(ic) = nnumber(ic) + 1
1414  enddo
1415 
1416  fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) + estrain(1:6)
1417  fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) + estress(1:6)
1418  !endif
1419  enddo
1420  enddo
1421 
1422  !C** calculate nodal stress and strain
1423  do i = 1, hecmesh%n_node
1424  if( nnumber(i) == 0 ) cycle
1425  fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1426  fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1427  !if( associated(tnstrain) )then
1428  ! tnstrain(6*(i-1)+1:6*(i-1)+6) = tnstrain(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1429  !endif
1430  enddo
1431 
1432  do nlyr = 1, ntot_lyr
1433  do i = 1, hecmesh%n_node
1434  if( nnumber(i) == 0 ) cycle
1435  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
1436  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1437  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
1438  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1439  fstrsolid%SHELL%LAYER(nlyr)%PLUS%MISES(i) = &
1440  & get_mises(fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6))
1441 
1442  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
1443  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1444  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
1445  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1446  fstrsolid%SHELL%LAYER(nlyr)%MINUS%MISES(i) = &
1447  & get_mises(fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6))
1448  enddo
1449  enddo
1450 
1451  !C** calculate von MISES stress
1452  do i = 1, hecmesh%n_node
1453  fstrsolid%MISES(i) = get_mises(fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6))
1454  enddo
1455  do i = 1, hecmesh%n_elem
1456  fstrsolid%EMISES(i) = get_mises(fstrsolid%ESTRESS(6*(i-1)+1:6*(i-1)+6))
1457  enddo
1458 
1459  !C** calculate Elemental Plastic Strain
1460  do i = 1, hecmesh%n_elem
1461  fstrsolid%EPLSTRAIN(i) = fstr_element_average_plstrain(fstrsolid%elements(i))
1462  enddo
1463 
1464  deallocate( nnumber )
1465 
1466  end subroutine fstr_nodalstress6d
1467 
1468  subroutine make_principal(fstrSOLID, hecMESH, RES)
1469  use hecmw_util
1470  use m_out
1471  use m_static_lib
1472 
1473  type(fstr_solid) :: fstrSOLID
1474  type(hecmwst_local_mesh) :: hecMESH
1475  type(fstr_solid_physic_val) :: RES
1476  integer(kind=kint) :: i, flag
1477  real(kind=kreal) :: tmat(3, 3), tvec(3), strain(6)
1478 
1479  flag=0
1480  if( fstrsolid%output_ctrl(3)%outinfo%on(19) .or. fstrsolid%output_ctrl(4)%outinfo%on(19) ) then
1481  if ( .not. associated(res%PSTRESS) ) then
1482  allocate(res%PSTRESS( 3*hecmesh%n_node ))
1483  endif
1484  flag=ior(flag,b'00000001')
1485  end if
1486  if( fstrsolid%output_ctrl(3)%outinfo%on(23) .or. fstrsolid%output_ctrl(4)%outinfo%on(23) ) then
1487  if ( .not. associated(res%PSTRESS_VECT) ) then
1488  allocate(res%PSTRESS_VECT( 3*hecmesh%n_node ,3))
1489  endif
1490  flag=ior(flag,b'00000010')
1491  end if
1492  if( fstrsolid%output_ctrl(3)%outinfo%on(21) .or. fstrsolid%output_ctrl(4)%outinfo%on(21) ) then
1493  if ( .not. associated(res%PSTRAIN) ) then
1494  allocate(res%PSTRAIN( 3*hecmesh%n_node ))
1495  endif
1496  flag=ior(flag,b'00000100')
1497  end if
1498  if( fstrsolid%output_ctrl(3)%outinfo%on(25) .or. fstrsolid%output_ctrl(4)%outinfo%on(25) ) then
1499  if ( .not. associated(res%PSTRAIN_VECT) ) then
1500  allocate(res%PSTRAIN_VECT( 3*hecmesh%n_node ,3))
1501  endif
1502  flag=ior(flag,b'00001000')
1503  end if
1504  if( fstrsolid%output_ctrl(3)%outinfo%on(20) .or. fstrsolid%output_ctrl(4)%outinfo%on(20) ) then
1505  if ( .not. associated(res%EPSTRESS) ) then
1506  allocate(res%EPSTRESS( 3*hecmesh%n_elem ))
1507  endif
1508  flag=ior(flag,b'00010000')
1509  end if
1510  if( fstrsolid%output_ctrl(3)%outinfo%on(24) .or. fstrsolid%output_ctrl(4)%outinfo%on(24) ) then
1511  if ( .not. associated(res%EPSTRESS_VECT) ) then
1512  allocate(res%EPSTRESS_VECT( 3*hecmesh%n_elem ,3))
1513  endif
1514  flag=ior(flag,b'00100000')
1515  end if
1516  if( fstrsolid%output_ctrl(3)%outinfo%on(22) .or. fstrsolid%output_ctrl(4)%outinfo%on(22) ) then
1517  if ( .not. associated(res%EPSTRAIN) ) then
1518  allocate(res%EPSTRAIN( 3*hecmesh%n_elem ))
1519  endif
1520  flag=ior(flag,b'01000000')
1521  end if
1522  if( fstrsolid%output_ctrl(3)%outinfo%on(26) .or. fstrsolid%output_ctrl(4)%outinfo%on(26) ) then
1523  if ( .not. associated(res%EPSTRAIN_VECT) ) then
1524  allocate(res%EPSTRAIN_VECT( 3*hecmesh%n_elem ,3))
1525  endif
1526  flag=ior(flag,b'10000000')
1527  end if
1528 
1529  if (iand(flag,b'00000011') /= 0) then
1530  do i = 1, hecmesh%n_node
1531  call get_principal(res%STRESS(6*i-5:6*i), tvec, tmat)
1532  if (iand(flag,b'00000001') /= 0) res%PSTRESS(3*(i-1)+1:3*(i-1)+3)=tvec
1533  if (iand(flag,b'00000010') /= 0) res%PSTRESS_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
1534  end do
1535  end if
1536  if (iand(flag,b'00001100') /= 0) then
1537  do i = 1, hecmesh%n_node
1538  strain(1:6) = res%STRAIN(6*i-5:6*i)
1539  strain(4:6) = 0.5d0*strain(4:6)
1540  call get_principal(strain, tvec, tmat)
1541  if (iand(flag,b'00000100') /= 0) res%PSTRAIN(3*(i-1)+1:3*(i-1)+3)=tvec
1542  if (iand(flag,b'00001000') /= 0) res%PSTRAIN_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
1543  end do
1544  end if
1545 
1546  if (iand(flag,b'00110000') /= 0) then
1547  do i = 1, hecmesh%n_elem
1548  call get_principal( res%ESTRESS(6*i-5:6*i), tvec, tmat)
1549  if (iand(flag,b'00010000') /= 0) res%EPSTRESS(3*(i-1)+1:3*(i-1)+3)=tvec
1550  if (iand(flag,b'00100000') /= 0) res%EPSTRESS_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
1551  end do
1552  end if
1553  if (iand(flag,b'11000000') /= 0) then
1554  do i = 1, hecmesh%n_elem
1555  strain(1:6) = res%ESTRAIN(6*i-5:6*i)
1556  strain(4:6) = 0.5d0*strain(4:6)
1557  call get_principal(strain, tvec, tmat)
1558  if (iand(flag,b'01000000') /= 0) res%EPSTRAIN(3*(i-1)+1:3*(i-1)+3)=tvec
1559  if (iand(flag,b'10000000') /= 0) res%EPSTRAIN_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
1560  end do
1561  end if
1562  end subroutine make_principal
1563 
1564 end module m_fstr_nodalstress
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:754
elementinfo::getquadpoint
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:539
elementinfo::getshellthicknessquadpoint
subroutine getshellthicknessquadpoint(etype, np, zeta)
Fetch the through-thickness quadrature coordinate of a shell element.
Definition: element.f90:504
elementinfo::numofshellthicknessquadpoints
integer function numofshellthicknessquadpoints(etype)
Obtains the number of through-thickness quadrature points of a shell element.
Definition: element.f90:490
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_fstr_nodalkinematics
Finite-rotation nodal kinematics for NLGEOM.
Definition: fstr_NodalKinematics.f90:14
m_fstr_nodalkinematics::fstr_ensure_finite_rotation_state
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...
Definition: fstr_NodalKinematics.f90:34
m_fstr_nodalstress
This module provides functions to calculation nodal stress.
Definition: fstr_NodalStress.f90:6
m_fstr_nodalstress::get_shell_layer_gauss_average
subroutine get_shell_layer_gauss_average(element, estrain, estress)
Definition: fstr_NodalStress.f90:756
m_fstr_nodalstress::set_shell_layer_surface_results
subroutine set_shell_layer_surface_results(element, fstrSOLID, icel)
Definition: fstr_NodalStress.f90:678
m_fstr_nodalstress::fstr_stress_add_shelllyr
subroutine fstr_stress_add_shelllyr(nn, fstrSOLID, icel, nodLOCAL, nlyr, strain, stress, flag)
Definition: fstr_NodalStress.f90:609
m_fstr_nodalstress::nodalstress_shellj2
subroutine nodalstress_shellj2(element, nn, fstrSOLID, icel, nodLOCAL, strain, stress, estrain, estress)
Recover nodal results from the converged shell material histories.
Definition: fstr_NodalStress.f90:799
m_fstr_nodalstress::fstr_nodalstress3d
subroutine fstr_nodalstress3d(hecMESH, fstrSOLID)
Calculate NODAL STRESS of solid elements.
Definition: fstr_NodalStress.f90:61
m_fstr_nodalstress::fstr_nodalstress6d
subroutine fstr_nodalstress6d(hecMESH, fstrSOLID)
Calculate NODAL STRESS of shell elements.
Definition: fstr_NodalStress.f90:1267
m_fstr_nodalstress::get_mises
real(kind=kreal) function get_mises(s)
Definition: fstr_NodalStress.f90:1005
m_fstr_nodalstress::fstr_nodalstress3d_c3d4_sesns
subroutine fstr_nodalstress3d_c3d4_sesns(hecMESH, fstrSOLID, nnumber, Nodal_STRAIN, Nodal_STRESS, Elemental_STRAIN, Elemental_STRESS)
Definition: fstr_NodalStress.f90:422
m_fstr_nodalstress::fstr_nodalstress2d
subroutine fstr_nodalstress2d(hecMESH, fstrSOLID)
Calculate NODAL STRESS of plane elements.
Definition: fstr_NodalStress.f90:1024
m_fstr_nodalstress::fstr_add_shell_mid
subroutine fstr_add_shell_mid(element, nlyr, nn, strain, stress, mid_strain, mid_stress)
Adds one layer surface's nodal values of an element, with half the layer's weight,...
Definition: fstr_NodalStress.f90:638
m_fstr_nodalstress::fstr_getavg_shell
subroutine fstr_getavg_shell(fstrSOLID, icel, estrain, estress)
The element's mid-surface values. The nodal ones come from fstr_add_shell_mid: the layers' nodal arra...
Definition: fstr_NodalStress.f90:653
m_fstr_nodalstress::get_shell_layer_surface_average
subroutine get_shell_layer_surface_average(element, ilayer, flag, estrain, estress, eplstrain)
Definition: fstr_NodalStress.f90:701
m_fstr_nodalstress::fstr_get_shell_gauss_output
subroutine fstr_get_shell_gauss_output(hecMESH, fstrSOLID, icel, k, strain, stress)
Read stored shell tensors or evaluate one elastic thickness point without changing history.
Definition: fstr_NodalStress.f90:16
m_fstr_nodalstress::make_principal
subroutine make_principal(fstrSOLID, hecMESH, RES)
Definition: fstr_NodalStress.f90:1469
m_fstr_nodalstress::search_idx_senes
integer(kind=kint) function search_idx_senes(irow, asect, nid, sid)
Definition: fstr_NodalStress.f90:403
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.F90:15
m_fstr::kel341sesns
integer(kind=kint), parameter kel341sesns
Definition: m_fstr.F90:82
m_fstr::get_coordsys
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords, icel)
This subroutine fetch coords defined by local coordinate system.
Definition: m_fstr.F90:1307
m_fstr::fstr_solid_phys_clear
subroutine fstr_solid_phys_clear(fstrSOLID)
Definition: m_fstr.F90:1386
m_out
This module manages step information.
Definition: m_out.f90:6
m_static_lib
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
mmechgauss::fstr_shell_output_point_count
integer(kind=kint) function fstr_shell_output_point_count(element)
Number of shell thickness-point tensors available for result output.
Definition: mechgauss.f90:52
mmechgauss::fstr_element_average_plstrain
real(kind=kreal) function fstr_element_average_plstrain(element)
Average equivalent plastic strain over the active quadrature history. Average shell history points eq...
Definition: mechgauss.f90:194
hecmw_util::hecmwst_local_mesh
Definition: hecmw_util_f.F90:246
m_fstr::fstr_solid_physic_val
Data for STATIC ANSLYSIS (fstrSOLID)
Definition: m_fstr.F90:203
m_fstr::fstr_solid
Definition: m_fstr.F90:230
mmechgauss::telement
All data should be recorded in every elements.
Definition: mechgauss.f90:34
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:15