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