FrontISTR  5.8.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, ntot_lyr)
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, ntot_lyr)
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, ntot_lyr)
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, ntot_lyr)
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  if( flag33 == 1 )then
247  do nlyr = 1, ntot_lyr
248  do i = 1, hecmesh%n_node
249  if( nnumber(i) == 0 ) cycle
250  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
251  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
252  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
253  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
254  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
255  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
256  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
257  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
258  enddo
259  enddo
260  endif
261 
262  !C** calculate von MISES stress
263  do i = 1, hecmesh%n_node
264  fstrsolid%MISES(i) = get_mises(fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6))
265  enddo
266  do i = 1, hecmesh%n_elem
267  fstrsolid%EMISES(i) = get_mises(fstrsolid%ESTRESS(6*(i-1)+1:6*(i-1)+6))
268  enddo
269 
270  !C** calculate Elemental Plastic Strain
271  do i = 1, hecmesh%n_elem
272  if (.not. associated(fstrsolid%elements(i)%gausses)) cycle
273  fstrsolid%EPLSTRAIN(i) = get_pl_estrain(fstrsolid%elements(i)%gausses)
274  enddo
275 
276  if( flag33 == 1 )then
277  if( fstrsolid%output_ctrl(3)%outinfo%on(27) .or. fstrsolid%output_ctrl(4)%outinfo%on(27) ) then
278  do nlyr = 1, ntot_lyr
279  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL%LAYER(nlyr)%PLUS)
280  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL%LAYER(nlyr)%MINUS)
281  enddo
282  endif
283  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL)
284  else
285  call make_principal(fstrsolid, hecmesh, fstrsolid%SOLID)
286  endif
287 
288  deallocate( nnumber )
289 
290  end subroutine fstr_nodalstress3d
291 
292  integer(kind=kint) function search_idx_senes( irow, asect, nid, sid )
293  integer(kind=kint), allocatable, intent(in) :: irow(:)
294  integer(kind=kint), allocatable, intent(in) :: asect(:)
295  integer(kind=kint), intent(in) :: nid
296  integer(kind=kint), intent(in) :: sid
297 
298  integer(kind=kint) :: i
299 
300  search_idx_senes = -1
301  do i=irow(nid-1)+1,irow(nid)
302  if( asect(i) == sid ) then
303  search_idx_senes = i
304  return
305  end if
306  end do
307 
308  end function
309 
310  subroutine fstr_nodalstress3d_c3d4_sesns( hecMESH, fstrSOLID, nnumber, &
311  Nodal_STRAIN, Nodal_STRESS, Elemental_STRAIN, Elemental_STRESS )
312  type(hecmwst_local_mesh),intent(in) :: hecMESH
313  type(fstr_solid),intent(inout) :: fstrSOLID
314  integer(kind=kint), allocatable, intent(inout) :: nnumber(:)
315  real(kind=kreal), pointer, intent(inout) :: nodal_strain(:)
316  real(kind=kreal), pointer, intent(inout) :: nodal_stress(:)
317  real(kind=kreal), pointer, intent(inout) :: elemental_strain(:)
318  real(kind=kreal), pointer, intent(inout) :: elemental_stress(:)
319 
320  integer(kind=kint) :: itype, iS, iE, jS, ic_type, icel, i, j, isect
321  integer(kind=kint) :: nsize, nid(2), idx(2), nd
322  integer(kind=kint) :: nnode, nlen
323  type(hecmwst_varray_int), allocatable :: nodal_sections(:)
324  real(kind=kreal) :: tmpval(6), hydval, nsecdup
325  integer(kind=kint), allocatable :: irow(:), jcol(:), asect(:)
326  real(kind=kreal), allocatable :: stress_hyd(:), strain_hyd(:)
327  real(kind=kreal), allocatable :: stress_dev(:)
328  real(kind=kreal) :: stress_hyd_ndave(6), strain_hyd_ndave(6)
329  real(kind=kreal) :: stress_dev_ndave(6), strain_dev_ndave(6)
330  real(kind=kreal), allocatable :: n_dup_dev(:), n_dup_hyd(:)
331  real(kind=kreal) :: edstrain(6), edstress(6)
332 
333  nnode = hecmesh%n_node
334  nsize = size(nodal_strain)
335 
336  ! create section info at node
337  call hecmw_varray_int_initialize_all( nodal_sections, nnode, 2 )
338  do itype = 1, hecmesh%n_elem_type
339  ic_type = hecmesh%elem_type_item(itype)
340  if( ic_type /= 341 ) cycle
341 
342  is = hecmesh%elem_type_index(itype-1) + 1
343  ie = hecmesh%elem_type_index(itype )
344 
345  do icel=is,ie
346  isect= hecmesh%section_ID(icel)
347  if( fstrsolid%sections(isect)%elemopt341 /= kel341sesns ) cycle
348  js = hecmesh%elem_node_index(icel-1)
349  do i=1,4
350  nd = hecmesh%elem_node_item(js+i)
351  call hecmw_varray_int_add_if_not_exits( nodal_sections(nd), isect )
352  end do
353  end do
354  enddo
355 
356  ! create CRS arrays of nodal stress/strain with different sections
357  allocate(irow(0:nnode))
358  irow(0) = 0
359  do i=1,nnode
360  irow(i) = irow(i-1)+hecmw_varray_int_get_nitem(nodal_sections(i))
361  end do
362  nlen = irow(nnode)
363 
364  allocate(asect(nlen))
365  do i=1,nnode
366  if( irow(i-1) == irow(i) ) cycle
367  call hecmw_varray_int_get_item_all( nodal_sections(i), asect(irow(i-1)+1:irow(i)) )
368  end do
369 
370  ! add stress/strain from smoothed elements
371  allocate(stress_hyd(6*nlen), strain_hyd(6*nlen))
372  allocate(stress_dev(6*nlen))
373  allocate(n_dup_dev(nlen),n_dup_hyd(nlen))
374 
375  stress_hyd(:) = 0.d0
376  strain_hyd(:) = 0.d0
377  stress_dev(:) = 0.d0
378  n_dup_hyd(:) = 0.d0
379  n_dup_dev(:) = 0.d0
380  do itype = 1, hecmesh%n_elem_type
381  ic_type = hecmesh%elem_type_item(itype)
382  if( ic_type /= 881 .and. ic_type /= 891 ) cycle
383 
384  is = hecmesh%elem_type_index(itype-1) + 1
385  ie = hecmesh%elem_type_index(itype )
386 
387  do icel=is,ie
388  js = hecmesh%elem_node_index(icel-1)
389  isect= hecmesh%section_ID(icel)
390  if( ic_type == 881 ) then
391  nid(1) = hecmesh%elem_node_item(js+1)
392  idx(1) = search_idx_senes( irow, asect, nid(1), isect )
393 
394  !strain
395  strain_hyd(6*idx(1)-5:6*idx(1)) = fstrsolid%elements(icel)%gausses(1)%strain_out(1:6)
396  !stress
397  stress_hyd(6*idx(1)-5:6*idx(1)) = fstrsolid%elements(icel)%gausses(1)%stress_out(1:6)
398  !number of duplication
399  n_dup_hyd(idx(1)) = n_dup_hyd(idx(1)) + 1.d0
400  else if( ic_type == 891 ) then
401  nid(1:2) = hecmesh%elem_node_item(js+1:js+2)
402  idx(1) = search_idx_senes( irow, asect, nid(1), isect )
403  idx(2) = search_idx_senes( irow, asect, nid(2), isect )
404 
405  !stress
406  tmpval(1:6) = fstrsolid%elements(icel)%gausses(1)%stress_out(1:6)
407  stress_dev(6*idx(1)-5:6*idx(1)) = stress_dev(6*idx(1)-5:6*idx(1)) + tmpval(1:6)
408  stress_dev(6*idx(2)-5:6*idx(2)) = stress_dev(6*idx(2)-5:6*idx(2)) + tmpval(1:6)
409  !number of duplication
410  n_dup_dev(idx(1)) = n_dup_dev(idx(1)) + 1.d0
411  n_dup_dev(idx(2)) = n_dup_dev(idx(2)) + 1.d0
412  end if
413  end do
414  enddo
415 
416  do i=1,nnode
417  if( irow(i-1) == irow(i) ) cycle
418  do j=irow(i-1)+1,irow(i)
419  if( n_dup_dev(j) < 1.0d-8 ) cycle
420  stress_dev(6*j-5:6*j) = stress_dev(6*j-5:6*j)/n_dup_dev(j)
421  end do
422  end do
423 
424  ! average at node for nodal output
425  do i=1,nnode
426  if( irow(i-1) == irow(i) ) cycle
427  strain_hyd_ndave(:) = 0.d0
428  stress_hyd_ndave(:) = 0.d0
429  stress_dev_ndave(:) = 0.d0
430  do j=irow(i-1)+1,irow(i)
431  strain_hyd_ndave(1:6) = strain_hyd_ndave(1:6) + strain_hyd(6*j-5:6*j)
432  stress_hyd_ndave(1:6) = stress_hyd_ndave(1:6) + stress_hyd(6*j-5:6*j)
433  stress_dev_ndave(1:6) = stress_dev_ndave(1:6) + stress_dev(6*j-5:6*j)
434  end do
435  nsecdup = dble(irow(i)-irow(i-1))
436  strain_hyd_ndave(1:6) = strain_hyd_ndave(1:6)/nsecdup
437  stress_hyd_ndave(1:6) = stress_hyd_ndave(1:6)/nsecdup
438  stress_dev_ndave(1:6) = stress_dev_ndave(1:6)/nsecdup
439 
440  if( nnumber(i) == 0 ) then
441  nodal_strain(6*i-5:6*i) = strain_hyd_ndave(1:6)
442  nodal_stress(6*i-5:6*i) = stress_hyd_ndave(1:6)+stress_dev_ndave(1:6)
443  else
444  nodal_strain(6*i-5:6*i) = 0.5d0*(nodal_strain(6*i-5:6*i)+strain_hyd_ndave(1:6))
445  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))
446  endif
447  end do
448 
449  ! ELEMENTAL STRAIN and STRESS
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  edstrain(1:6) = 0.d0
462  edstress(1:6) = 0.d0
463  do i=1,4
464  nd = hecmesh%elem_node_item(js+i)
465  idx(1) = search_idx_senes( irow, asect, hecmesh%elem_node_item(js+i), isect )
466  edstrain(1:6) = edstrain(1:6) + strain_hyd(6*idx(1)-5:6*idx(1))
467  edstress(1:6) = edstress(1:6) + stress_hyd(6*idx(1)-5:6*idx(1)) + stress_dev(6*idx(1)-5:6*idx(1))
468  end do
469  edstrain(1:6) = 0.25d0*edstrain(1:6)
470  edstress(1:6) = 0.25d0*edstress(1:6)
471 
472  elemental_strain(6*(icel-1)+1:6*(icel-1)+6) = elemental_strain(6*(icel-1)+1:6*(icel-1)+6) + edstrain(1:6)
473  elemental_stress(6*(icel-1)+1:6*(icel-1)+6) = elemental_stress(6*(icel-1)+1:6*(icel-1)+6) + edstress(1:6)
474 
475  fstrsolid%elements(icel)%gausses(1)%strain_out(1:6) = elemental_strain(6*(icel-1)+1:6*(icel-1)+6)
476  fstrsolid%elements(icel)%gausses(1)%stress_out(1:6) = elemental_stress(6*(icel-1)+1:6*(icel-1)+6)
477  end do
478  enddo
479 
480  deallocate(stress_hyd, strain_hyd)
481  deallocate(stress_dev)
482  deallocate(n_dup_dev, n_dup_hyd)
483 
484  end subroutine
485 
486  subroutine fstr_stress_add_shelllyr(nn,fstrSOLID,icel,nodLOCAL,nlyr,strain,stress,flag)
487  implicit none
488  type(fstr_solid) :: fstrsolid
489  integer(kind=kint) :: nodlocal(20)
490  integer(kind=kint) :: nn, i, j, k, m, nlyr, weight, icel, flag
491  real(kind=kreal) :: strain(nn, 6), stress(nn, 6)
492  type(fstr_solid_physic_val), pointer :: layer => null()
493 
494  do j = 1, nn
495  i = nodlocal(j)
496  m = nodlocal(j+nn)
497  if(flag == 1)then
498  layer => fstrsolid%SHELL%LAYER(nlyr)%PLUS
499  elseif(flag == -1)then
500  layer => fstrsolid%SHELL%LAYER(nlyr)%MINUS
501  endif
502  do k = 1, 6
503  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + strain(j,k)
504  layer%STRAIN(6*(m-1)+k) = layer%STRAIN(6*(m-1)+k) + strain(j,k)
505  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + stress(j,k)
506  layer%STRESS(6*(m-1)+k) = layer%STRESS(6*(m-1)+k) + stress(j,k)
507  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + strain(j,k)/nn
508  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + stress(j,k)/nn
509  enddo
510  enddo
511  end subroutine fstr_stress_add_shelllyr
512 
513  subroutine fstr_getavg_shell(nn,fstrSOLID,icel,nodLOCAL,strain,stress,estrain,estress)
514  implicit none
515  type (fstr_solid) :: fstrsolid
516  integer(kind=kint) :: nodlocal(20)
517  integer(kind=kint) :: nn, i, j, k, m, nlyr, icel, flag, ntot_lyr
518  real(kind=kreal) :: strain(nn,6), stress(nn,6), estrain(6), estress(6), weight
519  type(fstr_solid_physic_val), pointer :: layer => null()
520 
521  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
522  strain = 0.0d0
523  stress = 0.0d0
524  estrain = 0.0d0
525  estress = 0.0d0
526 
527  do nlyr = 1, ntot_lyr
528  layer => fstrsolid%SHELL%LAYER(nlyr)
529  weight = fstrsolid%elements(icel)%gausses(1)%pMaterial%shell_var(nlyr)%weight
530  do j = 1, nn
531  i = nodlocal(j)
532  do k = 1, 6
533  strain(j,k) = strain(j,k) &
534  & + weight*(0.5d0*layer%PLUS%STRAIN(6*(i-1)+k) + 0.5d0*layer%MINUS%STRAIN(6*(i-1)+k))
535  stress(j,k) = stress(j,k) &
536  & + weight*(0.5d0*layer%PLUS%STRESS(6*(i-1)+k) + 0.5d0*layer%MINUS%STRESS(6*(i-1)+k))
537  enddo
538  estrain(j) = estrain(j) &
539  & + weight*(0.5d0*layer%PLUS%ESTRAIN(6*(icel-1)+j) + 0.5d0*layer%MINUS%ESTRAIN(6*(icel-1)+j))
540  estress(j) = estress(j) &
541  & + weight*(0.5d0*layer%PLUS%ESTRESS(6*(icel-1)+j) + 0.5d0*layer%MINUS%ESTRESS(6*(icel-1)+j))
542  enddo
543  enddo
544  end subroutine fstr_getavg_shell
545 
546  !----------------------------------------------------------------------*
547  subroutine nodalstress_inv3( etype, ni, gausses, func, edstrain, edstress, tdstrain )
548  !----------------------------------------------------------------------*
549  use mmechgauss
550  integer(kind=kint) :: etype, ni
551  type(tgaussstatus) :: gausses(:)
552  real(kind=kreal) :: func(:, :), edstrain(:, :), edstress(:, :), tdstrain(:, :)
553  integer :: i, j, k, ic
554 
555  edstrain = 0.0d0
556  edstress = 0.0d0
557  tdstrain = 0.0d0
558 
559  if( etype == fe_hex8n ) then
560  do i = 1, ni
561  do j = 1, ni
562  do k = 1, 6
563  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
564  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
565  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
566  enddo
567  enddo
568  enddo
569  else if( etype == fe_tet10n ) then
570  do i = 1, ni
571  do j = 1, ni
572  do k = 1, 6
573  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
574  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
575  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
576  enddo
577  enddo
578  enddo
579  edstrain(5,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
580  edstress(5,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
581  tdstrain(5,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
582  edstrain(6,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
583  edstress(6,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
584  tdstrain(6,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
585  edstrain(7,1:6) = ( edstrain(3,1:6) + edstrain(1,1:6) ) / 2.0
586  edstress(7,1:6) = ( edstress(3,1:6) + edstress(1,1:6) ) / 2.0
587  tdstrain(7,1:6) = ( tdstrain(3,1:6) + tdstrain(1,1:6) ) / 2.0
588  edstrain(8,1:6) = ( edstrain(1,1:6) + edstrain(4,1:6) ) / 2.0
589  edstress(8,1:6) = ( edstress(1,1:6) + edstress(4,1:6) ) / 2.0
590  tdstrain(8,1:6) = ( tdstrain(1,1:6) + tdstrain(4,1:6) ) / 2.0
591  edstrain(9,1:6) = ( edstrain(2,1:6) + edstrain(4,1:6) ) / 2.0
592  edstress(9,1:6) = ( edstress(2,1:6) + edstress(4,1:6) ) / 2.0
593  tdstrain(9,1:6) = ( tdstrain(2,1:6) + tdstrain(4,1:6) ) / 2.0
594  edstrain(10,1:6) = ( edstrain(3,1:6) + edstrain(4,1:6) ) / 2.0
595  edstress(10,1:6) = ( edstress(3,1:6) + edstress(4,1:6) ) / 2.0
596  tdstrain(10,1:6) = ( tdstrain(3,1:6) + tdstrain(4,1:6) ) / 2.0
597  else if( etype == fe_prism15n ) then
598  do i = 1, ni
599  ic = 0
600  do j = 1, numofquadpoints(etype)
601  if( j==1 .or. j==2 .or. j==3 .or. j==7 .or. j==8 .or. j==9 ) then
602  ic = ic + 1
603  do k = 1, 6
604  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
605  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
606  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
607  enddo
608  endif
609  enddo
610  enddo
611  edstrain(7,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
612  edstress(7,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
613  tdstrain(7,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
614  edstrain(8,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
615  edstress(8,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
616  tdstrain(8,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
617  edstrain(9,1:6) = ( edstrain(3,1:6) + edstrain(1,1:6) ) / 2.0
618  edstress(9,1:6) = ( edstress(3,1:6) + edstress(1,1:6) ) / 2.0
619  tdstrain(9,1:6) = ( tdstrain(3,1:6) + tdstrain(1,1:6) ) / 2.0
620  edstrain(10,1:6) = ( edstrain(4,1:6) + edstrain(5,1:6) ) / 2.0
621  edstress(10,1:6) = ( edstress(4,1:6) + edstress(5,1:6) ) / 2.0
622  tdstrain(10,1:6) = ( tdstrain(4,1:6) + tdstrain(5,1:6) ) / 2.0
623  edstrain(11,1:6) = ( edstrain(5,1:6) + edstrain(6,1:6) ) / 2.0
624  edstress(11,1:6) = ( edstress(5,1:6) + edstress(6,1:6) ) / 2.0
625  tdstrain(11,1:6) = ( tdstrain(5,1:6) + tdstrain(6,1:6) ) / 2.0
626  edstrain(12,1:6) = ( edstrain(6,1:6) + edstrain(4,1:6) ) / 2.0
627  edstress(12,1:6) = ( edstress(6,1:6) + edstress(4,1:6) ) / 2.0
628  tdstrain(12,1:6) = ( tdstrain(6,1:6) + tdstrain(4,1:6) ) / 2.0
629  edstrain(13,1:6) = ( edstrain(1,1:6) + edstrain(4,1:6) ) / 2.0
630  edstress(13,1:6) = ( edstress(1,1:6) + edstress(4,1:6) ) / 2.0
631  tdstrain(13,1:6) = ( tdstrain(1,1:6) + tdstrain(4,1:6) ) / 2.0
632  edstrain(14,1:6) = ( edstrain(2,1:6) + edstrain(5,1:6) ) / 2.0
633  edstress(14,1:6) = ( edstress(2,1:6) + edstress(5,1:6) ) / 2.0
634  tdstrain(14,1:6) = ( tdstrain(2,1:6) + tdstrain(5,1:6) ) / 2.0
635  edstrain(15,1:6) = ( edstrain(3,1:6) + edstrain(6,1:6) ) / 2.0
636  edstress(15,1:6) = ( edstress(3,1:6) + edstress(6,1:6) ) / 2.0
637  tdstrain(15,1:6) = ( tdstrain(3,1:6) + tdstrain(6,1:6) ) / 2.0
638  else if( etype == fe_hex20n ) then
639  do i = 1, ni
640  ic = 0
641  do j = 1, numofquadpoints(etype)
642  if( j==1 .or. j==3 .or. j==7 .or. j==9 .or. &
643  j==19 .or. j==21 .or. j==25 .or. j==27 ) then
644  ic = ic + 1
645  do k = 1, 6
646  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
647  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
648  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
649  enddo
650  endif
651  enddo
652  enddo
653  edstrain(9,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
654  edstress(9,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
655  tdstrain(9,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
656  edstrain(10,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
657  edstress(10,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
658  tdstrain(10,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
659  edstrain(11,1:6) = ( edstrain(3,1:6) + edstrain(4,1:6) ) / 2.0
660  edstress(11,1:6) = ( edstress(3,1:6) + edstress(4,1:6) ) / 2.0
661  tdstrain(11,1:6) = ( tdstrain(3,1:6) + tdstrain(4,1:6) ) / 2.0
662  edstrain(12,1:6) = ( edstrain(4,1:6) + edstrain(1,1:6) ) / 2.0
663  edstress(12,1:6) = ( edstress(4,1:6) + edstress(1,1:6) ) / 2.0
664  tdstrain(12,1:6) = ( tdstrain(4,1:6) + tdstrain(1,1:6) ) / 2.0
665  edstrain(13,1:6) = ( edstrain(5,1:6) + edstrain(6,1:6) ) / 2.0
666  edstress(13,1:6) = ( edstress(5,1:6) + edstress(6,1:6) ) / 2.0
667  tdstrain(13,1:6) = ( tdstrain(5,1:6) + tdstrain(6,1:6) ) / 2.0
668  edstrain(14,1:6) = ( edstrain(6,1:6) + edstrain(7,1:6) ) / 2.0
669  edstress(14,1:6) = ( edstress(6,1:6) + edstress(7,1:6) ) / 2.0
670  tdstrain(14,1:6) = ( tdstrain(6,1:6) + tdstrain(7,1:6) ) / 2.0
671  edstrain(15,1:6) = ( edstrain(7,1:6) + edstrain(8,1:6) ) / 2.0
672  edstress(15,1:6) = ( edstress(7,1:6) + edstress(8,1:6) ) / 2.0
673  tdstrain(15,1:6) = ( tdstrain(7,1:6) + tdstrain(8,1:6) ) / 2.0
674  edstrain(16,1:6) = ( edstrain(8,1:6) + edstrain(5,1:6) ) / 2.0
675  edstress(16,1:6) = ( edstress(8,1:6) + edstress(5,1:6) ) / 2.0
676  tdstrain(16,1:6) = ( tdstrain(8,1:6) + tdstrain(5,1:6) ) / 2.0
677  edstrain(17,1:6) = ( edstrain(1,1:6) + edstrain(5,1:6) ) / 2.0
678  edstress(17,1:6) = ( edstress(1,1:6) + edstress(5,1:6) ) / 2.0
679  tdstrain(17,1:6) = ( tdstrain(1,1:6) + tdstrain(5,1:6) ) / 2.0
680  edstrain(18,1:6) = ( edstrain(2,1:6) + edstrain(6,1:6) ) / 2.0
681  edstress(18,1:6) = ( edstress(2,1:6) + edstress(6,1:6) ) / 2.0
682  tdstrain(18,1:6) = ( tdstrain(2,1:6) + tdstrain(6,1:6) ) / 2.0
683  edstrain(19,1:6) = ( edstrain(3,1:6) + edstrain(7,1:6) ) / 2.0
684  edstress(19,1:6) = ( edstress(3,1:6) + edstress(7,1:6) ) / 2.0
685  tdstrain(19,1:6) = ( tdstrain(3,1:6) + tdstrain(7,1:6) ) / 2.0
686  edstrain(20,1:6) = ( edstrain(4,1:6) + edstrain(8,1:6) ) / 2.0
687  edstress(20,1:6) = ( edstress(4,1:6) + edstress(8,1:6) ) / 2.0
688  tdstrain(20,1:6) = ( tdstrain(4,1:6) + tdstrain(8,1:6) ) / 2.0
689  endif
690  end subroutine nodalstress_inv3
691 
692  function get_mises(s)
693  implicit none
694  real(kind=kreal) :: get_mises, s(1:6)
695  real(kind=kreal) :: s11, s22, s33, s12, s23, s13, ps, smises
696 
697  s11 = s(1)
698  s22 = s(2)
699  s33 = s(3)
700  s12 = s(4)
701  s23 = s(5)
702  s13 = s(6)
703  ps = ( s11 + s22 + s33 ) / 3.0d0
704  smises = 0.5d0 * ( (s11-ps)**2 + (s22-ps)**2 + (s33-ps)**2 ) + s12**2 + s23**2 + s13**2
705  get_mises = dsqrt( 3.0d0 * smises )
706 
707  end function get_mises
708 
709  function get_pl_estrain(gausses)
710  implicit none
711  real(kind=kreal) :: get_pl_estrain
712  type(tgaussstatus) :: gausses(:)
713  integer(kind=kint) :: i
714 
715  get_pl_estrain = 0.d0
716  do i = 1, size(gausses)
717  get_pl_estrain = get_pl_estrain + gausses(i)%plstrain
718  enddo
719  get_pl_estrain = get_pl_estrain / size(gausses)
720 
721  end function get_pl_estrain
722 
724  !----------------------------------------------------------------------*
725  subroutine fstr_nodalstress2d( hecMESH, fstrSOLID )
726  !----------------------------------------------------------------------*
727  use m_static_lib
728  type (hecmwst_local_mesh) :: hecMESH
729  type (fstr_solid) :: fstrSOLID
730  real(kind=kreal), pointer :: tnstrain(:), testrain(:)
731  !C** local variables
732  integer(kind=kint) :: itype, icel, ic, is, iE, jS, i, j, ic_type, nn, ni, ID_area
733  real(kind=kreal) :: estrain(4), estress(4), tstrain(4), naturalcoord(4)
734  real(kind=kreal) :: edstrain(8,4), edstress(8,4), tdstrain(8,4)
735  real(kind=kreal) :: s11, s22, s33, s12, s23, s13, ps, smises
736  real(kind=kreal), allocatable :: func(:,:), inv_func(:,:)
737  integer(kind=kint), allocatable :: nnumber(:)
738 
739  tnstrain => fstrsolid%tnstrain
740  testrain => fstrsolid%testrain
741  call fstr_solid_phys_clear(fstrsolid)
742 
743  allocate( nnumber(hecmesh%n_node) )
744  if( .not. associated(fstrsolid%is_rot) ) allocate( fstrsolid%is_rot(hecmesh%n_node) )
745  nnumber = 0
746  fstrsolid%is_rot = 0
747 
748  !C +-------------------------------+
749  !C | according to ELEMENT TYPE |
750  !C +-------------------------------+
751  do itype = 1, hecmesh%n_elem_type
752  is = hecmesh%elem_type_index(itype-1) + 1
753  ie = hecmesh%elem_type_index(itype )
754  ic_type = hecmesh%elem_type_item(itype)
755  if( .not. hecmw_is_etype_surface(ic_type) ) cycle
756  !C** set number of nodes and shape function
757  nn = hecmw_get_max_node( ic_type )
758  ni = numofquadpoints( ic_type )
759  allocate( func(ni,nn), inv_func(nn,ni) )
760  if( ic_type == fe_tri6n ) then
761  ic = hecmw_get_max_node( fe_tri3n )
762  do i = 1, ni
763  call getquadpoint( ic_type, i, naturalcoord )
764  call getshapefunc( fe_tri3n, naturalcoord, func(i,1:ic) )
765  enddo
766  call inverse_func( ic, func, inv_func )
767  else if( ic_type == fe_quad4n ) then
768  do i = 1, ni
769  call getquadpoint( ic_type, i, naturalcoord )
770  call getshapefunc( ic_type, naturalcoord, func(i,1:nn) )
771  enddo
772  call inverse_func( ni, func, inv_func )
773  else if( ic_type == fe_quad8n ) then
774  ic = 0
775  do i = 1, ni
776  if( i==1 .or. i==3 .or. i==7 .or. i==9 ) then
777  ic = ic + 1
778  call getquadpoint( ic_type, i, naturalcoord )
779  call getshapefunc( fe_quad4n, naturalcoord, func(ic,1:4) )
780  endif
781  enddo
782  call inverse_func( ic, func, inv_func )
783  ni = ic
784  endif
785  !C** element loop
786  do icel = is, ie
787  js = hecmesh%elem_node_index(icel-1)
788  id_area = hecmesh%elem_ID(icel*2)
789  !--- calculate nodal stress and strain
790  if( ic_type == fe_tri6n .or. ic_type == fe_quad4n .or. ic_type == fe_quad8n ) then
791  call nodalstress_inv2( ic_type, ni, fstrsolid%elements(icel)%gausses, &
792  inv_func, edstrain(1:nn,1:4), edstress(1:nn,1:4), &
793  tdstrain(1:nn,1:4) )
794  else
795  call nodalstress_c2( ic_type, nn, fstrsolid%elements(icel)%gausses, &
796  edstrain(1:nn,1:4), edstress(1:nn,1:4) )
797  ! call NodalStress_C2( ic_type, nn, fstrSOLID%elements(icel)%gausses, &
798  ! edstrain(1:nn,1:4), edstress(1:nn,1:4), tdstrain(1:nn,1:4) )
799  endif
800  do j = 1, nn
801  ic = hecmesh%elem_node_item(js+j)
802  fstrsolid%STRAIN(3*ic-2) = fstrsolid%STRAIN(3*ic-2) + edstrain(j,1)
803  fstrsolid%STRAIN(3*ic-1) = fstrsolid%STRAIN(3*ic-1) + edstrain(j,2)
804  fstrsolid%STRAIN(3*ic-0) = fstrsolid%STRAIN(3*ic-0) + edstrain(j,3)
805  fstrsolid%STRESS(3*ic-2) = fstrsolid%STRESS(3*ic-2) + edstress(j,1)
806  fstrsolid%STRESS(3*ic-1) = fstrsolid%STRESS(3*ic-1) + edstress(j,2)
807  fstrsolid%STRESS(3*ic-0) = fstrsolid%STRESS(3*ic-0) + edstress(j,3)
808 
809  if( associated(tnstrain) ) then
810  tnstrain(3*ic-2) = tnstrain(3*ic-2) + tdstrain(j,1)
811  tnstrain(3*ic-1) = tnstrain(3*ic-1) + tdstrain(j,2)
812  tnstrain(3*ic ) = tnstrain(3*ic ) + tdstrain(j,3)
813  endif
814  nnumber(ic) = nnumber(ic) + 1
815  enddo
816  !--- calculate elemental stress and strain
817  ! if( ID_area == hecMESH%my_rank ) then
818  call elementstress_c2( ic_type, fstrsolid%elements(icel)%gausses, estrain, estress )
819  ! call ElementStress_C2( ic_type, fstrSOLID%elements(icel)%gausses, estrain, estress, tstrain )
820 
821  fstrsolid%ESTRAIN(3*icel-2) = estrain(1)
822  fstrsolid%ESTRAIN(3*icel-1) = estrain(2)
823  fstrsolid%ESTRAIN(3*icel-0) = estrain(3)
824  fstrsolid%ESTRESS(3*icel-2) = estress(1)
825  fstrsolid%ESTRESS(3*icel-1) = estress(2)
826  fstrsolid%ESTRESS(3*icel-0) = estress(3)
827 
828  !if( associated(testrain) ) then
829  ! testrain(3*icel-2) = tstrain(1)
830  ! testrain(3*icel-1) = tstrain(2)
831  ! testrain(3*icel ) = tstrain(3)
832  !endif
833  s11 = estress(1)
834  s22 = estress(2)
835  s12 = estress(3)
836  smises = 0.5d0 * ((s11-s22)**2+(s11)**2+(s22)**2) + 3*s12**2
837  fstrsolid%EMISES(icel) = sqrt( smises )
838  ! endif
839  enddo
840  deallocate( func, inv_func )
841  enddo
842 
843  !C** average over nodes
844  do i = 1, hecmesh%n_node
845  if( nnumber(i) == 0 ) cycle
846  fstrsolid%STRAIN(3*i-2:3*i-0) = fstrsolid%STRAIN(3*i-2:3*i-0) / nnumber(i)
847  fstrsolid%STRESS(3*i-2:3*i-0) = fstrsolid%STRESS(3*i-2:3*i-0) / nnumber(i)
848  if( associated(tnstrain) ) tnstrain(3*i-2:3*i) = tnstrain(3*i-2:3*i) / nnumber(i)
849  enddo
850  !C** calculate von MISES stress
851  do i = 1, hecmesh%n_node
852  s11 = fstrsolid%STRESS(3*i-2)
853  s22 = fstrsolid%STRESS(3*i-1)
854  s12 = fstrsolid%STRESS(3*i-0)
855  smises = 0.5d0 * ((s11-s22)**2+(s11)**2+(s22)**2) + 3*s12**2
856  fstrsolid%MISES(i) = sqrt( smises )
857  enddo
858 
859  deallocate( nnumber )
860  end subroutine fstr_nodalstress2d
861 
862  !----------------------------------------------------------------------*
863  subroutine nodalstress_inv2( etype, ni, gausses, func, edstrain, edstress, tdstrain )
864  !----------------------------------------------------------------------*
865  use mmechgauss
866  integer(kind=kint) :: etype, ni
867  type(tgaussstatus) :: gausses(:)
868  real(kind=kreal) :: func(:,:), edstrain(:,:), edstress(:,:), tdstrain(:,:)
869  integer :: i, j, k, ic
870 
871  edstrain = 0.0d0
872  edstress = 0.0d0
873  tdstrain = 0.0d0
874 
875  if( etype == fe_quad4n ) then
876  do i = 1, ni
877  do j = 1, ni
878  do k = 1, 4
879  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
880  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
881  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
882  enddo
883  enddo
884  enddo
885  else if( etype == fe_tri6n ) then
886  do i = 1, ni
887  do j = 1, ni
888  do k = 1, 4
889  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
890  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
891  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
892  enddo
893  enddo
894  enddo
895  edstrain(4,1:4) = ( edstrain(1,1:4) + edstrain(2,1:4) ) / 2.0
896  edstress(4,1:4) = ( edstress(1,1:4) + edstress(2,1:4) ) / 2.0
897  tdstrain(4,1:4) = ( tdstrain(1,1:4) + tdstrain(2,1:4) ) / 2.0
898  edstrain(5,1:4) = ( edstrain(2,1:4) + edstrain(3,1:4) ) / 2.0
899  edstress(5,1:4) = ( edstress(2,1:4) + edstress(3,1:4) ) / 2.0
900  tdstrain(5,1:4) = ( tdstrain(2,1:4) + tdstrain(3,1:4) ) / 2.0
901  edstrain(6,1:4) = ( edstrain(3,1:4) + edstrain(1,1:4) ) / 2.0
902  edstress(6,1:4) = ( edstress(3,1:4) + edstress(1,1:4) ) / 2.0
903  tdstrain(6,1:4) = ( tdstrain(3,1:4) + tdstrain(1,1:4) ) / 2.0
904  else if( etype == fe_quad8n ) then
905  do i = 1, ni
906  ic = 0
907  do j = 1, numofquadpoints(etype)
908  if( j==1 .or. j==3 .or. j==7 .or. j==9 ) then
909  ic = ic + 1
910  do k = 1, 4
911  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
912  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
913  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
914  enddo
915  endif
916  enddo
917  enddo
918  edstrain(5,1:4) = ( edstrain(1,1:4) + edstrain(2,1:4) ) / 2.0
919  edstress(5,1:4) = ( edstress(1,1:4) + edstress(2,1:4) ) / 2.0
920  tdstrain(5,1:4) = ( tdstrain(1,1:4) + tdstrain(2,1:4) ) / 2.0
921  edstrain(6,1:4) = ( edstrain(2,1:4) + edstrain(3,1:4) ) / 2.0
922  edstress(6,1:4) = ( edstress(2,1:4) + edstress(3,1:4) ) / 2.0
923  tdstrain(6,1:4) = ( tdstrain(2,1:4) + tdstrain(3,1:4) ) / 2.0
924  edstrain(7,1:4) = ( edstrain(3,1:4) + edstrain(4,1:4) ) / 2.0
925  edstress(7,1:4) = ( edstress(3,1:4) + edstress(4,1:4) ) / 2.0
926  tdstrain(7,1:4) = ( tdstrain(3,1:4) + tdstrain(4,1:4) ) / 2.0
927  edstrain(8,1:4) = ( edstrain(4,1:4) + edstrain(1,1:4) ) / 2.0
928  edstress(8,1:4) = ( edstress(4,1:4) + edstress(1,1:4) ) / 2.0
929  tdstrain(8,1:4) = ( tdstrain(4,1:4) + tdstrain(1,1:4) ) / 2.0
930  endif
931  end subroutine nodalstress_inv2
932 
933  !----------------------------------------------------------------------*
934  subroutine inverse_func( n, a, inv_a )
935  !----------------------------------------------------------------------*
936  integer(kind=kint) :: n
937  real(kind=kreal) :: a(:,:), inv_a(:,:)
938  integer(kind=kint) :: i, j, k
939  real(kind=kreal) :: buf
940 
941  do i = 1, n
942  do j = 1, n
943  if( i == j ) then
944  inv_a(i,j) = 1.0
945  else
946  inv_a(i,j) = 0.0
947  endif
948  enddo
949  enddo
950 
951  do i = 1, n
952  buf = 1.0 / a(i,i)
953  do j = 1, n
954  a(i,j) = a(i,j) * buf
955  inv_a(i,j) = inv_a(i,j) *buf
956  enddo
957  do j = 1, n
958  if( i /= j ) then
959  buf = a(j,i)
960  do k = 1, n
961  a(j,k) = a(j,k) - a(i,k) * buf
962  inv_a(j,k) = inv_a(j,k) - inv_a(i,k) * buf
963  enddo
964  endif
965  enddo
966  enddo
967  end subroutine inverse_func
968 
970  !----------------------------------------------------------------------*
971  subroutine fstr_nodalstress6d( hecMESH, fstrSOLID )
972  !----------------------------------------------------------------------*
973  use m_static_lib
974  type (hecmwST_local_mesh) :: hecMESH
975  type (fstr_solid) :: fstrSOLID
976  !C** local variables
977  integer(kind=kint) :: itype, icel, is, iE, jS, i, j, k, it, ic, ic_type, nn, isect, ihead, ID_area
978  integer(kind=kint) :: nodLOCAL(20), n_layer, ntot_lyr, nlyr, n_totlyr, com_total_layer, shellmatl
979  real(kind=kreal) :: ecoord(3,9), edisp(6,9), estrain(6), estress(6), ndstrain(9,6), ndstress(9,6)
980  real(kind=kreal) :: thick, thick_layer
981  real(kind=kreal) :: s11, s22, s33, s12, s23, s13, t11, t22, t33, t12, t23, t13, ps, smises, tmises
982  integer(kind=kint), allocatable :: nnumber(:)
983  type(fstr_solid_physic_val), pointer :: layer => null()
984 
985  call fstr_solid_phys_clear(fstrsolid)
986 
987  n_totlyr = fstrsolid%max_lyr
988 
989  allocate( nnumber(hecmesh%n_node) )
990  if( .not. associated(fstrsolid%is_rot) ) allocate( fstrsolid%is_rot(hecmesh%n_node) )
991  nnumber = 0
992  fstrsolid%is_rot = 0
993 
994  !C +-------------------------------+
995  !C | according to ELEMENT TYPE |
996  !C +-------------------------------+
997  do itype = 1, hecmesh%n_elem_type
998  is = hecmesh%elem_type_index(itype-1) + 1
999  ie = hecmesh%elem_type_index(itype )
1000  ic_type = hecmesh%elem_type_item(itype)
1001  if( .not. hecmw_is_etype_shell(ic_type) ) then
1002  ntot_lyr = 0
1003  cycle
1004  end if
1005  nn = hecmw_get_max_node( ic_type )
1006  !C** element loop
1007  do icel = is, ie
1008  js = hecmesh%elem_node_index(icel-1)
1009  id_area = hecmesh%elem_ID(icel*2)
1010  do j = 1, nn
1011  nodlocal(j) = hecmesh%elem_node_item(js+j)
1012  ecoord(1,j) = hecmesh%node(3*nodlocal(j)-2)
1013  ecoord(2,j) = hecmesh%node(3*nodlocal(j)-1)
1014  ecoord(3,j) = hecmesh%node(3*nodlocal(j) )
1015  edisp(1,j) = fstrsolid%unode(6*nodlocal(j)-5)
1016  edisp(2,j) = fstrsolid%unode(6*nodlocal(j)-4)
1017  edisp(3,j) = fstrsolid%unode(6*nodlocal(j)-3)
1018  edisp(4,j) = fstrsolid%unode(6*nodlocal(j)-2)
1019  edisp(5,j) = fstrsolid%unode(6*nodlocal(j)-1)
1020  edisp(6,j) = fstrsolid%unode(6*nodlocal(j) )
1021  enddo
1022  isect = hecmesh%section_ID(icel)
1023  ihead = hecmesh%section%sect_R_index(isect-1)
1024  thick = hecmesh%section%sect_R_item(ihead+1)
1025  !--- calculate elemental stress and strain
1026  if( ic_type == 731 .or. ic_type == 741 .or. ic_type == 743 ) then
1027  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
1028  do nlyr=1,ntot_lyr
1029  call elementstress_shell_mitc( ic_type, nn, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
1030  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), thick, 1.0d0, nlyr, ntot_lyr)
1031  do j = 1, nn
1032  i = nodlocal(j)
1033  layer => fstrsolid%SHELL%LAYER(nlyr)%PLUS
1034  do k = 1, 6
1035  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + ndstrain(j,k)
1036  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + ndstress(j,k)
1037  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + ndstrain(j,k)/nn
1038  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + ndstress(j,k)/nn
1039  enddo
1040  enddo
1041  !minus section
1042  call elementstress_shell_mitc( ic_type, nn, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
1043  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), thick,-1.0d0, nlyr, ntot_lyr)
1044  do j = 1, nn
1045  i = nodlocal(j)
1046  layer => fstrsolid%SHELL%LAYER(nlyr)%MINUS
1047  do k = 1, 6
1048  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + ndstrain(j,k)
1049  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + ndstress(j,k)
1050  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + ndstrain(j,k)/nn
1051  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + ndstress(j,k)/nn
1052  enddo
1053  enddo
1054  enddo
1055  call fstr_getavg_shell(nn,fstrsolid,icel,nodlocal,ndstrain(1:nn,1:6),ndstress(1:nn,1:6),estrain,estress)
1056  endif
1057 
1058  !if( ID_area == hecMESH%my_rank ) then
1059  !ADD VALUE and Count node
1060  do j = 1, nn
1061  ic = hecmesh%elem_node_item(js+j)
1062  fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) + ndstrain(j,1:6)
1063  fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) + ndstress(j,1:6)
1064  !if( associated(tnstrain) )then
1065  ! tnstrain(6*(ic-1)+1:6*(ic-1)+6) = tnstrain(6*(ic-1)+1:6*(ic-1)+6) + tdstrain(j,1:6)
1066  !endif
1067  nnumber(ic) = nnumber(ic) + 1
1068  enddo
1069 
1070  fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) + estrain(1:6)
1071  fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) + estress(1:6)
1072  !endif
1073  enddo
1074  enddo
1075 
1076  !C** calculate nodal stress and strain
1077  do i = 1, hecmesh%n_node
1078  if( nnumber(i) == 0 ) cycle
1079  fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1080  fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1081  !if( associated(tnstrain) )then
1082  ! tnstrain(6*(i-1)+1:6*(i-1)+6) = tnstrain(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1083  !endif
1084  enddo
1085 
1086  do nlyr = 1, ntot_lyr
1087  do i = 1, hecmesh%n_node
1088  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
1089  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1090  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
1091  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1092  fstrsolid%SHELL%LAYER(nlyr)%PLUS%MISES(i) = &
1093  & get_mises(fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6))
1094 
1095  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
1096  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1097  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
1098  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
1099  fstrsolid%SHELL%LAYER(nlyr)%MINUS%MISES(i) = &
1100  & get_mises(fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6))
1101  enddo
1102  enddo
1103 
1104  !C** calculate von MISES stress
1105  do i = 1, hecmesh%n_node
1106  fstrsolid%MISES(i) = get_mises(fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6))
1107  enddo
1108  do i = 1, hecmesh%n_elem
1109  fstrsolid%EMISES(i) = get_mises(fstrsolid%ESTRESS(6*(i-1)+1:6*(i-1)+6))
1110  enddo
1111 
1112  !C** calculate Elemental Plastic Strain
1113  do i = 1, hecmesh%n_elem
1114  if (.not. associated(fstrsolid%elements(i)%gausses)) cycle
1115  fstrsolid%EPLSTRAIN(i) = get_pl_estrain(fstrsolid%elements(i)%gausses)
1116  enddo
1117 
1118  deallocate( nnumber )
1119 
1120  end subroutine fstr_nodalstress6d
1121 
1122  subroutine make_principal(fstrSOLID, hecMESH, RES)
1123  use hecmw_util
1124  use m_out
1125  use m_static_lib
1126 
1127  type(fstr_solid) :: fstrSOLID
1128  type(hecmwst_local_mesh) :: hecMESH
1129  type(fstr_solid_physic_val) :: RES
1130  integer(kind=kint) :: i, flag
1131  real(kind=kreal) :: tmat(3, 3), tvec(3), strain(6)
1132 
1133  flag=ieor(flag,flag)
1134  if( fstrsolid%output_ctrl(3)%outinfo%on(19) .or. fstrsolid%output_ctrl(4)%outinfo%on(19) ) then
1135  if ( .not. associated(res%PSTRESS) ) then
1136  allocate(res%PSTRESS( 3*hecmesh%n_node ))
1137  endif
1138  flag=ior(flag,b'00000001')
1139  end if
1140  if( fstrsolid%output_ctrl(3)%outinfo%on(23) .or. fstrsolid%output_ctrl(4)%outinfo%on(23) ) then
1141  if ( .not. associated(res%PSTRESS_VECT) ) then
1142  allocate(res%PSTRESS_VECT( 3*hecmesh%n_node ,3))
1143  endif
1144  flag=ior(flag,b'00000010')
1145  end if
1146  if( fstrsolid%output_ctrl(3)%outinfo%on(21) .or. fstrsolid%output_ctrl(4)%outinfo%on(21) ) then
1147  if ( .not. associated(res%PSTRAIN) ) then
1148  allocate(res%PSTRAIN( 3*hecmesh%n_node ))
1149  endif
1150  flag=ior(flag,b'00000100')
1151  end if
1152  if( fstrsolid%output_ctrl(3)%outinfo%on(25) .or. fstrsolid%output_ctrl(4)%outinfo%on(25) ) then
1153  if ( .not. associated(res%PSTRAIN_VECT) ) then
1154  allocate(res%PSTRAIN_VECT( 3*hecmesh%n_node ,3))
1155  endif
1156  flag=ior(flag,b'00001000')
1157  end if
1158  if( fstrsolid%output_ctrl(3)%outinfo%on(20) .or. fstrsolid%output_ctrl(4)%outinfo%on(20) ) then
1159  if ( .not. associated(res%EPSTRESS) ) then
1160  allocate(res%EPSTRESS( 3*hecmesh%n_elem ))
1161  endif
1162  flag=ior(flag,b'00010000')
1163  end if
1164  if( fstrsolid%output_ctrl(3)%outinfo%on(24) .or. fstrsolid%output_ctrl(4)%outinfo%on(24) ) then
1165  if ( .not. associated(res%EPSTRESS_VECT) ) then
1166  allocate(res%EPSTRESS_VECT( 3*hecmesh%n_elem ,3))
1167  endif
1168  flag=ior(flag,b'00100000')
1169  end if
1170  if( fstrsolid%output_ctrl(3)%outinfo%on(22) .or. fstrsolid%output_ctrl(4)%outinfo%on(22) ) then
1171  if ( .not. associated(res%EPSTRAIN) ) then
1172  allocate(res%EPSTRAIN( 3*hecmesh%n_elem ))
1173  endif
1174  flag=ior(flag,b'01000000')
1175  end if
1176  if( fstrsolid%output_ctrl(3)%outinfo%on(26) .or. fstrsolid%output_ctrl(4)%outinfo%on(26) ) then
1177  if ( .not. associated(res%EPSTRAIN_VECT) ) then
1178  allocate(res%EPSTRAIN_VECT( 3*hecmesh%n_elem ,3))
1179  endif
1180  flag=ior(flag,b'10000000')
1181  end if
1182 
1183  if (iand(flag,b'00000011') /= 0) then
1184  do i = 1, hecmesh%n_node
1185  call get_principal(res%STRESS(6*i-5:6*i), tvec, tmat)
1186  if (iand(flag,b'00000001') /= 0) res%PSTRESS(3*(i-1)+1:3*(i-1)+3)=tvec
1187  if (iand(flag,b'00000010') /= 0) res%PSTRESS_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
1188  end do
1189  end if
1190  if (iand(flag,b'00001100') /= 0) then
1191  do i = 1, hecmesh%n_node
1192  strain(1:6) = res%STRAIN(6*i-5:6*i)
1193  strain(4:6) = 0.5d0*strain(4:6)
1194  call get_principal(strain, tvec, tmat)
1195  if (iand(flag,b'00000100') /= 0) res%PSTRAIN(3*(i-1)+1:3*(i-1)+3)=tvec
1196  if (iand(flag,b'00001000') /= 0) res%PSTRAIN_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
1197  end do
1198  end if
1199 
1200  if (iand(flag,b'00110000') /= 0) then
1201  do i = 1, hecmesh%n_elem
1202  call get_principal( res%ESTRESS(6*i-5:6*i), tvec, tmat)
1203  if (iand(flag,b'00010000') /= 0) res%EPSTRESS(3*(i-1)+1:3*(i-1)+3)=tvec
1204  if (iand(flag,b'00100000') /= 0) res%EPSTRESS_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
1205  end do
1206  end if
1207  if (iand(flag,b'11000000') /= 0) then
1208  do i = 1, hecmesh%n_elem
1209  strain(1:6) = res%ESTRAIN(6*i-5:6*i)
1210  strain(4:6) = 0.5d0*strain(4:6)
1211  call get_principal(strain, tvec, tmat)
1212  if (iand(flag,b'01000000') /= 0) res%EPSTRAIN(3*(i-1)+1:3*(i-1)+3)=tvec
1213  if (iand(flag,b'10000000') /= 0) res%EPSTRAIN_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
1214  end do
1215  end if
1216  end subroutine make_principal
1217 
1218 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_nodalstress
This module provides functions to calculation nodal stress.
Definition: fstr_NodalStress.f90:6
m_fstr_nodalstress::fstr_stress_add_shelllyr
subroutine fstr_stress_add_shelllyr(nn, fstrSOLID, icel, nodLOCAL, nlyr, strain, stress, flag)
Definition: fstr_NodalStress.f90:487
m_fstr_nodalstress::get_pl_estrain
real(kind=kreal) function get_pl_estrain(gausses)
Definition: fstr_NodalStress.f90:710
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:972
m_fstr_nodalstress::get_mises
real(kind=kreal) function get_mises(s)
Definition: fstr_NodalStress.f90:693
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:312
m_fstr_nodalstress::fstr_nodalstress2d
subroutine fstr_nodalstress2d(hecMESH, fstrSOLID)
Calculate NODAL STRESS of plane elements.
Definition: fstr_NodalStress.f90:726
m_fstr_nodalstress::make_principal
subroutine make_principal(fstrSOLID, hecMESH, RES)
Definition: fstr_NodalStress.f90:1123
m_fstr_nodalstress::search_idx_senes
integer(kind=kint) function search_idx_senes(irow, asect, nid, sid)
Definition: fstr_NodalStress.f90:293
m_fstr_nodalstress::fstr_getavg_shell
subroutine fstr_getavg_shell(nn, fstrSOLID, icel, nodLOCAL, strain, stress, estrain, estress)
Definition: fstr_NodalStress.f90:514
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:76
m_fstr::fstr_solid_phys_clear
subroutine fstr_solid_phys_clear(fstrSOLID)
Definition: m_fstr.f90:1156
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:234
m_fstr::fstr_solid_physic_val
Data for STATIC ANSLYSIS (fstrSOLID)
Definition: m_fstr.f90:213
m_fstr::fstr_solid
Definition: m_fstr.f90:240
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13