FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
9contains
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
993end 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 make_principal(fstrsolid, hecmesh, res)
subroutine fstr_getavg_shell(nn, fstrsolid, icel, nodlocal, strain, stress, estrain, estress)
subroutine fstr_stress_add_shelllyr(nn, fstrsolid, icel, nodlocal, nlyr, strain, stress, flag)
subroutine fstr_nodalstress2d(hecmesh, fstrsolid)
Calculate NODAL STRESS of plane elements.
subroutine fstr_nodalstress6d(hecmesh, fstrsolid)
Calculate NODAL STRESS of shell elements.
real(kind=kreal) function get_mises(s)
subroutine fstr_nodalstress3d(hecmesh, fstrsolid)
Calculate NODAL STRESS of solid elements.
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