FrontISTR  5.2.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 !-------------------------------------------------------------------------------
7  implicit none
8  private :: nodalstress_inv3, nodalstress_inv2, inverse_func
9 contains
10 
12  !----------------------------------------------------------------------*
13  subroutine fstr_nodalstress3d( hecMESH, fstrSOLID )
14  !----------------------------------------------------------------------*
15  use m_fstr
16  use m_static_lib
17  type(hecmwst_local_mesh) :: hecMESH
18  type(fstr_solid) :: fstrSOLID
19  real(kind=kreal), pointer :: tnstrain(:), testrain(:), yield_ratio(:)
20  integer(kind=kint), pointer :: is_rot(:)
21  !C** local variables
22  integer(kind=kint) :: itype, icel, ic, is, iE, jS, i, j, k, m, ic_type, nn, ni, ID_area
23  integer(kind=kint) :: nodlocal(20), ntemp
24  integer(kind=kint), allocatable :: nnumber(:)
25  real(kind=kreal) :: estrain(6), estress(6), naturalcoord(3)
26  real(kind=kreal) :: enqm(12)
27  real(kind=kreal) :: ndstrain(20,6), ndstress(20,6), tdstrain(20,6)
28  real(kind=kreal) :: ecoord(3, 20), edisp(60), tt(20), t0(20)
29  real(kind=kreal), allocatable :: func(:,:), inv_func(:,:)
30 
31  !C** Shell33 variables
32  integer(kind=kint) :: isect, ihead, ntot_lyr, nlyr, flag33, cid, truss
33  real(kind=kreal) :: thick, thick_lyr, dtot_lyr
34  call fstr_solid_phys_clear(fstrsolid)
35 
36  allocate( nnumber(hecmesh%n_node) )
37  if( .not. associated(fstrsolid%is_rot) ) allocate( fstrsolid%is_rot(hecmesh%n_node) )
38  !allocate( fstrSOLID%yield_ratio(hecMESH%n_elem) )
39  nnumber = 0
40  fstrsolid%is_rot = 0
41  !fstrSOLID%yield_ratio = 0.0d0
42 
43  tnstrain => fstrsolid%tnstrain
44  testrain => fstrsolid%testrain
45  is_rot => fstrsolid%is_rot
46  yield_ratio => fstrsolid%yield_ratio
47 
48  if( associated(tnstrain) ) tnstrain = 0.0d0
49 
50  !C** setting
51  ntot_lyr = fstrsolid%max_lyr
52  flag33 = fstrsolid%is_33shell
53  truss = fstrsolid%is_33beam
54 
55  !C +-------------------------------+
56  !C | according to ELEMENT TYPE |
57  !C +-------------------------------+
58  do itype = 1, hecmesh%n_elem_type
59  is = hecmesh%elem_type_index(itype-1) + 1
60  ie = hecmesh%elem_type_index(itype )
61  ic_type = hecmesh%elem_type_item(itype)
62  if( ic_type == fe_tet10nc ) ic_type = fe_tet10n
63  if( .not. (hecmw_is_etype_solid(ic_type) .or. ic_type == 781 &
64  & .or. ic_type == 761 .or. ic_type == fe_beam341 ) ) cycle
65  !C** set number of nodes and shape function
66  nn = hecmw_get_max_node( ic_type )
67  ni = numofquadpoints( ic_type )
68  allocate( func(ni,nn), inv_func(nn,ni) )
69  if( ic_type == fe_tet10n ) then
70  ic = hecmw_get_max_node( fe_tet4n )
71  do i = 1, ni
72  call getquadpoint( ic_type, i, naturalcoord )
73  call getshapefunc( fe_tet4n, naturalcoord, func(i,1:ic) )
74  enddo
75  call inverse_func( ic, func, inv_func )
76  else if( ic_type == fe_hex8n ) then
77  do i = 1, ni
78  call getquadpoint( ic_type, i, naturalcoord )
79  call getshapefunc( ic_type, naturalcoord, func(i,1:nn) )
80  enddo
81  call inverse_func( ni, func, inv_func )
82  else if( ic_type == fe_prism15n ) then
83  ic = 0
84  do i = 1, ni
85  if( i==1 .or. i==2 .or. i==3 .or. i==7 .or. i==8 .or. i==9 ) then
86  ic = ic + 1
87  call getquadpoint( ic_type, i, naturalcoord )
88  call getshapefunc( fe_prism6n, naturalcoord, func(ic,1:6) )
89  endif
90  enddo
91  call inverse_func( ic, func, inv_func )
92  ni = ic
93  else if( ic_type == fe_hex20n ) then
94  ic = 0
95  do i = 1, ni
96  if( i==1 .or. i==3 .or. i==7 .or. i==9 .or. &
97  i==19 .or. i==21 .or. i==25 .or. i==27 ) then
98  ic = ic + 1
99  call getquadpoint( ic_type, i, naturalcoord )
100  call getshapefunc( fe_hex8n, naturalcoord, func(ic,1:8) )
101  endif
102  enddo
103  call inverse_func( ic, func, inv_func )
104  ni = ic
105  endif
106  !C** element loop
107  do icel = is, ie
108  js = hecmesh%elem_node_index(icel-1)
109  id_area = hecmesh%elem_ID(icel*2)
110  isect= hecmesh%section_ID(icel)
111  ihead = hecmesh%section%sect_R_index(isect-1)
112  thick = hecmesh%section%sect_R_item(ihead+1)
113  !initialize
114  enqm = 0.0d0
115  estrain = 0.0d0
116  estress = 0.0d0
117  ndstrain = 0.0d0
118  ndstress = 0.0d0
119  !if( ID_area == hecMESH%my_rank ) then
120 
121  !--- calculate nodal and elemental value
122  if( ic_type == 641 ) then
123  do j = 1, 4
124  nodlocal(j) = hecmesh%elem_node_item(js+j)
125  ecoord(1:3,j) = hecmesh%node(3*nodlocal(j)-2:3*nodlocal(j))
126  edisp(3*j-2:3*j) = fstrsolid%unode(3*nodlocal(j)-2:3*nodlocal(j))
127  end do
128  ntemp = 0
129  if( associated( fstrsolid%temperature ) ) then
130  ntemp = 1
131  do j = 1, 4
132  nodlocal(j) = hecmesh%elem_node_item(js+j)
133  t0(j) = fstrsolid%last_temp( nodlocal(j) )
134  tt(j) = fstrsolid%temperature( nodlocal(j) )
135  end do
136  end if
137  call nodalstress_beam_641( ic_type, nn, ecoord, fstrsolid%elements(icel)%gausses, &
138  & hecmesh%section%sect_R_item(ihead+1:), edisp, &
139  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), tt(1:nn), t0(1:nn), ntemp )
140  call elementalstress_beam_641( fstrsolid%elements(icel)%gausses, estrain, estress, enqm )
141  fstrsolid%ENQM(icel*12-11:icel*12) = enqm(1:12)
142 
143 
144  elseif( ic_type == 781) then
145  do j = 1, 4
146  nodlocal(j ) = hecmesh%elem_node_item(js+j )
147  nodlocal(j+4) = hecmesh%elem_node_item(js+j+4)
148  is_rot(nodlocal(j+4)) = 1
149  ecoord(1:3,j ) = hecmesh%node(3*nodlocal(j )-2:3*nodlocal(j ))
150  ecoord(1:3,j+4) = hecmesh%node(3*nodlocal(j+4)-2:3*nodlocal(j+4))
151  edisp(6*j-5:6*j-3) = fstrsolid%unode(3*nodlocal(j )-2:3*nodlocal(j ))
152  edisp(6*j-2:6*j ) = fstrsolid%unode(3*nodlocal(j+4)-2:3*nodlocal(j+4))
153  enddo
154  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
155  do nlyr=1,ntot_lyr
156  call elementstress_shell_mitc( 741, 4, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
157  & ndstrain(1:4,1:6), ndstress(1:4,1:6), thick, 1.0d0, nlyr, ntot_lyr)
158  call fstr_stress_add_shelllyr(4,fstrsolid,icel,nodlocal,nlyr,ndstrain(1:4,1:6),ndstress(1:4,1:6),1)
159  !minus section
160  call elementstress_shell_mitc( 741, 4, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
161  & ndstrain(1:4,1:6), ndstress(1:4,1:6), thick,-1.0d0, nlyr, ntot_lyr)
162  call fstr_stress_add_shelllyr(4,fstrsolid,icel,nodlocal,nlyr,ndstrain(1:4,1:6),ndstress(1:4,1:6),-1)
163  enddo
164  call fstr_getavg_shell(4,fstrsolid,icel,nodlocal,ndstrain(1:4,1:6),ndstress(1:4,1:6),estrain,estress)
165 
166  elseif( ic_type == 761) then
167  do j = 1, 3
168  nodlocal(j ) = hecmesh%elem_node_item(js+j )
169  nodlocal(j+3) = hecmesh%elem_node_item(js+j+3)
170  is_rot(nodlocal(j+3)) = 1
171  ecoord(1:3,j ) = hecmesh%node(3*nodlocal(j )-2:3*nodlocal(j ))
172  ecoord(1:3,j+3) = hecmesh%node(3*nodlocal(j+3)-2:3*nodlocal(j+3))
173  edisp(6*j-5:6*j-3) = fstrsolid%unode(3*nodlocal(j )-2:3*nodlocal(j ))
174  edisp(6*j-2:6*j ) = fstrsolid%unode(3*nodlocal(j+3)-2:3*nodlocal(j+3))
175  enddo
176  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
177  do nlyr=1,ntot_lyr
178  call elementstress_shell_mitc( 731, 3, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
179  & ndstrain(1:3,1:6), ndstress(1:3,1:6), thick, 1.0d0, nlyr, ntot_lyr)
180  call fstr_stress_add_shelllyr(3,fstrsolid,icel,nodlocal,nlyr,ndstrain(1:3,1:6),ndstress(1:3,1:6),1)
181  !minus section
182  call elementstress_shell_mitc( 731, 3, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
183  & ndstrain(1:3,1:6), ndstress(1:3,1:6), thick,-1.0d0, nlyr, ntot_lyr)
184  call fstr_stress_add_shelllyr(3,fstrsolid,icel,nodlocal,nlyr,ndstrain(1:3,1:6),ndstress(1:3,1:6),-1)
185  enddo
186  call fstr_getavg_shell(3,fstrsolid,icel,nodlocal,ndstrain(1:3,1:6),ndstress(1:3,1:6),estrain,estress)
187 
188  else if( ic_type == 301 ) then
189  call nodalstress_c1( ic_type, nn, fstrsolid%elements(icel)%gausses, &
190  ndstrain(1:nn,1:6), ndstress(1:nn,1:6) )
191  call elementstress_c1( ic_type, fstrsolid%elements(icel)%gausses, estrain, estress )
192 
193  else if( ic_type == fe_tet10n .or. ic_type == fe_hex8n .or. &
194  ic_type == fe_prism15n .or. ic_type == fe_hex20n ) then
195  call nodalstress_inv3( ic_type, ni, fstrsolid%elements(icel)%gausses, &
196  inv_func, ndstrain(1:nn,1:6), ndstress(1:nn,1:6), &
197  tdstrain(1:nn,1:6) )
198  call elementstress_c3( ic_type, fstrsolid%elements(icel)%gausses, estrain, estress )
199 
200  else
201  call nodalstress_c3( ic_type, nn, fstrsolid%elements(icel)%gausses, &
202  ndstrain(1:nn,1:6), ndstress(1:nn,1:6) )
203  !call NodalStress_C3( ic_type, nn, fstrSOLID%elements(icel)%gausses, &
204  ! ndstrain(1:nn,1:6), ndstress(1:nn,1:6), tdstrain(1:nn,1:6) )
205  call elementstress_c3( ic_type, fstrsolid%elements(icel)%gausses, estrain, estress )
206 
207  endif
208 
209  !ADD VALUE and Count node
210  do j = 1, nn
211  ic = hecmesh%elem_node_item(js+j)
212  fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) + ndstrain(j,1:6)
213  fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) + ndstress(j,1:6)
214  if( associated(tnstrain) )then
215  tnstrain(6*(ic-1)+1:6*(ic-1)+6) = tnstrain(6*(ic-1)+1:6*(ic-1)+6) + tdstrain(j,1:6)
216  endif
217  nnumber(ic) = nnumber(ic) + 1
218  enddo
219 
220  fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) + estrain(1:6)
221  fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) + estress(1:6)
222 
223  !endif
224  enddo
225  deallocate( func, inv_func )
226  enddo
227 
228  !C** calculate nodal stress and strain
229  do i = 1, hecmesh%n_node
230  if( nnumber(i) == 0 ) cycle
231  fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
232  fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
233  if( associated(tnstrain) )then
234  tnstrain(6*(i-1)+1:6*(i-1)+6) = tnstrain(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
235  endif
236  enddo
237 
238  if( flag33 == 1 )then
239  do nlyr = 1, ntot_lyr
240  do i = 1, hecmesh%n_node
241  if( nnumber(i) == 0 ) cycle
242  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
243  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
244  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
245  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
246  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
247  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
248  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
249  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
250  enddo
251  enddo
252  endif
253 
254  !C** calculate von MISES stress
255  do i = 1, hecmesh%n_node
256  fstrsolid%MISES(i) = get_mises(fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6))
257  enddo
258  do i = 1, hecmesh%n_elem
259  fstrsolid%EMISES(i) = get_mises(fstrsolid%ESTRESS(6*(i-1)+1:6*(i-1)+6))
260  enddo
261 
262  if( flag33 == 1 )then
263  if( fstrsolid%output_ctrl(3)%outinfo%on(27) .or. fstrsolid%output_ctrl(4)%outinfo%on(27) ) then
264  do nlyr = 1, ntot_lyr
265  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL%LAYER(nlyr)%PLUS)
266  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL%LAYER(nlyr)%MINUS)
267  enddo
268  endif
269  call make_principal(fstrsolid, hecmesh, fstrsolid%SHELL)
270  else
271  call make_principal(fstrsolid, hecmesh, fstrsolid%SOLID)
272  endif
273 
274  deallocate( nnumber )
275 
276  end subroutine fstr_nodalstress3d
277 
278  subroutine fstr_stress_add_shelllyr(nn,fstrSOLID,icel,nodLOCAL,nlyr,strain,stress,flag)
279  use m_fstr
280  implicit none
281  type(fstr_solid) :: fstrsolid
282  integer(kind=kint) :: nodlocal(20)
283  integer(kind=kint) :: nn, i, j, k, m, nlyr, weight, icel, flag
284  real(kind=kreal) :: strain(nn, 6), stress(nn, 6)
285  type(fstr_solid_physic_val), pointer :: layer => null()
286 
287  do j = 1, nn
288  i = nodlocal(j)
289  m = nodlocal(j+nn)
290  if(flag == 1)then
291  layer => fstrsolid%SHELL%LAYER(nlyr)%PLUS
292  elseif(flag == -1)then
293  layer => fstrsolid%SHELL%LAYER(nlyr)%MINUS
294  endif
295  do k = 1, 6
296  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + strain(j,k)
297  layer%STRAIN(6*(m-1)+k) = layer%STRAIN(6*(m-1)+k) + strain(j,k)
298  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + stress(j,k)
299  layer%STRESS(6*(m-1)+k) = layer%STRESS(6*(m-1)+k) + stress(j,k)
300  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + strain(j,k)/nn
301  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + stress(j,k)/nn
302  enddo
303  enddo
304  end subroutine fstr_stress_add_shelllyr
305 
306  subroutine fstr_getavg_shell(nn,fstrSOLID,icel,nodLOCAL,strain,stress,estrain,estress)
307  use m_fstr
308  implicit none
309  type (fstr_solid) :: fstrsolid
310  integer(kind=kint) :: nodlocal(20)
311  integer(kind=kint) :: nn, i, j, k, m, nlyr, icel, flag, ntot_lyr
312  real(kind=kreal) :: strain(nn,6), stress(nn,6), estrain(6), estress(6), weight
313  type(fstr_solid_physic_val), pointer :: layer => null()
314 
315  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
316  strain = 0.0d0
317  stress = 0.0d0
318  estrain = 0.0d0
319  estress = 0.0d0
320 
321  do nlyr = 1, ntot_lyr
322  layer => fstrsolid%SHELL%LAYER(nlyr)
323  weight = fstrsolid%elements(icel)%gausses(1)%pMaterial%shell_var(nlyr)%weight
324  do j = 1, nn
325  i = nodlocal(j)
326  do k = 1, 6
327  strain(j,k) = strain(j,k) &
328  & + weight*(0.5d0*layer%PLUS%STRAIN(6*(i-1)+k) + 0.5d0*layer%MINUS%STRAIN(6*(i-1)+k))
329  stress(j,k) = stress(j,k) &
330  & + weight*(0.5d0*layer%PLUS%STRESS(6*(i-1)+k) + 0.5d0*layer%MINUS%STRESS(6*(i-1)+k))
331  enddo
332  estrain(j) = estrain(j) &
333  & + weight*(0.5d0*layer%PLUS%ESTRAIN(6*(icel-1)+j) + 0.5d0*layer%MINUS%ESTRAIN(6*(icel-1)+j))
334  estress(j) = estress(j) &
335  & + weight*(0.5d0*layer%PLUS%ESTRESS(6*(icel-1)+j) + 0.5d0*layer%MINUS%ESTRESS(6*(icel-1)+j))
336  enddo
337  enddo
338  end subroutine fstr_getavg_shell
339 
340  !----------------------------------------------------------------------*
341  subroutine nodalstress_inv3( etype, ni, gausses, func, edstrain, edstress, tdstrain )
342  !----------------------------------------------------------------------*
343  use m_fstr
344  use mmechgauss
345  integer(kind=kint) :: etype, ni
346  type(tgaussstatus) :: gausses(:)
347  real(kind=kreal) :: func(:, :), edstrain(:, :), edstress(:, :), tdstrain(:, :)
348  integer :: i, j, k, ic
349 
350  edstrain = 0.0d0
351  edstress = 0.0d0
352  tdstrain = 0.0d0
353 
354  if( etype == fe_hex8n ) then
355  do i = 1, ni
356  do j = 1, ni
357  do k = 1, 6
358  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
359  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
360  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
361  enddo
362  enddo
363  enddo
364  else if( etype == fe_tet10n ) then
365  do i = 1, ni
366  do j = 1, ni
367  do k = 1, 6
368  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
369  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
370  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
371  enddo
372  enddo
373  enddo
374  edstrain(5,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
375  edstress(5,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
376  tdstrain(5,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
377  edstrain(6,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
378  edstress(6,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
379  tdstrain(6,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
380  edstrain(7,1:6) = ( edstrain(3,1:6) + edstrain(1,1:6) ) / 2.0
381  edstress(7,1:6) = ( edstress(3,1:6) + edstress(1,1:6) ) / 2.0
382  tdstrain(7,1:6) = ( tdstrain(3,1:6) + tdstrain(1,1:6) ) / 2.0
383  edstrain(8,1:6) = ( edstrain(1,1:6) + edstrain(4,1:6) ) / 2.0
384  edstress(8,1:6) = ( edstress(1,1:6) + edstress(4,1:6) ) / 2.0
385  tdstrain(8,1:6) = ( tdstrain(1,1:6) + tdstrain(4,1:6) ) / 2.0
386  edstrain(9,1:6) = ( edstrain(2,1:6) + edstrain(4,1:6) ) / 2.0
387  edstress(9,1:6) = ( edstress(2,1:6) + edstress(4,1:6) ) / 2.0
388  tdstrain(9,1:6) = ( tdstrain(2,1:6) + tdstrain(4,1:6) ) / 2.0
389  edstrain(10,1:6) = ( edstrain(3,1:6) + edstrain(4,1:6) ) / 2.0
390  edstress(10,1:6) = ( edstress(3,1:6) + edstress(4,1:6) ) / 2.0
391  tdstrain(10,1:6) = ( tdstrain(3,1:6) + tdstrain(4,1:6) ) / 2.0
392  else if( etype == fe_prism15n ) then
393  do i = 1, ni
394  ic = 0
395  do j = 1, numofquadpoints(etype)
396  if( j==1 .or. j==2 .or. j==3 .or. j==7 .or. j==8 .or. j==9 ) then
397  ic = ic + 1
398  do k = 1, 6
399  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
400  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
401  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
402  enddo
403  endif
404  enddo
405  enddo
406  edstrain(7,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
407  edstress(7,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
408  tdstrain(7,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
409  edstrain(8,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
410  edstress(8,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
411  tdstrain(8,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
412  edstrain(9,1:6) = ( edstrain(3,1:6) + edstrain(1,1:6) ) / 2.0
413  edstress(9,1:6) = ( edstress(3,1:6) + edstress(1,1:6) ) / 2.0
414  tdstrain(9,1:6) = ( tdstrain(3,1:6) + tdstrain(1,1:6) ) / 2.0
415  edstrain(10,1:6) = ( edstrain(4,1:6) + edstrain(5,1:6) ) / 2.0
416  edstress(10,1:6) = ( edstress(4,1:6) + edstress(5,1:6) ) / 2.0
417  tdstrain(10,1:6) = ( tdstrain(4,1:6) + tdstrain(5,1:6) ) / 2.0
418  edstrain(11,1:6) = ( edstrain(5,1:6) + edstrain(6,1:6) ) / 2.0
419  edstress(11,1:6) = ( edstress(5,1:6) + edstress(6,1:6) ) / 2.0
420  tdstrain(11,1:6) = ( tdstrain(5,1:6) + tdstrain(6,1:6) ) / 2.0
421  edstrain(12,1:6) = ( edstrain(6,1:6) + edstrain(4,1:6) ) / 2.0
422  edstress(12,1:6) = ( edstress(6,1:6) + edstress(4,1:6) ) / 2.0
423  tdstrain(12,1:6) = ( tdstrain(6,1:6) + tdstrain(4,1:6) ) / 2.0
424  edstrain(13,1:6) = ( edstrain(1,1:6) + edstrain(4,1:6) ) / 2.0
425  edstress(13,1:6) = ( edstress(1,1:6) + edstress(4,1:6) ) / 2.0
426  tdstrain(13,1:6) = ( tdstrain(1,1:6) + tdstrain(4,1:6) ) / 2.0
427  edstrain(14,1:6) = ( edstrain(2,1:6) + edstrain(5,1:6) ) / 2.0
428  edstress(14,1:6) = ( edstress(2,1:6) + edstress(5,1:6) ) / 2.0
429  tdstrain(14,1:6) = ( tdstrain(2,1:6) + tdstrain(5,1:6) ) / 2.0
430  edstrain(15,1:6) = ( edstrain(3,1:6) + edstrain(6,1:6) ) / 2.0
431  edstress(15,1:6) = ( edstress(3,1:6) + edstress(6,1:6) ) / 2.0
432  tdstrain(15,1:6) = ( tdstrain(3,1:6) + tdstrain(6,1:6) ) / 2.0
433  else if( etype == fe_hex20n ) then
434  do i = 1, ni
435  ic = 0
436  do j = 1, numofquadpoints(etype)
437  if( j==1 .or. j==3 .or. j==7 .or. j==9 .or. &
438  j==19 .or. j==21 .or. j==25 .or. j==27 ) then
439  ic = ic + 1
440  do k = 1, 6
441  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
442  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
443  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
444  enddo
445  endif
446  enddo
447  enddo
448  edstrain(9,1:6) = ( edstrain(1,1:6) + edstrain(2,1:6) ) / 2.0
449  edstress(9,1:6) = ( edstress(1,1:6) + edstress(2,1:6) ) / 2.0
450  tdstrain(9,1:6) = ( tdstrain(1,1:6) + tdstrain(2,1:6) ) / 2.0
451  edstrain(10,1:6) = ( edstrain(2,1:6) + edstrain(3,1:6) ) / 2.0
452  edstress(10,1:6) = ( edstress(2,1:6) + edstress(3,1:6) ) / 2.0
453  tdstrain(10,1:6) = ( tdstrain(2,1:6) + tdstrain(3,1:6) ) / 2.0
454  edstrain(11,1:6) = ( edstrain(3,1:6) + edstrain(4,1:6) ) / 2.0
455  edstress(11,1:6) = ( edstress(3,1:6) + edstress(4,1:6) ) / 2.0
456  tdstrain(11,1:6) = ( tdstrain(3,1:6) + tdstrain(4,1:6) ) / 2.0
457  edstrain(12,1:6) = ( edstrain(4,1:6) + edstrain(1,1:6) ) / 2.0
458  edstress(12,1:6) = ( edstress(4,1:6) + edstress(1,1:6) ) / 2.0
459  tdstrain(12,1:6) = ( tdstrain(4,1:6) + tdstrain(1,1:6) ) / 2.0
460  edstrain(13,1:6) = ( edstrain(5,1:6) + edstrain(6,1:6) ) / 2.0
461  edstress(13,1:6) = ( edstress(5,1:6) + edstress(6,1:6) ) / 2.0
462  tdstrain(13,1:6) = ( tdstrain(5,1:6) + tdstrain(6,1:6) ) / 2.0
463  edstrain(14,1:6) = ( edstrain(6,1:6) + edstrain(7,1:6) ) / 2.0
464  edstress(14,1:6) = ( edstress(6,1:6) + edstress(7,1:6) ) / 2.0
465  tdstrain(14,1:6) = ( tdstrain(6,1:6) + tdstrain(7,1:6) ) / 2.0
466  edstrain(15,1:6) = ( edstrain(7,1:6) + edstrain(8,1:6) ) / 2.0
467  edstress(15,1:6) = ( edstress(7,1:6) + edstress(8,1:6) ) / 2.0
468  tdstrain(15,1:6) = ( tdstrain(7,1:6) + tdstrain(8,1:6) ) / 2.0
469  edstrain(16,1:6) = ( edstrain(8,1:6) + edstrain(5,1:6) ) / 2.0
470  edstress(16,1:6) = ( edstress(8,1:6) + edstress(5,1:6) ) / 2.0
471  tdstrain(16,1:6) = ( tdstrain(8,1:6) + tdstrain(5,1:6) ) / 2.0
472  edstrain(17,1:6) = ( edstrain(1,1:6) + edstrain(5,1:6) ) / 2.0
473  edstress(17,1:6) = ( edstress(1,1:6) + edstress(5,1:6) ) / 2.0
474  tdstrain(17,1:6) = ( tdstrain(1,1:6) + tdstrain(5,1:6) ) / 2.0
475  edstrain(18,1:6) = ( edstrain(2,1:6) + edstrain(6,1:6) ) / 2.0
476  edstress(18,1:6) = ( edstress(2,1:6) + edstress(6,1:6) ) / 2.0
477  tdstrain(18,1:6) = ( tdstrain(2,1:6) + tdstrain(6,1:6) ) / 2.0
478  edstrain(19,1:6) = ( edstrain(3,1:6) + edstrain(7,1:6) ) / 2.0
479  edstress(19,1:6) = ( edstress(3,1:6) + edstress(7,1:6) ) / 2.0
480  tdstrain(19,1:6) = ( tdstrain(3,1:6) + tdstrain(7,1:6) ) / 2.0
481  edstrain(20,1:6) = ( edstrain(4,1:6) + edstrain(8,1:6) ) / 2.0
482  edstress(20,1:6) = ( edstress(4,1:6) + edstress(8,1:6) ) / 2.0
483  tdstrain(20,1:6) = ( tdstrain(4,1:6) + tdstrain(8,1:6) ) / 2.0
484  endif
485  end subroutine nodalstress_inv3
486 
487  function get_mises(s)
488  use m_fstr
489  implicit none
490  real(kind=kreal) :: get_mises, s(1:6)
491  real(kind=kreal) :: s11, s22, s33, s12, s23, s13, ps, smises
492 
493  s11 = s(1)
494  s22 = s(2)
495  s33 = s(3)
496  s12 = s(4)
497  s23 = s(5)
498  s13 = s(6)
499  ps = ( s11 + s22 + s33 ) / 3.0d0
500  smises = 0.5d0 * ( (s11-ps)**2 + (s22-ps)**2 + (s33-ps)**2 ) + s12**2 + s23**2 + s13**2
501  get_mises = dsqrt( 3.0d0 * smises )
502 
503  end function get_mises
504 
506  !----------------------------------------------------------------------*
507  subroutine fstr_nodalstress2d( hecMESH, fstrSOLID )
508  !----------------------------------------------------------------------*
509  use m_fstr
510  use m_static_lib
511  type (hecmwST_local_mesh) :: hecMESH
512  type (fstr_solid) :: fstrSOLID
513  real(kind=kreal), pointer :: tnstrain(:), testrain(:)
514  !C** local variables
515  integer(kind=kint) :: itype, icel, ic, is, iE, jS, i, j, ic_type, nn, ni, ID_area
516  real(kind=kreal) :: estrain(4), estress(4), tstrain(4), naturalcoord(4)
517  real(kind=kreal) :: edstrain(8,4), edstress(8,4), tdstrain(8,4)
518  real(kind=kreal) :: s11, s22, s33, s12, s23, s13, ps, smises
519  real(kind=kreal), allocatable :: func(:,:), inv_func(:,:)
520  integer(kind=kint), allocatable :: nnumber(:)
521 
522  tnstrain => fstrsolid%tnstrain
523  testrain => fstrsolid%testrain
524  call fstr_solid_phys_clear(fstrsolid)
525 
526  allocate( nnumber(hecmesh%n_node) )
527  if( .not. associated(fstrsolid%is_rot) ) allocate( fstrsolid%is_rot(hecmesh%n_node) )
528  nnumber = 0
529  fstrsolid%is_rot = 0
530 
531  !C +-------------------------------+
532  !C | according to ELEMENT TYPE |
533  !C +-------------------------------+
534  do itype = 1, hecmesh%n_elem_type
535  is = hecmesh%elem_type_index(itype-1) + 1
536  ie = hecmesh%elem_type_index(itype )
537  ic_type = hecmesh%elem_type_item(itype)
538  if( .not. hecmw_is_etype_surface(ic_type) ) cycle
539  !C** set number of nodes and shape function
540  nn = hecmw_get_max_node( ic_type )
541  ni = numofquadpoints( ic_type )
542  allocate( func(ni,nn), inv_func(nn,ni) )
543  if( ic_type == fe_tri6n ) then
544  ic = hecmw_get_max_node( fe_tri3n )
545  do i = 1, ni
546  call getquadpoint( ic_type, i, naturalcoord )
547  call getshapefunc( fe_tri3n, naturalcoord, func(i,1:ic) )
548  enddo
549  call inverse_func( ic, func, inv_func )
550  else if( ic_type == fe_quad4n ) then
551  do i = 1, ni
552  call getquadpoint( ic_type, i, naturalcoord )
553  call getshapefunc( ic_type, naturalcoord, func(i,1:nn) )
554  enddo
555  call inverse_func( ni, func, inv_func )
556  else if( ic_type == fe_quad8n ) then
557  ic = 0
558  do i = 1, ni
559  if( i==1 .or. i==3 .or. i==7 .or. i==9 ) then
560  ic = ic + 1
561  call getquadpoint( ic_type, i, naturalcoord )
562  call getshapefunc( fe_quad4n, naturalcoord, func(ic,1:4) )
563  endif
564  enddo
565  call inverse_func( ic, func, inv_func )
566  ni = ic
567  endif
568  !C** element loop
569  do icel = is, ie
570  js = hecmesh%elem_node_index(icel-1)
571  id_area = hecmesh%elem_ID(icel*2)
572  !--- calculate nodal stress and strain
573  if( ic_type == fe_tri6n .or. ic_type == fe_quad4n .or. ic_type == fe_quad8n ) then
574  call nodalstress_inv2( ic_type, ni, fstrsolid%elements(icel)%gausses, &
575  inv_func, edstrain(1:nn,1:4), edstress(1:nn,1:4), &
576  tdstrain(1:nn,1:4) )
577  else
578  call nodalstress_c2( ic_type, nn, fstrsolid%elements(icel)%gausses, &
579  edstrain(1:nn,1:4), edstress(1:nn,1:4) )
580  ! call NodalStress_C2( ic_type, nn, fstrSOLID%elements(icel)%gausses, &
581  ! edstrain(1:nn,1:4), edstress(1:nn,1:4), tdstrain(1:nn,1:4) )
582  endif
583  do j = 1, nn
584  ic = hecmesh%elem_node_item(js+j)
585  fstrsolid%STRAIN(3*ic-2) = fstrsolid%STRAIN(3*ic-2) + edstrain(j,1)
586  fstrsolid%STRAIN(3*ic-1) = fstrsolid%STRAIN(3*ic-1) + edstrain(j,2)
587  fstrsolid%STRAIN(3*ic-0) = fstrsolid%STRAIN(3*ic-0) + edstrain(j,3)
588  fstrsolid%STRESS(3*ic-2) = fstrsolid%STRESS(3*ic-2) + edstress(j,1)
589  fstrsolid%STRESS(3*ic-1) = fstrsolid%STRESS(3*ic-1) + edstress(j,2)
590  fstrsolid%STRESS(3*ic-0) = fstrsolid%STRESS(3*ic-0) + edstress(j,3)
591 
592  if( associated(tnstrain) ) then
593  tnstrain(3*ic-2) = tnstrain(3*ic-2) + tdstrain(j,1)
594  tnstrain(3*ic-1) = tnstrain(3*ic-1) + tdstrain(j,2)
595  tnstrain(3*ic ) = tnstrain(3*ic ) + tdstrain(j,3)
596  endif
597  nnumber(ic) = nnumber(ic) + 1
598  enddo
599  !--- calculate elemental stress and strain
600  ! if( ID_area == hecMESH%my_rank ) then
601  call elementstress_c2( ic_type, fstrsolid%elements(icel)%gausses, estrain, estress )
602  ! call ElementStress_C2( ic_type, fstrSOLID%elements(icel)%gausses, estrain, estress, tstrain )
603 
604  fstrsolid%ESTRAIN(3*icel-2) = estrain(1)
605  fstrsolid%ESTRAIN(3*icel-1) = estrain(2)
606  fstrsolid%ESTRAIN(3*icel-0) = estrain(3)
607  fstrsolid%ESTRESS(3*icel-2) = estress(1)
608  fstrsolid%ESTRESS(3*icel-1) = estress(2)
609  fstrsolid%ESTRESS(3*icel-0) = estress(3)
610 
611  !if( associated(testrain) ) then
612  ! testrain(3*icel-2) = tstrain(1)
613  ! testrain(3*icel-1) = tstrain(2)
614  ! testrain(3*icel ) = tstrain(3)
615  !endif
616  s11 = estress(1)
617  s22 = estress(2)
618  s12 = estress(3)
619  smises = 0.5d0 * ((s11-s22)**2+(s11)**2+(s22)**2) + 3*s12**2
620  fstrsolid%EMISES(icel) = sqrt( smises )
621  ! endif
622  enddo
623  deallocate( func, inv_func )
624  enddo
625 
626  !C** average over nodes
627  do i = 1, hecmesh%n_node
628  if( nnumber(i) == 0 ) cycle
629  fstrsolid%STRAIN(3*i-2:3*i-0) = fstrsolid%STRAIN(3*i-2:3*i-0) / nnumber(i)
630  fstrsolid%STRESS(3*i-2:3*i-0) = fstrsolid%STRESS(3*i-2:3*i-0) / nnumber(i)
631  if( associated(tnstrain) ) tnstrain(3*i-2:3*i) = tnstrain(3*i-2:3*i) / nnumber(i)
632  enddo
633  !C** calculate von MISES stress
634  do i = 1, hecmesh%n_node
635  s11 = fstrsolid%STRESS(3*i-2)
636  s22 = fstrsolid%STRESS(3*i-1)
637  s12 = fstrsolid%STRESS(3*i-0)
638  smises = 0.5d0 * ((s11-s22)**2+(s11)**2+(s22)**2) + 3*s12**2
639  fstrsolid%MISES(i) = sqrt( smises )
640  enddo
641 
642  deallocate( nnumber )
643  end subroutine fstr_nodalstress2d
644 
645  !----------------------------------------------------------------------*
646  subroutine nodalstress_inv2( etype, ni, gausses, func, edstrain, edstress, tdstrain )
647  !----------------------------------------------------------------------*
648  use m_fstr
649  use mmechgauss
650  integer(kind=kint) :: etype, ni
651  type(tgaussstatus) :: gausses(:)
652  real(kind=kreal) :: func(:,:), edstrain(:,:), edstress(:,:), tdstrain(:,:)
653  integer :: i, j, k, ic
654 
655  edstrain = 0.0d0
656  edstress = 0.0d0
657  tdstrain = 0.0d0
658 
659  if( etype == fe_quad4n ) then
660  do i = 1, ni
661  do j = 1, ni
662  do k = 1, 4
663  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
664  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
665  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
666  enddo
667  enddo
668  enddo
669  else if( etype == fe_tri6n ) then
670  do i = 1, ni
671  do j = 1, ni
672  do k = 1, 4
673  edstrain(i,k) = edstrain(i,k) + func(i,j) * gausses(j)%strain_out(k)
674  edstress(i,k) = edstress(i,k) + func(i,j) * gausses(j)%stress_out(k)
675  ! tdstrain(i,k) = tdstrain(i,k) + func(i,j) * gausses(j)%tstrain(k)
676  enddo
677  enddo
678  enddo
679  edstrain(4,1:4) = ( edstrain(1,1:4) + edstrain(2,1:4) ) / 2.0
680  edstress(4,1:4) = ( edstress(1,1:4) + edstress(2,1:4) ) / 2.0
681  tdstrain(4,1:4) = ( tdstrain(1,1:4) + tdstrain(2,1:4) ) / 2.0
682  edstrain(5,1:4) = ( edstrain(2,1:4) + edstrain(3,1:4) ) / 2.0
683  edstress(5,1:4) = ( edstress(2,1:4) + edstress(3,1:4) ) / 2.0
684  tdstrain(5,1:4) = ( tdstrain(2,1:4) + tdstrain(3,1:4) ) / 2.0
685  edstrain(6,1:4) = ( edstrain(3,1:4) + edstrain(1,1:4) ) / 2.0
686  edstress(6,1:4) = ( edstress(3,1:4) + edstress(1,1:4) ) / 2.0
687  tdstrain(6,1:4) = ( tdstrain(3,1:4) + tdstrain(1,1:4) ) / 2.0
688  else if( etype == fe_quad8n ) then
689  do i = 1, ni
690  ic = 0
691  do j = 1, numofquadpoints(etype)
692  if( j==1 .or. j==3 .or. j==7 .or. j==9 ) then
693  ic = ic + 1
694  do k = 1, 4
695  edstrain(i,k) = edstrain(i,k) + func(i,ic) * gausses(j)%strain_out(k)
696  edstress(i,k) = edstress(i,k) + func(i,ic) * gausses(j)%stress_out(k)
697  ! tdstrain(i,k) = tdstrain(i,k) + func(i,ic) * gausses(j)%tstrain(k)
698  enddo
699  endif
700  enddo
701  enddo
702  edstrain(5,1:4) = ( edstrain(1,1:4) + edstrain(2,1:4) ) / 2.0
703  edstress(5,1:4) = ( edstress(1,1:4) + edstress(2,1:4) ) / 2.0
704  tdstrain(5,1:4) = ( tdstrain(1,1:4) + tdstrain(2,1:4) ) / 2.0
705  edstrain(6,1:4) = ( edstrain(2,1:4) + edstrain(3,1:4) ) / 2.0
706  edstress(6,1:4) = ( edstress(2,1:4) + edstress(3,1:4) ) / 2.0
707  tdstrain(6,1:4) = ( tdstrain(2,1:4) + tdstrain(3,1:4) ) / 2.0
708  edstrain(7,1:4) = ( edstrain(3,1:4) + edstrain(4,1:4) ) / 2.0
709  edstress(7,1:4) = ( edstress(3,1:4) + edstress(4,1:4) ) / 2.0
710  tdstrain(7,1:4) = ( tdstrain(3,1:4) + tdstrain(4,1:4) ) / 2.0
711  edstrain(8,1:4) = ( edstrain(4,1:4) + edstrain(1,1:4) ) / 2.0
712  edstress(8,1:4) = ( edstress(4,1:4) + edstress(1,1:4) ) / 2.0
713  tdstrain(8,1:4) = ( tdstrain(4,1:4) + tdstrain(1,1:4) ) / 2.0
714  endif
715  end subroutine nodalstress_inv2
716 
717  !----------------------------------------------------------------------*
718  subroutine inverse_func( n, a, inv_a )
719  !----------------------------------------------------------------------*
720  use m_fstr
721  integer(kind=kint) :: n
722  real(kind=kreal) :: a(:,:), inv_a(:,:)
723  integer(kind=kint) :: i, j, k
724  real(kind=kreal) :: buf
725 
726  do i = 1, n
727  do j = 1, n
728  if( i == j ) then
729  inv_a(i,j) = 1.0
730  else
731  inv_a(i,j) = 0.0
732  endif
733  enddo
734  enddo
735 
736  do i = 1, n
737  buf = 1.0 / a(i,i)
738  do j = 1, n
739  a(i,j) = a(i,j) * buf
740  inv_a(i,j) = inv_a(i,j) *buf
741  enddo
742  do j = 1, n
743  if( i /= j ) then
744  buf = a(j,i)
745  do k = 1, n
746  a(j,k) = a(j,k) - a(i,k) * buf
747  inv_a(j,k) = inv_a(j,k) - inv_a(i,k) * buf
748  enddo
749  endif
750  enddo
751  enddo
752  end subroutine inverse_func
753 
755  !----------------------------------------------------------------------*
756  subroutine fstr_nodalstress6d( hecMESH, fstrSOLID )
757  !----------------------------------------------------------------------*
758  use m_fstr
759  use m_static_lib
760  type (hecmwST_local_mesh) :: hecMESH
761  type (fstr_solid) :: fstrSOLID
762  !C** local variables
763  integer(kind=kint) :: itype, icel, is, iE, jS, i, j, k, it, ic, ic_type, nn, isect, ihead, ID_area
764  integer(kind=kint) :: nodLOCAL(20), n_layer, ntot_lyr, nlyr, n_totlyr, com_total_layer, shellmatl
765  real(kind=kreal) :: ecoord(3,9), edisp(6,9), estrain(6), estress(6), ndstrain(9,6), ndstress(9,6)
766  real(kind=kreal) :: thick, thick_layer
767  real(kind=kreal) :: s11, s22, s33, s12, s23, s13, t11, t22, t33, t12, t23, t13, ps, smises, tmises
768  integer(kind=kint), allocatable :: nnumber(:)
769  type(fstr_solid_physic_val), pointer :: layer => null()
770 
771  call fstr_solid_phys_clear(fstrsolid)
772 
773  n_totlyr = fstrsolid%max_lyr
774 
775  allocate( nnumber(hecmesh%n_node) )
776  if( .not. associated(fstrsolid%is_rot) ) allocate( fstrsolid%is_rot(hecmesh%n_node) )
777  nnumber = 0
778  fstrsolid%is_rot = 0
779 
780  !C +-------------------------------+
781  !C | according to ELEMENT TYPE |
782  !C +-------------------------------+
783  do itype = 1, hecmesh%n_elem_type
784  is = hecmesh%elem_type_index(itype-1) + 1
785  ie = hecmesh%elem_type_index(itype )
786  ic_type = hecmesh%elem_type_item(itype)
787  if( .not. hecmw_is_etype_shell(ic_type) ) then
788  ntot_lyr = 0
789  cycle
790  end if
791  nn = hecmw_get_max_node( ic_type )
792  !C** element loop
793  do icel = is, ie
794  js = hecmesh%elem_node_index(icel-1)
795  id_area = hecmesh%elem_ID(icel*2)
796  do j = 1, nn
797  nodlocal(j) = hecmesh%elem_node_item(js+j)
798  ecoord(1,j) = hecmesh%node(3*nodlocal(j)-2)
799  ecoord(2,j) = hecmesh%node(3*nodlocal(j)-1)
800  ecoord(3,j) = hecmesh%node(3*nodlocal(j) )
801  edisp(1,j) = fstrsolid%unode(6*nodlocal(j)-5)
802  edisp(2,j) = fstrsolid%unode(6*nodlocal(j)-4)
803  edisp(3,j) = fstrsolid%unode(6*nodlocal(j)-3)
804  edisp(4,j) = fstrsolid%unode(6*nodlocal(j)-2)
805  edisp(5,j) = fstrsolid%unode(6*nodlocal(j)-1)
806  edisp(6,j) = fstrsolid%unode(6*nodlocal(j) )
807  enddo
808  isect = hecmesh%section_ID(icel)
809  ihead = hecmesh%section%sect_R_index(isect-1)
810  thick = hecmesh%section%sect_R_item(ihead+1)
811  !--- calculate elemental stress and strain
812  if( ic_type == 731 .or. ic_type == 741 .or. ic_type == 743 ) then
813  ntot_lyr = fstrsolid%elements(icel)%gausses(1)%pMaterial%totallyr
814  do nlyr=1,ntot_lyr
815  call elementstress_shell_mitc( ic_type, nn, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
816  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), thick, 1.0d0, nlyr, ntot_lyr)
817  do j = 1, nn
818  i = nodlocal(j)
819  layer => fstrsolid%SHELL%LAYER(nlyr)%PLUS
820  do k = 1, 6
821  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + ndstrain(j,k)
822  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + ndstress(j,k)
823  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + ndstrain(j,k)/nn
824  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + ndstress(j,k)/nn
825  enddo
826  enddo
827  !minus section
828  call elementstress_shell_mitc( ic_type, nn, 6, ecoord, fstrsolid%elements(icel)%gausses, edisp, &
829  & ndstrain(1:nn,1:6), ndstress(1:nn,1:6), thick,-1.0d0, nlyr, ntot_lyr)
830  do j = 1, nn
831  i = nodlocal(j)
832  layer => fstrsolid%SHELL%LAYER(nlyr)%MINUS
833  do k = 1, 6
834  layer%STRAIN(6*(i-1)+k) = layer%STRAIN(6*(i-1)+k) + ndstrain(j,k)
835  layer%STRESS(6*(i-1)+k) = layer%STRESS(6*(i-1)+k) + ndstress(j,k)
836  layer%ESTRAIN(6*(icel-1)+k) = layer%ESTRAIN(6*(icel-1)+k) + ndstrain(j,k)/nn
837  layer%ESTRESS(6*(icel-1)+k) = layer%ESTRESS(6*(icel-1)+k) + ndstress(j,k)/nn
838  enddo
839  enddo
840  enddo
841  call fstr_getavg_shell(nn,fstrsolid,icel,nodlocal,ndstrain(1:nn,1:6),ndstress(1:nn,1:6),estrain,estress)
842  endif
843 
844  !if( ID_area == hecMESH%my_rank ) then
845  !ADD VALUE and Count node
846  do j = 1, nn
847  ic = hecmesh%elem_node_item(js+j)
848  fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRAIN(6*(ic-1)+1:6*(ic-1)+6) + ndstrain(j,1:6)
849  fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) = fstrsolid%STRESS(6*(ic-1)+1:6*(ic-1)+6) + ndstress(j,1:6)
850  !if( associated(tnstrain) )then
851  ! tnstrain(6*(ic-1)+1:6*(ic-1)+6) = tnstrain(6*(ic-1)+1:6*(ic-1)+6) + tdstrain(j,1:6)
852  !endif
853  nnumber(ic) = nnumber(ic) + 1
854  enddo
855 
856  fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRAIN(6*(icel-1)+1:6*(icel-1)+6) + estrain(1:6)
857  fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) = fstrsolid%ESTRESS(6*(icel-1)+1:6*(icel-1)+6) + estress(1:6)
858  !endif
859  enddo
860  enddo
861 
862  !C** calculate nodal stress and strain
863  do i = 1, hecmesh%n_node
864  if( nnumber(i) == 0 ) cycle
865  fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
866  fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) = fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
867  !if( associated(tnstrain) )then
868  ! tnstrain(6*(i-1)+1:6*(i-1)+6) = tnstrain(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
869  !endif
870  enddo
871 
872  do nlyr = 1, ntot_lyr
873  do i = 1, hecmesh%n_node
874  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
875  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
876  fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
877  & fstrsolid%SHELL%LAYER(nlyr)%PLUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
878  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) = &
879  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRAIN(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
880  fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) = &
881  & fstrsolid%SHELL%LAYER(nlyr)%MINUS%STRESS(6*(i-1)+1:6*(i-1)+6) / nnumber(i)
882  enddo
883  enddo
884 
885  !C** calculate von MISES stress
886  do i = 1, hecmesh%n_node
887  fstrsolid%MISES(i) = get_mises(fstrsolid%STRESS(6*(i-1)+1:6*(i-1)+6))
888  enddo
889  do i = 1, hecmesh%n_elem
890  fstrsolid%EMISES(i) = get_mises(fstrsolid%ESTRESS(6*(i-1)+1:6*(i-1)+6))
891  enddo
892  deallocate( nnumber )
893 
894  end subroutine fstr_nodalstress6d
895 
896  subroutine make_principal(fstrSOLID, hecMESH, RES)
897  use hecmw_util
898  use m_fstr
899  use m_out
900  use m_static_lib
901 
902  type(fstr_solid) :: fstrSOLID
903  type(hecmwst_local_mesh) :: hecMESH
904  type(fstr_solid_physic_val) :: RES
905  integer(kind=kint) :: i, flag
906  real(kind=kreal) :: tmat(3, 3), tvec(3), strain(6)
907 
908  flag=ieor(flag,flag)
909  if( fstrsolid%output_ctrl(3)%outinfo%on(19) .or. fstrsolid%output_ctrl(4)%outinfo%on(19) ) then
910  if ( .not. associated(res%PSTRESS) ) then
911  allocate(res%PSTRESS( 3*hecmesh%n_node ))
912  endif
913  flag=ior(flag,b'00000001')
914  end if
915  if( fstrsolid%output_ctrl(3)%outinfo%on(23) .or. fstrsolid%output_ctrl(4)%outinfo%on(23) ) then
916  if ( .not. associated(res%PSTRESS_VECT) ) then
917  allocate(res%PSTRESS_VECT( 3*hecmesh%n_node ,3))
918  endif
919  flag=ior(flag,b'00000010')
920  end if
921  if( fstrsolid%output_ctrl(3)%outinfo%on(21) .or. fstrsolid%output_ctrl(4)%outinfo%on(21) ) then
922  if ( .not. associated(res%PSTRAIN) ) then
923  allocate(res%PSTRAIN( 3*hecmesh%n_node ))
924  endif
925  flag=ior(flag,b'00000100')
926  end if
927  if( fstrsolid%output_ctrl(3)%outinfo%on(25) .or. fstrsolid%output_ctrl(4)%outinfo%on(25) ) then
928  if ( .not. associated(res%PSTRAIN_VECT) ) then
929  allocate(res%PSTRAIN_VECT( 3*hecmesh%n_node ,3))
930  endif
931  flag=ior(flag,b'00001000')
932  end if
933  if( fstrsolid%output_ctrl(3)%outinfo%on(20) .or. fstrsolid%output_ctrl(4)%outinfo%on(20) ) then
934  if ( .not. associated(res%EPSTRESS) ) then
935  allocate(res%EPSTRESS( 3*hecmesh%n_elem ))
936  endif
937  flag=ior(flag,b'00010000')
938  end if
939  if( fstrsolid%output_ctrl(3)%outinfo%on(24) .or. fstrsolid%output_ctrl(4)%outinfo%on(24) ) then
940  if ( .not. associated(res%EPSTRESS_VECT) ) then
941  allocate(res%EPSTRESS_VECT( 3*hecmesh%n_elem ,3))
942  endif
943  flag=ior(flag,b'00100000')
944  end if
945  if( fstrsolid%output_ctrl(3)%outinfo%on(22) .or. fstrsolid%output_ctrl(4)%outinfo%on(22) ) then
946  if ( .not. associated(res%EPSTRAIN) ) then
947  allocate(res%EPSTRAIN( 3*hecmesh%n_elem ))
948  endif
949  flag=ior(flag,b'01000000')
950  end if
951  if( fstrsolid%output_ctrl(3)%outinfo%on(26) .or. fstrsolid%output_ctrl(4)%outinfo%on(26) ) then
952  if ( .not. associated(res%EPSTRAIN_VECT) ) then
953  allocate(res%EPSTRAIN_VECT( 3*hecmesh%n_elem ,3))
954  endif
955  flag=ior(flag,b'10000000')
956  end if
957 
958  if (iand(flag,b'00000011') /= 0) then
959  do i = 1, hecmesh%n_node
960  call get_principal(res%STRESS(6*i-5:6*i), tvec, tmat)
961  if (iand(flag,b'00000001') /= 0) res%PSTRESS(3*(i-1)+1:3*(i-1)+3)=tvec
962  if (iand(flag,b'00000010') /= 0) res%PSTRESS_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
963  end do
964  end if
965  if (iand(flag,b'00001100') /= 0) then
966  do i = 1, hecmesh%n_node
967  strain(1:6) = res%STRAIN(6*i-5:6*i)
968  strain(4:6) = 0.5d0*strain(4:6)
969  call get_principal(strain, tvec, tmat)
970  if (iand(flag,b'00000100') /= 0) res%PSTRAIN(3*(i-1)+1:3*(i-1)+3)=tvec
971  if (iand(flag,b'00001000') /= 0) res%PSTRAIN_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
972  end do
973  end if
974 
975  if (iand(flag,b'00110000') /= 0) then
976  do i = 1, hecmesh%n_elem
977  call get_principal( res%ESTRESS(6*i-5:6*i), tvec, tmat)
978  if (iand(flag,b'00010000') /= 0) res%EPSTRESS(3*(i-1)+1:3*(i-1)+3)=tvec
979  if (iand(flag,b'00100000') /= 0) res%EPSTRESS_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
980  end do
981  end if
982  if (iand(flag,b'11000000') /= 0) then
983  do i = 1, hecmesh%n_elem
984  strain(1:6) = res%ESTRAIN(6*i-5:6*i)
985  strain(4:6) = 0.5d0*strain(4:6)
986  call get_principal(strain, tvec, tmat)
987  if (iand(flag,b'01000000') /= 0) res%EPSTRAIN(3*(i-1)+1:3*(i-1)+3)=tvec
988  if (iand(flag,b'10000000') /= 0) res%EPSTRAIN_VECT(3*(i-1)+1:3*(i-1)+3,1:3)=tmat
989  end do
990  end if
991  end subroutine make_principal
992 
993 end module m_fstr_nodalstress
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal
This module provides functions to caluclation nodal stress.
subroutine fstr_stress_add_shelllyr(nn, fstrSOLID, icel, nodLOCAL, nlyr, strain, stress, flag)
subroutine fstr_nodalstress3d(hecMESH, fstrSOLID)
Calculate NODAL STRESS of solid elements.
subroutine fstr_nodalstress6d(hecMESH, fstrSOLID)
Calculate NODAL STRESS of shell elements.
real(kind=kreal) function get_mises(s)
subroutine fstr_nodalstress2d(hecMESH, fstrSOLID)
Calculate NODAL STRESS of plane elements.
subroutine make_principal(fstrSOLID, hecMESH, RES)
subroutine fstr_getavg_shell(nn, fstrSOLID, icel, nodLOCAL, strain, stress, estrain, estress)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
subroutine fstr_solid_phys_clear(fstrSOLID)
Definition: m_fstr.f90:1094
This module manages step infomation.
Definition: m_out.f90:6
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
Data for STATIC ANSLYSIS (fstrSOLID)
Definition: m_fstr.f90:193
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13