OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
h3d_shell_tensor.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| h3d_shell_tensor ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
25!||--- called by ------------------------------------------------------
26!|| genh3d ../engine/source/output/h3d/h3d_results/genh3d.F
27!||--- calls -----------------------------------------------------
28!|| h3d_write_sh_tensor ../engine/source/output/h3d/h3d_results/h3d_write_sh_tensor.F
29!|| h3d_write_sh_tensor_array ../engine/source/output/h3d/h3d_results/h3d_write_sh_tensor_array.F
30!|| layini ../engine/source/elements/shell/coque/layini.F
31!|| roto_tens2d ../engine/source/materials/tools/roto_tens2d.F
32!|| roto_tens2d_aniso ../engine/source/materials/tools/roto_tens2d_aniso.F
33!|| sh3_tstrain ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
34!|| sh4_tstrain ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
35!|| uroto_tens2d ../engine/source/materials/tools/uroto_tens2d.F
36!|| uroto_tens2d_aniso ../engine/source/materials/tools/uroto_tens2d_aniso.F
37!||--- uses -----------------------------------------------------
38!|| drape_mod ../engine/share/modules/drape_mod.F
39!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
40!|| element_mod ../common_source/modules/elements/element_mod.F90
41!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
42!|| names_and_titles_mod ../common_source/modules/names_and_titles_mod.F
43!|| stack_mod ../engine/share/modules/stack_mod.F
44!||====================================================================
45 SUBROUTINE h3d_shell_tensor(ELBUF_TAB,SHELL_TENSOR,IPARG ,ITENS ,INVERT,NELCUT,
46 2 EL2FA ,NBF ,TENS ,EPSDOT,IADP ,
47 3 NBF_L ,NBPART,IADG ,X ,IXC ,
48 4 IGEO ,IXTG ,IPM ,STACK,ID_ELEM ,ITY_ELEM ,INFO1,
49 5 INFO2 ,IS_WRITTEN_SHELL,IPARTC ,IPARTTG ,LAYER_INPUT ,IPT_INPUT ,
50 6 PLY_INPUT,GAUSS_INPUT,IUVAR_INPUT,H3D_PART, KEYWORD,D ,
51 7 ID ,BUFMAT ,MAT_PARAM,GEO, DRAPE_SH4N, DRAPE_SH3N, DRAPEG)
52C-----------------------------------------------
53C M o d u l e s
54C-----------------------------------------------
55 USE elbufdef_mod
56 USE stack_mod
57 USE matparam_def_mod
58 USE names_and_titles_mod, ONLY: ncharline100
59 USE drape_mod
60 use element_mod , only : nixc,nixtg
61C-----------------------------------------------
62C I m p l i c i t T y p e s
63C-----------------------------------------------
64#include "implicit_f.inc"
65C-----------------------------------------------
66C C o m m o n B l o c k s
67C-----------------------------------------------
68#include "mvsiz_p.inc"
69C-----------------------------------------------
70#include "com01_c.inc"
71#include "com04_c.inc"
72#include "param_c.inc"
73#include "tabsiz_c.inc"
74C-----------------------------------------------
75C D u m m y A r g u m e n t s
76C-----------------------------------------------
77 INTEGER IPARG(NPARG,*),ITENS,INVERT(*),IUVAR_INPUT,
78 . EL2FA(*),IXC(NIXC,*), IGEO(NPROPGI,*),
79 . NELCUT,NBF,IADP(*),NBF_L,NBPART,IADG(NSPMD,*),
80 . IXTG(NIXTG,*),IPM(NPROPMI,*),ID_ELEM(*),ITY_ELEM(*),
81 . INFO1,INFO2,IS_WRITTEN_SHELL(*),IPARTC(*),IPARTTG(*),H3D_PART(*),
82 . LAYER_INPUT ,IPT_INPUT,GAUSS_INPUT,PLY_INPUT,ID
83 my_real, INTENT(IN),TARGET :: bufmat(sbufmat)
84 my_real
85 . tens(3,*),epsdot(6,*),x(3,*),shell_tensor(3,*),d(3,*)
86 my_real, INTENT(IN) :: geo(npropg,numgeo)
87 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP), TARGET :: ELBUF_TAB
88 TYPE (STACK_PLY) :: STACK
89 TYPE (MATPARAM_STRUCT_) ,DIMENSION(NUMMAT) ,INTENT(IN) :: MAT_PARAM
90 CHARACTER(LEN=NCHARLINE100):: KEYWORD
91 TYPE (DRAPE_) , INTENT(IN) :: DRAPE_SH4N(NUMELC_DRAPE), DRAPE_SH3N(NUMELTG_DRAPE)
92 TYPE (DRAPEG_), INTENT(IN) :: DRAPEG
93C-----------------------------------------------
94C L o c a l V a r i a b l e s
95C-----------------------------------------------
96 my_real :: a1,a2,a3,thk,chard,factor,factor_n,zshift
97 my_real :: value(5)
98 INTEGER I,J,K,N,NG,NEL,NFT,ITY,NPT,MPT,IPT,NBFUNCT,NCHARD,MLW,
99 . ILAY,IR,IS,IT,NPTR,NPTS,NPTT,NLAY,NPG,IPLY,IDRAPE,
100 . IPID,NS1,NS2,ISTRE,IADBUF,NUPARAM,IMAT,NNI,N0,
101 . IHBE,IREP,ISROT,IVISC,IGTYP,ISUBSTACK,
102 . id_ply,ipang,ippos,ipthk,offset,iselect,mat_orth,
103 . ixlay,ixfem,laynpt_max,numel_drape,sedrape,nlay_max,
104 . ipt_all,islice,pts,ipg,lens,mpt0
105 INTEGER NN1,NN2,NN3,NN4,NN5,NN6,NN7,NN8,NN9,NN10
106 INTEGER PID(MVSIZ),MAT(MVSIZ),IOK_PART(MVSIZ),JJ(15)
107 my_real ,DIMENSION(3,MVSIZ) :: STRAIN
108 my_real ,DIMENSION(4*MVSIZ) :: XN,YN,ZN,DXN,DYN,DZN
109 my_real ,DIMENSION(:,:) , ALLOCATABLE :: SIGE,SIGM,EPSM
110
111 TYPE(buf_lay_) ,POINTER :: BUFLY
112 TYPE(g_bufel_) ,POINTER :: GBUF
113 TYPE(l_bufel_) ,POINTER :: LBUF
114 my_real, DIMENSION(:) ,POINTER :: uparam,dir_a,dir_b
115 !
116 INTEGER, DIMENSION(:) , ALLOCATABLE :: MATLY !! MATLY(MVSIZ*LAY_MAX)
117 my_real, DIMENSION(:) , ALLOCATABLE :: thkly !! THKLY(MVSIZ*LAY_MAX*LAYNPT_MAX)
118 my_real, DIMENSION(:,:), ALLOCATABLE :: posly,thk_ly
119C-----------------------------------------------
120 ! material orthotropy defined in MATPARAM data structure
121 ! MATPARAM%ORTHOTROPY
122 ! -> 1 : ISOTROPIC
123 ! -> 2 : IRTHOTROPIC
124 ! -> 3 : ANISOTROPIC
125
126 offset = 0
127 value(1:5) = zero
128 iselect = 1
129 id_ply = 0
130 npg = 1
131
132 nn1 = 1
133 nn2 = nn1
134 nn3 = nn2
135 nn4 = nn3 + numelq
136 nn5 = nn4 + numelc
137 nn6 = nn5 + numeltg
138 nn7 = nn6
139 nn8 = nn7
140 nn9 = nn8
141 nn10= nn9
142C
143 DO i=1,numelc+numeltg
144 is_written_shell(i) = 0
145 ENDDO
146C
147 DO 490 ng=1,ngroup
148C IF(ANIM_K == 0.AND.IPARG(8,NG) == 1)GOTO 490
149 mlw = iparg(1,ng)
150 nel = iparg(2,ng)
151 nft = iparg(3,ng)
152 ity = iparg(5,ng)
153 igtyp = iparg(38,ng)
154 isrot = iparg(41,ng)
155 ixfem = iparg(54,ng)
156 isubstack = iparg(71,ng)
157 idrape = elbuf_tab(ng)%IDRAPE
158 npt = iabs(iparg(6,ng))
159 iok_part(1:nel) = 0
160!
161 DO i=1,15
162 jj(i) = nel*(i-1)
163 ENDDO
164!
165 IF (mlw /= 13) THEN
166C-----------------------------------------------
167C QUAD
168C-----------------------------------------------
169 IF(ity == 2)THEN
170 DO i=1,nel
171 n = i + nft
172 shell_tensor(1,offset+nft+i) = zero
173 shell_tensor(2,offset+nft+i) = zero
174 shell_tensor(3,offset+nft+i) = zero
175 ENDDO
176C-----------------------------------------------
177C COQUES
178C-----------------------------------------------
179 ELSEIF (ity == 3 .OR. ity == 7) THEN
180 gbuf => elbuf_tab(ng)%GBUF
181 nptr = elbuf_tab(ng)%NPTR
182 npts = elbuf_tab(ng)%NPTS
183 nlay = elbuf_tab(ng)%NLAY
184 npg = nptr*npts
185C
186 ihbe = iparg(23,ng)
187 IF (ity == 3) THEN
188 n0 = 0
189 nni = nn4
190 IF (ihbe == 11) npg = 4
191 ELSE
192 n0 = numelc
193 nni = nn5
194 IF (ihbe == 30) npg = 3 !dkt18
195 ENDIF
196c
197 IF (ity == 3) offset = 0
198 IF (ity == 7) offset = numelc
199c
200 DO i=1,nel
201 IF (ity == 3) THEN
202 id_elem(offset+nft+i) = ixc(nixc,nft+i)
203 ity_elem(offset+nft+i) = 3
204 IF( h3d_part(ipartc(nft+i)) == 1) iok_part(i) = 1
205 ELSEIF (ity == 7) THEN
206 id_elem(offset+nft+i) = ixtg(nixtg,nft+i)
207 ity_elem(offset+nft+i) = 7
208 IF( h3d_part(iparttg(nft+i)) == 1) iok_part(i) = 1
209 ENDIF
210 ENDDO
211c
212 IF (mlw == 0) GOTO 490
213C
214 a1 = zero
215 a2 = zero
216 a3 = zero
217 istre = 1
218 npt = iabs(iparg(6,ng))
219 mpt = max(1,npt)
220 mpt0 = mpt
221 IF (npt == 0) mpt = 0
222C
223 IF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16 .OR. igtyp == 17) THEN
224 npt = 1
225 mpt = npt
226 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
227 IF(layer_input == -2) THEN
228 npt= elbuf_tab(ng)%BUFLY(1)%NPTT
229 ELSEIF(layer_input == -3) THEN
230 npt= elbuf_tab(ng)%BUFLY(nlay)%NPTT
231 ELSEIF(layer_input > 0 .AND. layer_input <= nlay) THEN
232 npt= elbuf_tab(ng)%BUFLY(layer_input)%NPTT
233 ENDIF
234 IF (ply_input > 0) THEN
235 DO j=1,nlay
236 id_ply = 0
237 IF (igtyp == 51) THEN
238 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
239 ELSEIF (igtyp == 52) THEN
240 id_ply = ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
241 ENDIF
242 IF (id_ply == ply_input ) THEN
243 npt= elbuf_tab(ng)%BUFLY(j)%NPTT
244 EXIT
245 ENDIF
246 ENDDO
247 ENDIF
248 mpt = max(1,npt)
249 ENDIF
250c
251 ilay = layer_input
252 ipt = ipt_input
253 iply = ply_input
254 IF (ilay == -2) ilay = 1 ! Lower
255 IF (ilay == -3) ilay = nlay ! upper
256 IF (ipt == -2) ipt = 1 ! Lower
257 IF (igtyp == 51 .OR. igtyp == 52) THEN
258 IF (ipt == -3 .AND. ilay > 0) ipt = max(1,elbuf_tab(ng)%BUFLY(ilay)%NPTT) ! Upper
259 ELSE
260 IF (ipt == -3) ipt = max(1,npt) ! Upper
261 ENDIF
262C------------------------
263C STRESS
264C------------------------
265 IF (keyword == 'TENS/STRESS/MEMB' .OR.
266 . keyword == 'TENS/STRESS/BEND' .OR.
267 . keyword == 'TENS/STRESS' .OR.
268 . keyword == 'TENS/STRAIN' .OR.
269 . keyword == 'TENS/MSTRAIN' ) THEN
270 IF (ity == 3) THEN
271 ipid = ixc(6,nft+1)
272 DO i=1,nel
273 mat(i)=ixc(1,nft+i)
274 pid(i)=ixc(6,nft+i)
275 ENDDO
276 ELSE ! ITY == 7
277 ipid = ixtg(5,nft+1)
278 DO i=1,nel
279 mat(i)=ixtg(1,nft+i)
280 pid(i)=ixtg(5,nft+i)
281 ENDDO
282 ENDIF
283c
284 irep = igeo(6,ipid)
285 zshift = geo(199, ipid)
286 ENDIF
287 IF( keyword == 'TENS/STRAIN' .OR. keyword == 'TENS/MSTRAIN') THEN
288 !Npt_max
289 laynpt_max = 1
290 ixlay = 0
291 IF(igtyp == 51 .OR. igtyp == 52) THEN
292 DO ilay=1,nlay
293 laynpt_max = max(laynpt_max ,elbuf_tab(ng)%BUFLY(ilay)%NPTT)
294 ENDDO
295 ENDIF
296 nlay_max = max(nlay,npt)
297 ALLOCATE(matly(mvsiz*nlay_max), thkly(mvsiz*nlay_max*laynpt_max),
298 . posly(mvsiz,nlay_max*laynpt_max),thk_ly(nel,nlay_max*laynpt_max))
299 matly = 0
300 thkly = zero
301 posly = zero
302 thk_ly = zero
303c
304 ! computing position of slice or Ply
305 IF(ity == 7) THEN
306 numel_drape = numeltg_drape
307 sedrape = stdrape
308 CALL layini(
309 . elbuf_tab(ng),1 ,nel ,geo ,igeo ,
310 . mat ,pid ,thkly ,matly ,posly ,
311 . igtyp ,ixfem ,ixlay ,nlay ,mpt0 ,
312 . isubstack ,stack ,drape_sh3n ,nft ,gbuf%THK,
313 . nel ,thk_ly ,drapeg%INDX_SH3N ,sedrape,numel_drape )
314 ELSE ! ITY = 3
315 numel_drape = numelc_drape
316 sedrape = scdrape
317 CALL layini(
318 . elbuf_tab(ng),1 ,nel ,geo ,igeo ,
319 . mat ,pid ,thkly ,matly ,posly ,
320 . igtyp ,ixfem ,ixlay ,nlay ,mpt0 ,
321 . isubstack ,stack ,drape_sh4n ,nft ,gbuf%THK ,
322 . nel ,thk_ly ,drapeg%INDX_SH4N,sedrape ,numel_drape )
323 ENDIF
324 ENDIF
325C--------------------------------------------------
326 IF (keyword == 'TENS/STRESS/MEMB') THEN
327C--------------------------------------------------
328 DO i=1,nel
329 value(1) = gbuf%FOR(jj(1)+i)
330 value(2) = gbuf%FOR(jj(2)+i)
331 value(3) = gbuf%FOR(jj(3)+i)
332 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
333 . VALUE)
334 ENDDO
335C--------------------------------------------------
336 ELSEIF (keyword == 'TENS/STRESS/BEND') THEN
337C--------------------------------------------------
338 DO i=1,nel
339 value(1) = gbuf%MOM(jj(1)+i)
340 value(2) = gbuf%MOM(jj(2)+i)
341 value(3) = gbuf%MOM(jj(3)+i)
342 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
343 . VALUE)
344 ENDDO
345
346C--------------------------------------------------
347 ELSEIF (keyword == 'TENS/STRESS') THEN ! stress output in element coordinate system
348C--------------------------------------------------
349c
350 iselect = 0
351 ALLOCATE (sige(nel,3))
352 sige(1:nel,1:3) = zero
353c
354 IF (mpt == 0) THEN ! ILAYER=NULL, NPT=NULL
355 iselect = 1
356 IF (ipt_input == -2 ) THEN
357 factor = -six ! lower
358 factor_n = one + six*zshift
359 ELSEIF (ipt_input == -3) THEN
360 factor = six ! upper
361 factor_n = one - six*zshift
362 ELSE
363 factor = zero ! mem
364 factor_n = one
365 END IF
366 DO i=1,nel
367 sige(i,1) = gbuf%FOR(jj(1)+i)*factor_n + gbuf%MOM(jj(1)+i) * factor
368 sige(i,2) = gbuf%FOR(jj(2)+i)*factor_n + gbuf%MOM(jj(2)+i) * factor
369 sige(i,3) = gbuf%FOR(jj(3)+i)*factor_n + gbuf%MOM(jj(3)+i) * factor
370 ENDDO
371c
372 ELSE IF (ilay == -1 .AND. iply == -1 .AND. ipt == -1) THEN
373 iselect = 1
374 DO i=1,nel
375 sige(i,1) = gbuf%FOR(jj(1)+i)
376 sige(i,2) = gbuf%FOR(jj(2)+i)
377 sige(i,3) = gbuf%FOR(jj(3)+i)
378 ENDDO
379c
380 ELSEIF (ilay == -1 .AND. iply > 0 .AND. ipt > 0) THEN
381c /TENS/STRESS/PLY=.../NPT=...
382 DO j=1,nlay ! shells type 17,19,51 and 52 only
383 IF (igtyp == 17 .OR. igtyp == 19 .OR. igtyp == 51) THEN
384 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
385 ELSE IF (igtyp == 52) THEN
386 id_ply = ply_info(1,stack%IGEO(2+j,isubstack) - numstack)
387 END IF
388 ilay = j
389 IF (id_ply == iply .AND. ipt <= elbuf_tab(ng)%BUFLY(ilay)%NPTT) THEN
390 imat = elbuf_tab(ng)%BUFLY(ilay)%IMAT
391 ivisc = mat_param(imat)%IVISC
392 iselect = 1
393 DO i=1,nel
394 DO ir=1,nptr
395 DO is=1,npts
396 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(ir,is,ipt)
397 sige(i,1) = sige(i,1) + lbuf%SIG(jj(1) + i) / npg
398 sige(i,2) = sige(i,2) + lbuf%SIG(jj(2) + i) / npg
399 sige(i,3) = sige(i,3) + lbuf%SIG(jj(3) + i) / npg
400 ENDDO
401 ENDDO
402 ENDDO
403 IF (ivisc > 0) THEN
404 DO i=1,nel
405 DO ir=1,nptr
406 DO is=1,npts
407 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(ir,is,ipt)
408 sige(i,1) = sige(i,1) + lbuf%VISC(jj(1) + i) / npg
409 sige(i,2) = sige(i,2) + lbuf%VISC(jj(2) + i) / npg
410 sige(i,3) = sige(i,3) + lbuf%VISC(jj(3) + i) / npg
411 ENDDO
412 ENDDO
413 ENDDO
414 END IF
415 mat_orth = mat_param(imat)%ORTHOTROPY
416 IF (mat_orth > 0) THEN
417 IF (idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52) ) THEN
418 dir_a => elbuf_tab(ng)%BUFLY(ilay)%LBUF_DIR(ipt)%DIRA
419 dir_b => elbuf_tab(ng)%BUFLY(ilay)%LBUF_DIR(ipt)%DIRB
420 ELSE
421 dir_a => elbuf_tab(ng)%BUFLY(ilay)%DIRA
422 dir_b => elbuf_tab(ng)%BUFLY(ilay)%DIRB
423 ENDIF
424 END IF
425 IF (mat_orth == 2) THEN
426 CALL uroto_tens2d(nel,sige,dir_a)
427 ELSE IF (mat_orth == 3) THEN ! anisotropic (law 58,158 only)
428 CALL uroto_tens2d_aniso(nel,sige,dir_a,dir_b)
429 END IF
430 EXIT
431 ENDIF ! ID_PLY & ipt
432 ENDDO ! NLAY
433c
434 ELSEIF (ilay > 0 .AND. ilay <= nlay .AND. iply == -1) THEN
435 ! /TENS/STRESS/LAYER=
436 IF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16) THEN ! shells type 10,11,16 only
437 iselect = 1
438 imat = elbuf_tab(ng)%BUFLY(ilay)%IMAT
439 ivisc = mat_param(imat)%IVISC
440 DO i=1,nel
441 DO ir=1,nptr
442 DO is=1,npts
443 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(ir,is,1)
444 sige(i,1) = sige(i,1) + lbuf%SIG(jj(1) + i) / npg
445 sige(i,2) = sige(i,2) + lbuf%SIG(jj(2) + i) / npg
446 sige(i,3) = sige(i,3) + lbuf%SIG(jj(3) + i) / npg
447 ENDDO
448 ENDDO
449 ENDDO
450 IF (ivisc > 0) THEN
451 DO i=1,nel
452 DO ir=1,nptr
453 DO is=1,npts
454 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(ir,is,1)
455 sige(i,1) = sige(i,1) + lbuf%VISC(jj(1) + i) / npg
456 sige(i,2) = sige(i,2) + lbuf%VISC(jj(2) + i) / npg
457 sige(i,3) = sige(i,3) + lbuf%VISC(jj(3) + i) / npg
458 ENDDO
459 ENDDO
460 ENDDO
461 END IF
462 mat_orth = mat_param(imat)%ORTHOTROPY
463 IF (mat_orth > 0) THEN
464 dir_a => elbuf_tab(ng)%BUFLY(ilay)%DIRA
465 dir_b => elbuf_tab(ng)%BUFLY(ilay)%DIRB
466 END IF
467 IF (mat_orth == 2) THEN ! standard orthotropic rotation
468 CALL uroto_tens2d(nel,sige,dir_a)
469 ELSE IF (mat_orth == 3) THEN ! anisotropic (law 58,158 only)
470 CALL uroto_tens2d_aniso(nel,sige,dir_a,dir_b)
471 END IF
472 ENDIF
473
474 ELSEIF (ipt > 0 .AND. ilay ==-1 .AND. iply == -1) THEN ! shells type 1,9 only
475c /TENS/STRESS/NPT=
476 IF (igtyp == 1 .OR. igtyp == 9) THEN
477 IF (ipt <= elbuf_tab(ng)%BUFLY(1)%NPTT) THEN
478 iselect = 1
479 imat = elbuf_tab(ng)%BUFLY(1)%IMAT
480 ivisc = mat_param(imat)%IVISC
481 DO i=1,nel
482 DO ir=1,nptr
483 DO is=1,npts
484 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,ipt)
485 sige(i,1) = sige(i,1) + lbuf%SIG(jj(1) + i) / npg
486 sige(i,2) = sige(i,2) + lbuf%SIG(jj(2) + i) / npg
487 sige(i,3) = sige(i,3) + lbuf%SIG(jj(3) + i) / npg
488 ENDDO
489 ENDDO
490 ENDDO
491 IF (ivisc > 0) THEN
492 DO i=1,nel
493 DO ir=1,nptr
494 DO is=1,npts
495 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,ipt)
496 sige(i,1) = sige(i,1) + lbuf%VISC(jj(1) + i) / npg
497 sige(i,2) = sige(i,2) + lbuf%VISC(jj(2) + i) / npg
498 sige(i,3) = sige(i,3) + lbuf%VISC(jj(3) + i) / npg
499 ENDDO
500 ENDDO
501 ENDDO
502 END IF
503 mat_orth = mat_param(imat)%ORTHOTROPY
504 IF (mat_orth == 2) THEN
505 dir_a => elbuf_tab(ng)%BUFLY(1)%DIRA
506 CALL uroto_tens2d(nel,sige,dir_a)
507 END IF
508 ENDIF
509 ENDIF
510 ENDIF
511c---
512 IF (iselect == 1) THEN
513 CALL h3d_write_sh_tensor_array(
514 . iok_part ,iselect ,nel ,offset ,nft ,
515 . is_written_shell,shell_tensor,sige )
516 END IF
517c
518 DEALLOCATE (sige)
519c--------------------------------------------------
520 ELSEIF (keyword == 'TENS/MSTRESS') THEN ! stress output in material (orthotropic) coordinates
521c--------------------------------------------------
522 iselect = 0
523 ALLOCATE (sigm(nel,3))
524 sigm(1:nel,1:3) = zero
525c
526 IF (ilay == -1 .AND. iply > 0 .AND. ipt > 0) THEN
527c /TENS/MSTRESS/PLY=.../NPT=...
528 DO j=1,nlay ! shells type 17,19,51 and 52 only
529 IF (igtyp == 17 .OR. igtyp == 19 .OR. igtyp == 51) THEN
530 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
531 ELSE IF (igtyp == 52) THEN
532 id_ply = ply_info(1,stack%IGEO(2+j,isubstack) - numstack)
533 END IF
534
535 IF (id_ply == iply) THEN
536 ilay = j
537 imat = elbuf_tab(ng)%BUFLY(ilay)%IMAT
538 ivisc = mat_param(imat)%IVISC
539 IF (ipt <= elbuf_tab(ng)%BUFLY(ilay)%NPTT) THEN
540 iselect = 1
541 DO i=1,nel
542 DO ir=1,nptr
543 DO is=1,npts
544 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(ir,is,ipt)
545 sigm(i,1) = sigm(i,1) + lbuf%SIG(jj(1) + i) / npg
546 sigm(i,2) = sigm(i,2) + lbuf%SIG(jj(2) + i) / npg
547 sigm(i,3) = sigm(i,3) + lbuf%SIG(jj(3) + i) / npg
548 ENDDO
549 ENDDO
550 ENDDO
551 IF (ivisc > 0) THEN
552 DO i=1,nel
553 DO ir=1,nptr
554 DO is=1,npts
555 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(ir,is,ipt)
556 sigm(i,1) = sigm(i,1) + lbuf%VISC(jj(1) + i) / npg
557 sigm(i,2) = sigm(i,2) + lbuf%VISC(jj(2) + i) / npg
558 sigm(i,3) = sigm(i,3) + lbuf%VISC(jj(3) + i) / npg
559 ENDDO
560 ENDDO
561 ENDDO
562 END IF
563 ENDIF
564 ENDIF
565 ENDDO ! NLAY
566c
567 ELSEIF (ilay > 0 .AND. ilay <= nlay .AND. iply == -1 .AND. ipt == -1) THEN
568c /TENS/MSTRESS/LAYER=
569 IF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16) THEN ! shells type 10,11,16 only
570 iselect = 1
571 imat = elbuf_tab(ng)%BUFLY(ilay)%IMAT
572 ivisc = mat_param(imat)%IVISC
573 DO i=1,nel
574 DO ir=1,nptr
575 DO is=1,npts
576 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(ir,is,1)
577 sigm(i,1) = sigm(i,1) + lbuf%SIG(jj(1) + i) / npg
578 sigm(i,2) = sigm(i,2) + lbuf%SIG(jj(2) + i) / npg
579 sigm(i,3) = sigm(i,3) + lbuf%SIG(jj(3) + i) / npg
580 ENDDO
581 ENDDO
582 ENDDO
583 IF (ivisc > 0) THEN
584 DO i=1,nel
585 DO ir=1,nptr
586 DO is=1,npts
587 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(ir,is,ipt)
588 sigm(i,1) = sigm(i,1) + lbuf%VISC(jj(1) + i) / npg
589 sigm(i,2) = sigm(i,2) + lbuf%VISC(jj(2) + i) / npg
590 sigm(i,3) = sigm(i,3) + lbuf%VISC(jj(3) + i) / npg
591 ENDDO
592 ENDDO
593 ENDDO
594 END IF
595 ENDIF
596
597 ELSEIF (ipt > 0 .AND. ilay ==-1 .AND. iply == -1) THEN ! shells type 1,9 only
598c /TENS/MSTRESS/NPT=
599 IF (igtyp == 1 .OR. igtyp == 9) THEN
600 IF (ipt <= elbuf_tab(ng)%BUFLY(1)%NPTT) THEN
601 iselect = 1
602 imat = elbuf_tab(ng)%BUFLY(ilay)%IMAT
603 ivisc = mat_param(imat)%IVISC
604 DO i=1,nel
605 DO ir=1,nptr
606 DO is=1,npts
607 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,ipt)
608 sigm(i,1) = sigm(i,1) + lbuf%SIG(jj(1) + i) / npg
609 sigm(i,2) = sigm(i,2) + lbuf%SIG(jj(2) + i) / npg
610 sigm(i,3) = sigm(i,3) + lbuf%SIG(jj(3) + i) / npg
611 ENDDO
612 ENDDO
613 ENDDO
614 IF (ivisc > 0) THEN
615 DO i=1,nel
616 DO ir=1,nptr
617 DO is=1,npts
618 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,ipt)
619 sigm(i,1) = sigm(i,1) + lbuf%VISC(jj(1) + i) / npg
620 sigm(i,2) = sigm(i,2) + lbuf%VISC(jj(2) + i) / npg
621 sigm(i,3) = sigm(i,3) + lbuf%VISC(jj(3) + i) / npg
622 ENDDO
623 ENDDO
624 ENDDO
625 END IF
626
627 ENDIF
628 ENDIF
629 ENDIF
630c---
631 CALL h3d_write_sh_tensor_array(
632 . iok_part ,iselect ,nel ,offset ,nft ,
633 . is_written_shell,shell_tensor,sigm )
634
635 DEALLOCATE (sigm)
636C--------------------------------------------------
637 ELSE IF (keyword == 'TENS/STRAIN/MEMB') THEN
638C--------------------------------------------------
639 DO i=1,nel
640 n = i + nft
641 thk = gbuf%THK(i)
642 j = el2fa(nni+n)
643 value(1) = gbuf%STRA(jj(1)+i)
644 value(2) = gbuf%STRA(jj(2)+i)
645 value(3) = gbuf%STRA(jj(3)+i)
646 value(3) = value(3) * half
647 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
648 . shell_tensor,i,offset,nft,VALUE)
649 ENDDO
650C--------------------------------------------------
651 ELSEIF (keyword == 'TENS/STRAIN/BEND') THEN ! bend
652C--------------------------------------------------
653 DO i=1,nel
654 n = i + nft
655 thk = gbuf%THK(i)
656 j = el2fa(nni+n)
657 value(1) = gbuf%STRA(jj(6)+i) * thk
658 value(2) = gbuf%STRA(jj(7)+i) * thk
659 value(3) = gbuf%STRA(jj(8)+i) * thk
660 value(3) = value(3) * half
661 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
662 . shell_tensor,i,offset,nft,VALUE)
663 ENDDO
664C--------------------------------------------------
665 ELSEIF (keyword == 'TENS/STRAIN') THEN ! strain tensor output in element coordinate system
666C--------------------------------------------------
667 IF (mpt == 0) THEN ! ILAYER=NULL, NPT=NULL
668 iselect = 1
669 DO i=1,nel
670 IF (ipt == 1) THEN
671 factor = (zshift-half)*gbuf%THK(i)
672 ELSE
673 factor = (zshift+half)*gbuf%THK(i)
674 ENDIF
675 value(1) = gbuf%STRA(jj(1)+i) + factor*gbuf%STRA(jj(6)+i)
676 value(2) = gbuf%STRA(jj(2)+i) + factor*gbuf%STRA(jj(7)+i)
677 value(3) = gbuf%STRA(jj(3)+i) + factor*gbuf%STRA(jj(8)+i)
678 value(3) = value(3) * half
679 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
680 . shell_tensor,i,offset,nft,VALUE)
681 ENDDO
682c
683 ELSE IF (ilay == -1 .AND. iply == -1 .AND. ipt == -1) THEN
684 DO i=1,nel
685 value(1) = gbuf%STRA(jj(1)+i)
686 value(2) = gbuf%STRA(jj(2)+i)
687 value(3) = gbuf%STRA(jj(3)+i)
688 value(3) = value(3) * half
689 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
690 . shell_tensor,i,offset,nft,VALUE)
691 ENDDO
692c
693c PLY=IPLY NPT=IPT
694
695 ELSE IF (iply > 0 .AND. ipt > 0) THEN
696 IF (igtyp == 17 .OR. igtyp == 19 .OR. igtyp == 51 .OR. igtyp == 52) THEN
697 ipang = 1
698 ipthk = ipang + nlay
699 ippos = ipthk + nlay
700 ipt_all = 0
701 DO j=1,nlay
702 bufly => elbuf_tab(ng)%BUFLY(j)
703 nptt = bufly%NPTT
704 IF (igtyp == 17 .OR. igtyp == 51) THEN
705 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
706 ELSEIF (igtyp == 52) THEN
707 id_ply = ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
708 ENDIF
709 IF (id_ply == iply .AND. ipt <= nptt) THEN
710 islice = ipt_all + ipt
711 IF(npg > 1) THEN
712 lens = nel*gbuf%G_STRPG/npg
713 DO i=1,nel
714 n = i + nft
715 thk = gbuf%THK(i)
716 factor = posly(i,islice)
717 value(1:3) = zero
718 DO ipg = 1, npg
719 pts = (ipg -1)*lens
720 value(1) = value(1)+ gbuf%STRPG(pts + jj(1)+i) + factor*gbuf%STRPG(pts + jj(6)+i) * thk
721 value(2) = value(2)+ gbuf%STRPG(pts + jj(2)+i) + factor*gbuf%STRPG(pts + jj(7)+i) * thk
722 value(3) = value(3)+ gbuf%STRPG(pts + jj(3)+i) + factor*gbuf%STRPG(pts + jj(8)+i) * thk
723 ENDDO
724 value(1:3) = value(1:3)/npg
725 value(3) = value(3) * half
726 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
727 . shell_tensor,i,offset,nft,VALUE)
728 ENDDO
729 ELSE
730 DO i=1,nel
731 n = i + nft
732 thk = gbuf%THK(i)
733 factor = posly(i,islice)
734 value(1) = gbuf%STRA(jj(1)+i) + factor*gbuf%STRA(jj(6)+i) * thk
735 value(2) = gbuf%STRA(jj(2)+i) + factor*gbuf%STRA(jj(7)+i) * thk
736 value(3) = gbuf%STRA(jj(3)+i) + factor*gbuf%STRA(jj(8)+i) * thk
737 value(3) = value(3) * half
738 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
739 . shell_tensor,i,offset,nft,VALUE)
740 ENDDO
741 ENDIF ! NPG
742 ENDIF
743 ipt_all = ipt_all + nptt
744 ENDDO
745 ENDIF
746c
747 ELSEIF (ilay > 0 .AND. ilay <= nlay .AND. iply == -1) THEN
748c ILAYER=IL PLY=NULL NPT=NULL
749 IF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16) THEN
750 IF (npg > 1) THEN
751 lens = nel*gbuf%G_STRPG/npg
752 DO i=1,nel
753 n = i + nft
754 thk = gbuf%THK(i)
755 factor = posly(i,ilay)
756 j = el2fa(nni+n)
757 value(1:3) = zero
758 DO ipg = 1, npg
759 pts = (ipg -1)*lens
760 value(1) = value(1)+ gbuf%STRPG(pts + jj(1)+i) + factor*gbuf%STRPG(pts + jj(6)+i) * thk
761 value(2) = value(2)+ gbuf%STRPG(pts + jj(2)+i) + factor*gbuf%STRPG(pts + jj(7)+i) * thk
762 value(3) = value(3)+ gbuf%STRPG(pts + jj(3)+i) + factor*gbuf%STRPG(pts + jj(8)+i) * thk
763 ENDDO
764 value(1) = value(1)/npg
765 value(2) = value(2)/npg
766 value(3) = value(3)/npg
767 value(3) = value(3) * half
768 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
769 . shell_tensor,i,offset,nft,VALUE)
770 ENDDO
771 ELSE
772 DO i=1,nel
773 n = i + nft
774 thk = gbuf%THK(i)
775 factor = posly(i,ilay)
776 j = el2fa(nni+n)
777 value(1) = gbuf%STRA(jj(1)+i) + factor*gbuf%STRA(jj(6)+i) * thk
778 value(2) = gbuf%STRA(jj(2)+i) + factor*gbuf%STRA(jj(7)+i) * thk
779 value(3) = gbuf%STRA(jj(3)+i) + factor*gbuf%STRA(jj(8)+i) * thk
780 value(3) = value(3) * half
781 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
782 . shell_tensor,i,offset,nft,VALUE)
783 ENDDO
784 ENDIF ! NPG
785 ENDIF
786c
787 ELSEIF (ipt <= mpt .AND. ipt > 0) THEN
788c NPT=IPT PLY=NULL ILAYER=NULL
789 IF (igtyp == 1 .OR. igtyp == 9) THEN
790 IF(npg > 1) THEN
791 lens = nel*gbuf%G_STRPG/npg
792 DO i=1,nel
793 n = i + nft
794 thk = gbuf%THK(i)
795 factor = posly(i,ipt)
796 j = el2fa(nni+n)
797 value(1:3) = zero
798 DO ipg =1,npg
799 pts = (ipg -1)*lens
800 value(1) = value(1)+ gbuf%STRPG(pts + jj(1)+i) + factor*gbuf%STRPG(pts + jj(6)+i) * thk
801 value(2) = value(2)+ gbuf%STRPG(pts + jj(2)+i) + factor*gbuf%STRPG(pts + jj(7)+i) * thk
802 value(3) = value(3)+ gbuf%STRPG(pts + jj(3)+i) + factor*gbuf%STRPG(pts + jj(8)+i) * thk
803 ENDDO
804 value(1:3) = value(1:3)/npg
805 value(3) = value(3) * half
806 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
807 . shell_tensor,i,offset,nft,VALUE)
808 ENDDO
809 ELSE
810 DO i=1,nel
811 n = i + nft
812 thk = gbuf%THK(i)
813 factor = posly(i,ipt)
814 j = el2fa(nni+n)
815 value(1) = gbuf%STRA(jj(1)+i) + factor*gbuf%STRA(jj(6)+i) * thk
816 value(2) = gbuf%STRA(jj(2)+i) + factor*gbuf%STRA(jj(7)+i) * thk
817 value(3) = gbuf%STRA(jj(3)+i) + factor*gbuf%STRA(jj(8)+i) * thk
818 value(3) = value(3) * half
819 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
820 . shell_tensor,i,offset,nft,VALUE)
821 ENDDO
822 ENDIF ! NPG
823 ENDIF
824 ENDIF
825 DEALLOCATE(matly, thkly,posly,thk_ly)
826C--------------------------------------------------
827 ELSEIF (keyword == 'TENS/MSTRAIN') THEN ! strain tensor output in material coordinate system
828C--------------------------------------------------
829 ALLOCATE (epsm(nel,3))
830 epsm(:,:) = zero
831c
832 IF (iply > 0 .AND. ipt > 0) THEN
833
834 IF (igtyp == 17 .OR. igtyp == 19 .OR. igtyp == 51 .OR. igtyp == 52) THEN
835 ipang = 1
836 ipthk = ipang + nlay
837 ippos = ipthk + nlay
838 ipt_all = 0
839 DO j=1,nlay
840 bufly => elbuf_tab(ng)%BUFLY(j)
841 nptt = bufly%NPTT
842 IF (igtyp == 17 .OR. igtyp == 19 .OR. igtyp == 51) THEN
843 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
844 ELSEIF (igtyp == 52) THEN
845 id_ply = ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
846 ENDIF
847 IF (id_ply == iply .AND. ipt <= nptt) THEN
848 ilay = j
849 islice = ipt_all + ipt
850 IF(npg > 1) THEN
851 lens = nel*gbuf%G_STRPG/npg
852 DO i=1,nel
853 thk = gbuf%THK(i)
854 factor = posly(i,islice)
855 epsm(i,1:3) = zero
856 DO ipg = 1, npg
857 pts = (ipg -1)*lens
858 epsm(i,1) = epsm(i,1)+ gbuf%STRPG(pts + jj(1)+i) + factor*gbuf%STRPG(pts + jj(6)+i) * thk
859 epsm(i,2) = epsm(i,2)+ gbuf%STRPG(pts + jj(2)+i) + factor*gbuf%STRPG(pts + jj(7)+i) * thk
860 epsm(i,3) = epsm(i,3)+ gbuf%STRPG(pts + jj(3)+i) + factor*gbuf%STRPG(pts + jj(8)+i) * thk
861 ENDDO
862 epsm(i,1) = epsm(i,1)/npg
863 epsm(i,2) = epsm(i,2)/npg
864 epsm(i,3) = half*epsm(i,3)/npg
865 ENDDO
866 ELSE
867 DO i=1,nel
868 thk = gbuf%THK(i)
869 factor = posly(i,islice)
870 epsm(i,1) = gbuf%STRA(jj(1)+i) + factor*gbuf%STRA(jj(6)+i) * thk
871 epsm(i,2) = gbuf%STRA(jj(2)+i) + factor*gbuf%STRA(jj(7)+i) * thk
872 epsm(i,3) = gbuf%STRA(jj(3)+i) + factor*gbuf%STRA(jj(8)+i) * thk
873 epsm(i,3) = epsm(i,3) * half
874 ENDDO
875 ENDIF ! NPG
876 imat = elbuf_tab(ng)%BUFLY(ilay)%IMAT
877 mat_orth = mat_param(imat)%ORTHOTROPY
878 IF (mat_orth > 0) THEN
879 IF (idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52) ) THEN
880 dir_a => elbuf_tab(ng)%BUFLY(ilay)%LBUF_DIR(ipt)%DIRA
881 dir_b => elbuf_tab(ng)%BUFLY(ilay)%LBUF_DIR(ipt)%DIRB
882 ELSE
883 dir_a => elbuf_tab(ng)%BUFLY(ilay)%DIRA
884 dir_b => elbuf_tab(ng)%BUFLY(ilay)%DIRB
885 ENDIF
886 END IF
887 IF (mat_orth == 2) THEN
888 CALL roto_tens2d(nel,epsm,dir_a)
889 ELSE IF (mat_orth == 3) THEN ! anisotropic (law 58,158 only)
890 CALL roto_tens2d_aniso(nel,epsm,dir_a,dir_b)
891 END IF
892 ENDIF
893 ipt_all = ipt_all + nptt
894 ENDDO
895 ENDIF
896c
897 ELSEIF (ilay > 0 .AND. ilay <= nlay .AND. iply == -1 .AND. ipt == -1) THEN
898c PLY=NULL ILAYER=IL NPT=NULL
899 IF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16) THEN
900 IF(npg > 1) THEN
901 lens = nel*gbuf%G_STRPG/npg
902 DO i=1,nel
903 thk = gbuf%THK(i)
904 factor = posly(i,ilay)
905 epsm(i,1:3) = zero
906 DO ipg = 1, npg
907 pts = (ipg -1)*lens
908 epsm(i,1) = epsm(i,1)+ gbuf%STRPG(pts + jj(1)+i) + factor*gbuf%STRPG(pts + jj(6)+i) * thk
909 epsm(i,2) = epsm(i,2)+ gbuf%STRPG(pts + jj(2)+i) + factor*gbuf%STRPG(pts + jj(7)+i) * thk
910 epsm(i,3) = epsm(i,3)+ gbuf%STRPG(pts + jj(3)+i) + factor*gbuf%STRPG(pts + jj(8)+i) * thk
911 ENDDO
912 epsm(i,1) = epsm(i,1)/npg
913 epsm(i,2) = epsm(i,2)/npg
914 epsm(i,3) = half*epsm(i,3)/npg
915 ENDDO
916 ELSE
917 DO i=1,nel
918 thk = gbuf%THK(i)
919 factor = posly(i,ilay)
920 epsm(i,1) = gbuf%STRA(jj(1)+i) + factor*gbuf%STRA(jj(6)+i) * thk
921 epsm(i,2) = gbuf%STRA(jj(2)+i) + factor*gbuf%STRA(jj(7)+i) * thk
922 epsm(i,3) = gbuf%STRA(jj(3)+i) + factor*gbuf%STRA(jj(8)+i) * thk
923 epsm(i,3) = epsm(i,3) * half
924 ENDDO
925 ENDIF
926 imat = elbuf_tab(ng)%BUFLY(ilay)%IMAT
927 mat_orth = mat_param(imat)%ORTHOTROPY
928 IF (mat_orth > 0) THEN
929 dir_a => elbuf_tab(ng)%BUFLY(ilay)%DIRA
930 dir_b => elbuf_tab(ng)%BUFLY(ilay)%DIRB
931 END IF
932 IF (mat_orth == 2) THEN
933 CALL roto_tens2d(nel,epsm,dir_a)
934 ELSE IF (mat_orth == 3) THEN ! anisotropic (law 58,158 only)
935 CALL roto_tens2d_aniso(nel,epsm,dir_a,dir_b)
936 END IF
937 ENDIF
938c
939 ELSEIF (ipt > 0 .AND. ipt <= mpt .AND. iply == -1 .AND. ilay == -1) THEN
940c PLY=NULL ILAYER=NULL NPT=IPT
941 IF (igtyp == 1 .OR. igtyp == 9) THEN
942 IF(npg > 1) THEN
943 lens = nel*gbuf%G_STRPG/npg
944 DO i=1,nel
945 thk = gbuf%THK(i)
946 factor = posly(i,ipt)
947 epsm(i,1:3) = zero
948 DO ipg = 1, npg
949 pts = (ipg -1)*lens
950 epsm(i,1) = epsm(i,1)+ gbuf%STRPG(pts + jj(1)+i) + factor*gbuf%STRPG(pts + jj(6)+i) * thk
951 epsm(i,2) = epsm(i,2)+ gbuf%STRPG(pts + jj(2)+i) + factor*gbuf%STRPG(pts + jj(7)+i) * thk
952 epsm(i,3) = epsm(i,3)+ gbuf%STRPG(pts + jj(3)+i) + factor*gbuf%STRPG(pts + jj(8)+i) * thk
953 ENDDO
954 epsm(i,1) = epsm(i,1)/npg
955 epsm(i,2) = epsm(i,2)/npg
956 epsm(i,3) = half*epsm(i,3)/npg
957 ENDDO
958 ELSE
959 DO i=1,nel
960 thk = gbuf%THK(i)
961 factor = posly(i,ipt)
962 epsm(i,1) = gbuf%STRA(jj(1)+i) + factor*gbuf%STRA(jj(6)+i) * thk
963 epsm(i,2) = gbuf%STRA(jj(2)+i) + factor*gbuf%STRA(jj(7)+i) * thk
964 epsm(i,3) = gbuf%STRA(jj(3)+i) + factor*gbuf%STRA(jj(8)+i) * thk
965 epsm(i,3) = epsm(i,3) * half
966 ENDDO
967 ENDIF
968 imat = elbuf_tab(ng)%BUFLY(1)%IMAT
969 mat_orth = mat_param(imat)%ORTHOTROPY
970 IF (mat_orth == 2) THEN
971 dir_a => elbuf_tab(ng)%BUFLY(1)%DIRA
972 CALL roto_tens2d(nel,epsm,dir_a)
973 END IF
974 ENDIF
975 ENDIF
976c---
977 CALL h3d_write_sh_tensor_array(
978 . iok_part ,iselect ,nel ,offset ,nft ,
979 . is_written_shell,shell_tensor,epsm )
980
981 DEALLOCATE (epsm)
982 DEALLOCATE(matly, thkly,posly,thk_ly)
983C--------------------------------------------------
984 ELSEIF (keyword == 'TENS/EPSDOT/MEMB') THEN
985C--------------------------------------------------
986 a1 = one
987 a2 = zero
988 DO i=1,nel
989 thk = gbuf%THK(i)
990 value(1) = a1*epsdot(1,i+nft+offset) + a2*epsdot(4,i+nft+offset)*thk
991 value(2) = a1*epsdot(2,i+nft+offset) + a2*epsdot(5,i+nft+offset)*thk
992 value(3) = (a1*epsdot(3,i+nft+offset) + a2*epsdot(6,i+nft+offset)*thk)* half
993 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
994 . VALUE)
995 ENDDO
996C--------------------------------------------------
997 ELSEIF (keyword == 'TENS/EPSDOT/BEND') THEN
998C--------------------------------------------------
999 DO i=1,nel
1000 thk = gbuf%THK(i)
1001 value(1) = epsdot(4,i+nft+offset)
1002 value(2) = epsdot(5,i+nft+offset)
1003 value(3) = epsdot(6,i+nft+offset) * half
1004 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1005 . VALUE)
1006 ENDDO
1007C--------------------------------------------------
1008 ELSEIF (keyword == 'TENS/EPSDOT') THEN
1009C--------------------------------------------------
1010c ILAYER=NULL NPT=NULL
1011 IF ( ilay == -1 .AND. ipt == -1 .AND. iply == -1) THEN
1012 a1 = one
1013 a2 = zero
1014 DO i=1,nel
1015 thk = gbuf%THK(i)
1016 value(1) = a1*epsdot(1,i+nft+offset) + a2*epsdot(4,i+nft+offset)*thk
1017 value(2) = a1*epsdot(2,i+nft+offset) + a2*epsdot(5,i+nft+offset)*thk
1018 value(3) = (a1*epsdot(3,i+nft+offset) + a2*epsdot(6,i+nft+offset)*thk)* half
1019 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
1020 . shell_tensor,i,offset,nft,VALUE)
1021 ENDDO
1022c PLY=IPLY NPT=IPT
1023 ELSEIF ( iply > 0 .AND. ipt <= mpt .AND. ipt > 0 ) THEN
1024 IF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
1025 ipang = 1
1026 ipthk = ipang + nlay
1027 ippos = ipthk + nlay
1028 DO j=1,nlay
1029 IF (igtyp == 17 .OR. igtyp == 51) THEN
1030 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
1031 ELSEIF (igtyp == 52) THEN
1032 id_ply = ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
1033 ENDIF
1034 bufly => elbuf_tab(ng)%BUFLY(j)
1035 nptt = bufly%NPTT
1036 IF (id_ply == iply .AND. ipt <= nptt) THEN
1037 a2 = stack%GEO(ippos+j,isubstack)+
1038 . half*(((2*ipt-one)/nptt)-one) *
1039 . stack%GEO(ipthk+j,isubstack)
1040 DO i=1,nel
1041 thk = gbuf%THK(i)
1042 value(1) = epsdot(1,i+nft+offset) + a2*epsdot(4,i+nft+offset)*thk
1043 value(2) = epsdot(2,i+nft+offset) + a2*epsdot(5,i+nft+offset)*thk
1044 value(3) =(epsdot(3,i+nft+offset) + a2*epsdot(6,i+nft+offset)*thk)* half
1045 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,
1046 . shell_tensor,i,offset,nft,VALUE)
1047 ENDDO
1048 ENDIF
1049 ENDDO
1050 ENDIF
1051
1052c PLY=NULL ILAYER=ILAY NPT=IPT
1053 ELSEIF (iply == -1 .AND. ilay <= nlay .AND. ilay > 0 .AND. ipt <= mpt .AND. ipt > 0 ) THEN
1054 IF (igtyp == 51 .OR. igtyp == 52) THEN
1055 a1 = zero
1056 a2 = zero
1057 ns1 = 8
1058 ns2 = 8
1059 ipang = 1
1060 ipthk = ipang + nlay
1061 ippos = ipthk + nlay
1062 IF (igtyp == 17 .OR. igtyp == 51) THEN
1063 id_ply = igeo(1,stack%IGEO(2+ilay,isubstack))
1064 ELSEIF (igtyp == 52) THEN
1065 id_ply = ply_info(1,stack%IGEO(2+ilay,isubstack)-numstack)
1066 ENDIF
1067 bufly => elbuf_tab(ng)%BUFLY(ilay)
1068 nptt = bufly%NPTT
1069 IF (ipt <= nptt) THEN
1070 a1 = one
1071 a2 = stack%GEO(ippos+ilay,isubstack)+
1072 . half*(((2*ipt-one)/nptt)-one) *
1073 . stack%GEO(ipthk+ilay,isubstack)
1074 DO i=1,nel
1075 n = i + nft
1076 thk = gbuf%THK(i)
1077 value(1) = a1*epsdot(1,i+nft+offset) + a2*epsdot(4,i+nft+offset)*thk
1078 value(2) = a1*epsdot(2,i+nft+offset) + a2*epsdot(5,i+nft+offset)*thk
1079 value(3) = (a1*epsdot(3,i+nft+offset) + a2*epsdot(6,i+nft+offset)*thk)* half
1080 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1081 . VALUE)
1082 ENDDO
1083 ENDIF
1084 ENDIF
1085c PLY=NULL ILAYER=IL NPT=NULL
1086 ELSEIF (iply == -1 .AND. ilay <= nlay .AND. ilay > 0 .AND. ipt == -1 ) THEN
1087 IF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16 .OR. igtyp == 17) THEN
1088 a1 = one
1089 a2 = half*(((2*ilay-one)/nlay)-one)
1090 DO i=1,nel
1091 n = i + nft
1092 thk = gbuf%THK(i)
1093 value(1) = a1*epsdot(1,i+nft+offset) + a2*epsdot(4,i+nft+offset)*thk
1094 value(2) = a1*epsdot(2,i+nft+offset) + a2*epsdot(5,i+nft+offset)*thk
1095 value(3) = (a1*epsdot(3,i+nft+offset) + a2*epsdot(6,i+nft+offset)*thk)* half
1096 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1097 . VALUE)
1098 ENDDO
1099 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
1100 a1 = one
1101 a2 = stack%GEO(ippos+ilay,isubstack)+
1102 . half*(((2*ipt-one)/nptt)-one) *
1103 . stack%GEO(ipthk+ilay,isubstack)
1104 DO i=1,nel
1105 n = i + nft
1106 thk = gbuf%THK(i)
1107 value(1) = a1*epsdot(1,i+nft+offset) + a2*epsdot(4,i+nft+offset)*thk
1108 value(2) = a1*epsdot(2,i+nft+offset) + a2*epsdot(5,i+nft+offset)*thk
1109 value(3) = (a1*epsdot(3,i+nft+offset) + a2*epsdot(6,i+nft+offset)*thk)* half
1110 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1111 . VALUE)
1112 ENDDO
1113 ENDIF
1114c PLY=NULL ILAYER=NULL NPT=IPT
1115 ELSEIF ( ipt <= mpt .AND. ipt > 0) THEN
1116 a1 = one
1117 a2 = half*(((2*ipt-one)/mpt)-one)
1118 IF (igtyp == 1 .OR. igtyp == 9) THEN
1119 DO i=1,nel
1120 thk = gbuf%THK(i)
1121 value(1) = a1*epsdot(1,i+nft+offset) + a2*epsdot(4,i+nft+offset)*thk
1122 value(2) = a1*epsdot(2,i+nft+offset) + a2*epsdot(5,i+nft+offset)*thk
1123 value(3) = (a1*epsdot(3,i+nft+offset) + a2*epsdot(6,i+nft+offset)*thk)* half
1124 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1125 . VALUE)
1126 ENDDO
1127 ENDIF
1128 ENDIF
1129C--------------------------------------------------
1130 ELSE IF (keyword == 'TENS/STRAIN_ENG') THEN
1131C--------------------------------------------------
1132 IF (ity == 3 ) THEN !4n
1133 DO i=1,nel
1134 n = i + nft
1135 nni = ixc(2,n)
1136 j = 4*(i-1) +1
1137 xn(j)=x(1,nni)
1138 yn(j)=x(2,nni)
1139 zn(j)=x(3,nni)
1140 dxn(j)=d(1,nni)
1141 dyn(j)=d(2,nni)
1142 dzn(j)=d(3,nni)
1143 nni = ixc(3,n)
1144 xn(j+1)=x(1,nni)
1145 yn(j+1)=x(2,nni)
1146 zn(j+1)=x(3,nni)
1147 dxn(j+1)=d(1,nni)
1148 dyn(j+1)=d(2,nni)
1149 dzn(j+1)=d(3,nni)
1150 nni = ixc(4,n)
1151 xn(j+2)=x(1,nni)
1152 yn(j+2)=x(2,nni)
1153 zn(j+2)=x(3,nni)
1154 dxn(j+2)=d(1,nni)
1155 dyn(j+2)=d(2,nni)
1156 dzn(j+2)=d(3,nni)
1157 nni = ixc(5,n)
1158 xn(j+3)=x(1,nni)
1159 yn(j+3)=x(2,nni)
1160 zn(j+3)=x(3,nni)
1161 dxn(j+3)=d(1,nni)
1162 dyn(j+3)=d(2,nni)
1163 dzn(j+3)=d(3,nni)
1164 strain(1:3,i)=zero
1165 ENDDO
1166 CALL sh4_tstrain(xn,yn,zn,dxn,dyn,dzn,strain,nel)
1167 DO i=1,nel
1168 value(1:3)= strain(1:3,i)
1169 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1170 . VALUE)
1171 ENDDO
1172 ELSEIF (ity == 7) THEN
1173 DO i=1,nel
1174 n = i + nft
1175 nni = ixtg(2,n)
1176 j = 3*(i-1) +1
1177 xn(j)=x(1,nni)
1178 yn(j)=x(2,nni)
1179 zn(j)=x(3,nni)
1180 dxn(j)=d(1,nni)
1181 dyn(j)=d(2,nni)
1182 dzn(j)=d(3,nni)
1183 nni = ixtg(3,n)
1184 xn(j+1)=x(1,nni)
1185 yn(j+1)=x(2,nni)
1186 zn(j+1)=x(3,nni)
1187 dxn(j+1)=d(1,nni)
1188 dyn(j+1)=d(2,nni)
1189 dzn(j+1)=d(3,nni)
1190 nni = ixtg(4,n)
1191 xn(j+2)=x(1,nni)
1192 yn(j+2)=x(2,nni)
1193 zn(j+2)=x(3,nni)
1194 dxn(j+2)=d(1,nni)
1195 dyn(j+2)=d(2,nni)
1196 dzn(j+2)=d(3,nni)
1197 strain(1:3,i)=zero
1198 ENDDO
1199 CALL sh3_tstrain(xn,yn,zn,dxn,dyn,dzn,strain,nel,ihbe)
1200 DO i=1,nel
1201 value(1:3)= strain(1:3,i)
1202 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1203 . VALUE)
1204 ENDDO
1205 END IF
1206C--------------------------------------------------
1207C-----------------------------------------------
1208 ELSEIF (keyword == 'TENS/STRESS/TMAX') THEN
1209C---------- -------------------------------------
1210 DO i=1,nel
1211 value(1:3) = gbuf%TM_SIG1(jj(1:3) + i)
1212 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1213 . VALUE)
1214 ENDDO
1215C---------- -------------------------------------
1216 ELSEIF (keyword == 'TENS/STRESS/TMIN') THEN
1217C---------- -------------------------------------
1218 DO i=1,nel
1219 value(1:3) = gbuf%TM_SIG3(jj(1:3) + i)
1220 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1221 . VALUE)
1222 ENDDO
1223C---------- -------------------------------------
1224 ELSEIF (keyword == 'TENS/STRAIN/TMAX') THEN
1225C---------- -------------------------------------
1226 DO i=1,nel
1227 value(1:3) = gbuf%TM_STRA1(jj(1:3) + i)
1228 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1229 . VALUE)
1230 ENDDO
1231C---------- -------------------------------------
1232 ELSEIF (keyword == 'TENS/STRAIN/TMIN') THEN
1233C---------- -------------------------------------
1234 DO i=1,nel
1235 value(1:3) = gbuf%TM_STRA3(jj(1:3) + i)
1236 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,
1237 . VALUE)
1238 ENDDO
1239C--------------------------------------------------
1240 ELSEIF (keyword == 'TENS/BSTRESS') THEN
1241C--------------------------------------------------
1242 IF (mlw == 87) THEN
1243 imat = ixc(1,nft+1)
1244 iadbuf = ipm(7,imat)
1245 nuparam= ipm(9,imat)
1246 uparam => bufmat(iadbuf:iadbuf+nuparam-1)
1247 nbfunct = uparam(25)
1248 nchard = 34 + 2*nbfunct + 22
1249 chard = uparam(nchard)
1250 ELSEIF (mlw == 36) THEN
1251 imat = ixc(1,nft+1)
1252 iadbuf = ipm(7,imat)
1253 nuparam= ipm(9,imat)
1254 uparam => bufmat(iadbuf:iadbuf+nuparam-1)
1255 nbfunct = uparam(1)
1256 nchard = 2*nbfunct + 14
1257 chard = uparam(nchard)
1258 ENDIF
1259 IF ( ilay == -1 .AND. ipt == -1 .AND. iply == -1) THEN !global value = mean on gauss points IPs and layers
1260 IF(id == -1) THEN ! sum of all backstresses
1261 IF(mlw == 36 .AND. chard > zero) THEN
1262 DO i=1,nel
1263 DO j=1,nlay
1264 bufly => elbuf_tab(ng)%BUFLY(j)
1265 nptt = bufly%NPTT
1266 DO ir=1,nptr
1267 DO is=1,npts
1268 DO it=1,nptt
1269 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1270 DO k=1,3
1271 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg/nptt/nlay
1272 ENDDO !K
1273 ENDDO !IT
1274 ENDDO !IS
1275 ENDDO !IR
1276 ENDDO !Jlay
1277 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1278 ENDDO !I
1279 ELSEIF(mlw == 78) THEN
1280 DO i=1,nel
1281 DO j=1,nlay
1282 bufly => elbuf_tab(ng)%BUFLY(j)
1283 nptt = bufly%NPTT
1284 DO ir=1,nptr
1285 DO is=1,npts
1286 DO it=1,nptt
1287 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1288 DO k=1,3
1289 value(k) = value(k) + (lbuf%SIGA(jj(k) + i)+lbuf%SIGB(jj(k) + i))/npg/nptt/nlay
1290 ENDDO !K
1291 ENDDO !IT
1292 ENDDO !IS
1293 ENDDO !IR
1294 ENDDO !Jlay
1295 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1296 ENDDO !I
1297 ELSEIF(mlw == 87 .AND. chard > zero) THEN
1298 !SIGBXX(I) = SIGB(I,1) + SIGB(I,4) + SIGB(I,7) + SIGB(I,10)
1299 !SIGBYY(I) = SIGB(I,2) + SIGB(I,5) + SIGB(I,8) + SIGB(I,11)
1300 !SIGBXY(I) = SIGB(I,3) + SIGB(I,6) + SIGB(I,9) + SIGB(I,12)
1301 DO i=1,nel
1302 DO j=1,nlay
1303 bufly => elbuf_tab(ng)%BUFLY(j)
1304 nptt = bufly%NPTT
1305 DO ir=1,nptr
1306 DO is=1,npts
1307 DO it=1,nptt
1308 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1309 DO k=1,3
1310 value(k) = value(k) + (lbuf%SIGB(jj(k) + i )
1311 . +lbuf%SIGB(jj(k+3) + i )
1312 . +lbuf%SIGB(jj(k+6) + i )
1313 . +lbuf%SIGB(jj(k+9) + i ))/npg/nptt/nlay
1314 ENDDO !K
1315 ENDDO !it
1316 ENDDO !IS
1317 ENDDO !IR
1318 ENDDO !Jlay
1319 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1320 ENDDO !I
1321 ENDIF ! MLW==
1322 ELSEIF(id > 0) THEN !!
1323 IF(mlw == 36.AND. chard > zero) THEN ! necessarily id = 1 there is only one bs
1324 DO i=1,nel
1325 DO j=1,nlay
1326 bufly => elbuf_tab(ng)%BUFLY(j)
1327 nptt = bufly%NPTT
1328 DO ir=1,nptr
1329 DO is=1,npts
1330 DO it=1,nptt
1331 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1332 DO k=1,3
1333 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg/nptt/nlay
1334 ENDDO !K
1335 ENDDO !IT
1336 ENDDO !IS
1337 ENDDO !IR
1338 ENDDO !Jlay
1339 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1340 ENDDO !I
1341 ELSEIF(mlw == 78) THEN
1342 IF(id == 1) THEN
1343 DO i=1,nel
1344 DO j=1,nlay
1345 bufly => elbuf_tab(ng)%BUFLY(j)
1346 nptt = bufly%NPTT
1347 DO ir=1,nptr
1348 DO is=1,npts
1349 DO it=1,nptt
1350 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1351 DO k=1,3
1352 value(k) = VALUE(k) + lbuf%SIGA(jj(k) + i) /npg/nptt/nlay
1353 ENDDO !K
1354 ENDDO !IT
1355 ENDDO !IS
1356 ENDDO !IR
1357 ENDDO !Jlay
1358 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1359 ENDDO !I
1360 ELSEIF(id == 2) THEN
1361 DO i=1,nel
1362 DO j=1,nlay
1363 bufly => elbuf_tab(ng)%BUFLY(j)
1364 nptt = bufly%NPTT
1365 DO ir=1,nptr
1366 DO is=1,npts
1367 DO it=1,nptt
1368 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1369 DO k=1,3
1370 value(k) = value(k) + lbuf%SIGB(jj(k) + i) /npg/nptt/nlay
1371 ENDDO !K
1372 ENDDO !IT
1373 ENDDO !IS
1374 ENDDO !IR
1375 ENDDO !Jlay
1376 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1377 ENDDO !I
1378 ELSEIF(id == 3) THEN
1379 DO i=1,nel
1380 DO j=1,nlay
1381 bufly => elbuf_tab(ng)%BUFLY(j)
1382 nptt = bufly%NPTT
1383 DO ir=1,nptr
1384 DO is=1,npts
1385 DO it=1,nptt
1386 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1387 DO k=1,3
1388 value(k) = value(k) + lbuf%SIGC(jj(k) + i) /npg/nptt/nlay
1389 ENDDO !K
1390 ENDDO !IT
1391 ENDDO !IS
1392 ENDDO !ir
1393 ENDDO !Jlay
1394 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1395 ENDDO !I
1396 ENDIF !ID ==
1397 ELSEIF(mlw == 87.AND. chard > zero) THEN
1398 IF(id == 1) THEN
1399 DO i=1,nel
1400 DO j=1,nlay
1401 bufly => elbuf_tab(ng)%BUFLY(j)
1402 nptt = bufly%NPTT
1403 DO ir=1,nptr
1404 DO is=1,npts
1405 DO it=1,nptt
1406 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1407 DO k=1,3
1408 value(k) = value(k) + lbuf%SIGB(jj(k) + i ) /npg/nptt/nlay
1409 ENDDO !K
1410 ENDDO !IT
1411 ENDDO !IS
1412 ENDDO !IR
1413 ENDDO !Jlay
1414 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1415 ENDDO !I
1416 ELSEIF(id == 2) THEN
1417 DO i=1,nel
1418 DO j=1,nlay
1419 bufly => elbuf_tab(ng)%BUFLY(j)
1420 nptt = bufly%NPTT
1421 DO ir=1,nptr
1422 DO is=1,npts
1423 DO it=1,nptt
1424 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1425 DO k=1,3
1426 value(k) = value(k) + lbuf%SIGB(jj(k+3) + i) /npg/nptt/nlay
1427 ENDDO !K
1428 ENDDO !IT
1429 ENDDO !IS
1430 ENDDO !IR
1431 ENDDO !Jlay
1432 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1433 ENDDO !I
1434 ELSEIF(id == 3) THEN
1435 DO i=1,nel
1436 DO j=1,nlay
1437 bufly => elbuf_tab(ng)%BUFLY(j)
1438 nptt = bufly%NPTT
1439 DO ir=1,nptr
1440 DO is=1,npts
1441 DO it=1,nptt
1442 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1443 DO k=1,3
1444 value(k) = value(k) + lbuf%SIGB(jj(k+6) + i) /npg/nptt/nlay
1445 ENDDO !K
1446 ENDDO !IT
1447 ENDDO !IS
1448 ENDDO !IR
1449 ENDDO !Jlay
1450 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1451 ENDDO !I
1452 ELSEIF( id == 4) THEN
1453 DO i=1,nel
1454 DO j=1,nlay
1455 bufly => elbuf_tab(ng)%BUFLY(j)
1456 nptt = bufly%NPTT
1457 DO ir=1,nptr
1458 DO is=1,npts
1459 DO it=1,nptt
1460 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1461 DO k=1,3
1462 value(k) = value(k) + lbuf%SIGB(jj(k+9) + i)/npg/nptt/nlay
1463 ENDDO !K
1464 ENDDO !IT
1465 ENDDO !IS
1466 ENDDO !IR
1467 ENDDO !Jlay
1468 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE )
1469 ENDDO !I
1470 ENDIF !(ID == )
1471 endif!(MLW ==
1472 ENDIF !(ID >0 )
1473C--------------------------------------------------
1474c Ply = iPly npt = IPT for prop with Ply 17-51-52
1475C--------------------------------------------------
1476 ELSEIF ( iply > 0 .AND. ipt <= mpt .AND. ipt > 0 ) THEN
1477 DO j=1,nlay
1478 id_ply = 0
1479 IF (igtyp == 17 .OR. igtyp == 51) THEN
1480 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
1481 ELSEIF (igtyp == 52) THEN
1482 id_ply = ply_info(1,stack%IGEO(2+j,isubstack) - numstack)
1483 ENDIF
1484 IF (id_ply == iply) THEN
1485 bufly => elbuf_tab(ng)%BUFLY(j)
1486 !--------
1487 ! LAW36
1488 !--------
1489 IF (mlw == 36 .AND.( id == -1 .OR. id == 1).AND. chard > zero) THEN
1490 DO i=1,nel
1491 DO ir=1,nptr
1492 DO is=1,npts
1493 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1494 DO k=1,3
1495 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg
1496 ENDDO !k
1497 ENDDO !IS
1498 ENDDO !IR
1499 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1500 ENDDO ! I
1501 !--------
1502 ! LAW78
1503 !--------
1504 ELSEIF (mlw == 78) THEN
1505 IF(id == -1) THEN ! somme of all backstresses
1506 DO i=1,nel
1507 DO ir=1,nptr
1508 DO is=1,npts
1509 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1510 DO k=1,3
1511 value(k) = value(k) + (lbuf%SIGA(jj(k) + i)+lbuf%SIGB(jj(k) + i))/npg
1512 ENDDO !k
1513 ENDDO !IS
1514 ENDDO !IR
1515 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1516 ENDDO ! I
1517 ELSEIF(id ==1 ) THEN
1518 DO i=1,nel
1519 DO ir=1,nptr
1520 DO is=1,npts
1521 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1522 DO k=1,3
1523 value(k) = value(k) + lbuf%SIGA(jj(k) + i)/npg
1524 ENDDO !k
1525 ENDDO !IS
1526 ENDDO !IR
1527 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1528 ENDDO ! I
1529 ELSEIF(id ==2 ) THEN
1530 DO i=1,nel
1531 DO ir=1,nptr
1532 DO is=1,npts
1533 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1534 DO k=1,3
1535 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg
1536 ENDDO !k
1537 ENDDO !IS
1538 ENDDO !IR
1539 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1540 ENDDO ! I
1541
1542 ELSEIF(id ==3 ) THEN
1543 DO i=1,nel
1544 DO ir=1,nptr
1545 DO is=1,npts
1546 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1547 DO k=1,3
1548 VALUE(k) = value(k) + lbuf%SIGC(jj(k) + i)/npg
1549 ENDDO !k
1550 ENDDO !IS
1551 ENDDO !IR
1552 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1553 ENDDO ! I
1554 ENDIF !ID == -1
1555 !--------
1556 ! LAW87
1557 !--------
1558 ELSEIF( mlw == 87 .AND. chard > zero) THEN
1559 IF(id == -1) THEN ! somme of all backstresses
1560 DO i=1,nel
1561 DO ir=1,nptr
1562 DO is=1,npts
1563 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1564 DO k=1,3
1565 value(k) = value(k) + (lbuf%SIGB(jj(k) + i )
1566 . +lbuf%SIGB(jj(k+3) + i )
1567 . +lbuf%SIGB(jj(k+6) + i )
1568 . +lbuf%SIGB(jj(k+9) + i ))/npg
1569
1570 ENDDO !k
1571 ENDDO !IS
1572 ENDDO !IR
1573 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1574 ENDDO ! I
1575 ELSEIF(id ==1 ) THEN
1576 DO i=1,nel
1577 DO ir=1,nptr
1578 DO is=1,npts
1579 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1580 DO k=1,3
1581 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg
1582 ENDDO !k
1583 ENDDO !IS
1584 ENDDO !IR
1585 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1586 ENDDO ! I
1587 ELSEIF(id ==2 ) THEN
1588 DO i=1,nel
1589 DO ir=1,nptr
1590 DO is=1,npts
1591 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1592 DO k=1,3
1593 value(k) = value(k) + lbuf%SIGB(jj(k+3) + i)/npg
1594 ENDDO !k
1595 ENDDO !IS
1596 ENDDO !IR
1597 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1598 ENDDO ! I
1599 ELSEIF(id ==3 ) THEN
1600 DO i=1,nel
1601 DO ir=1,nptr
1602 DO is=1,npts
1603 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1604 DO k=1,3
1605 value(k) = value(k) + lbuf%SIGB(jj(k+6) + i)/npg
1606 ENDDO !k
1607 ENDDO !IS
1608 ENDDO !IR
1609 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1610 ENDDO ! I
1611 ELSEIF(id ==4 ) THEN
1612 DO i=1,nel
1613 DO ir=1,nptr
1614 DO is=1,npts
1615 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1616 DO k=1,3
1617 value(k) = value(k) + lbuf%SIGB(jj(k+9) + i)/npg
1618 ENDDO !k
1619 ENDDO !IS
1620 ENDDO !IR
1621 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1622 ENDDO ! I
1623 ENDIF !ID == -1
1624 ENDIF !(MLW ==
1625 END IF !(ID_PLY == IPLY
1626 ENDDO !JLAY
1627C--------------------------------------------------
1628c PLY=NULL ILAYER=IL NPT=IPT
1629C--------------------------------------------------
1630 ELSEIF (ilay > 0 .AND. ilay <= nlay .AND. ipt <= mpt .AND. ipt > 0 ) THEN
1631 j = ilay
1632 IF(igtyp == 9) j = 1
1633 bufly => elbuf_tab(ng)%BUFLY(j)
1634 !--------
1635 ! LAW36
1636 !--------
1637 IF (mlw == 36.AND. (id==-1 . or .id==1).AND. chard > zero) THEN
1638 DO i=1,nel
1639 DO ir=1,nptr
1640 DO is=1,npts
1641 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1642 DO k=1,3
1643 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg
1644 ENDDO !k
1645 ENDDO !IS
1646 ENDDO !IR
1647 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1648 ENDDO ! I
1649 !--------
1650 ! LAW78
1651 !--------
1652 ELSEIF (mlw == 78) THEN
1653 IF(id == -1) THEN ! somme of all backstresses
1654 DO i=1,nel
1655 DO ir=1,nptr
1656 DO is=1,npts
1657 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1658 DO k=1,3
1659 value(k) = value(k) + (lbuf%SIGA(jj(k) + i)+lbuf%SIGB(jj(k) + i))/npg
1660 ENDDO !k
1661 ENDDO !IS
1662 ENDDO !IR
1663 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1664 ENDDO ! I
1665 ELSEIF(id ==1 ) THEN
1666 DO i=1,nel
1667 DO ir=1,nptr
1668 DO is=1,npts
1669 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1670 DO k=1,3
1671 value(k) = value(k) + lbuf%SIGA(jj(k) + i)/npg
1672 ENDDO !k
1673 ENDDO !IS
1674 ENDDO !IR
1675 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1676 ENDDO ! I
1677 ELSEIF(id ==2 ) THEN
1678 DO i=1,nel
1679 DO ir=1,nptr
1680 DO is=1,npts
1681 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1682 DO k=1,3
1683 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg
1684 ENDDO !k
1685 ENDDO !IS
1686 ENDDO !IR
1687 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1688 ENDDO ! I
1689 ELSEIF(id ==3 ) THEN
1690 DO i=1,nel
1691 DO ir=1,nptr
1692 DO is=1,npts
1693 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1694 DO k=1,3
1695 value(k) = value(k) + lbuf%SIGC(jj(k) + i)/npg
1696 ENDDO !k
1697 ENDDO !IS
1698 ENDDO !IR
1699 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1700 ENDDO ! I
1701 ENDIF !ID == -1
1702 !--------
1703 ! LAW87
1704 !--------
1705 ELSEIF( mlw == 87 .AND. chard > zero) THEN
1706 IF(id == -1) THEN ! somme of all backstresses
1707 DO i=1,nel
1708 DO ir=1,nptr
1709 DO is=1,npts
1710 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1711 DO k=1,3
1712 value(k) = value(k) + (lbuf%SIGB(jj(k) + i )
1713 . +lbuf%SIGB(jj(k+3) + i )
1714 . +lbuf%SIGB(jj(k+6) + i )
1715 . +lbuf%SIGB(jj(k+9) + i ))/npg
1716
1717 ENDDO !k
1718 ENDDO !IS
1719 ENDDO !IR
1720 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1721 ENDDO ! I
1722 ELSEIF(id ==1 ) THEN
1723 DO i=1,nel
1724 DO ir=1,nptr
1725 DO is=1,npts
1726 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1727 DO k=1,3
1728 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg
1729 ENDDO !k
1730 ENDDO !IS
1731 ENDDO !IR
1732 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1733 ENDDO ! I
1734 ELSEIF(id ==2 ) THEN
1735 DO i=1,nel
1736 DO ir=1,nptr
1737 DO is=1,npts
1738 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1739 DO k=1,3
1740 value(k) = value(k) + lbuf%SIGB(jj(k+3) + i)/npg
1741 ENDDO !k
1742 ENDDO !IS
1743 ENDDO !IR
1744 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1745 ENDDO ! I
1746 ELSEIF(id ==3 ) THEN
1747 DO i=1,nel
1748 DO ir=1,nptr
1749 DO is=1,npts
1750 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1751 DO k=1,3
1752 value(k) = value(k) + lbuf%SIGB(jj(k+6) + i)/npg
1753 ENDDO !k
1754 ENDDO !IS
1755 ENDDO !IR
1756 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1757 ENDDO ! I
1758 ELSEIF(id ==4 ) THEN
1759 DO i=1,nel
1760 DO ir=1,nptr
1761 DO is=1,npts
1762 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1763 DO k=1,3
1764 value(k) = value(k) + lbuf%SIGB(jj(k+9) + i)/npg
1765 ENDDO !k
1766 ENDDO !IS
1767 ENDDO !IR
1768 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1769 ENDDO ! I
1770 ENDIF !ID == -1
1771 ENDIF !(MLW ==
1772C--------------------------------------------------
1773c PLY=NULL ILAYER=IL NPT=NULL ! prop 9-10-11 have layers and not PLY
1774C--------------------------------------------------
1775 ELSEIF (iply == -1 .AND. ilay <= nlay .AND. ilay > 0 .AND. ipt == -1 ) THEN
1776 IF (igtyp == 9 .OR.igtyp == 10 .OR. igtyp == 11 ) THEN
1777 ! output in orthotopic frame / elementary
1778 j = ilay
1779 IF(igtyp == 9) j = 1
1780 bufly => elbuf_tab(ng)%BUFLY(j)
1781 nptt = bufly%NPTT
1782 !--------
1783 ! LAW36
1784 !--------
1785 IF (mlw == 36.AND. (id==-1 .OR. id==1) .AND. chard > zero) THEN ! only one Bstress
1786 DO i=1,nel
1787 DO ir=1,nptr
1788 DO is=1,npts
1789 DO it=1,nptt
1790 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1791 DO k=1,3
1792 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg/nptt
1793 ENDDO !K
1794 enddo!IT
1795 ENDDO !IS
1796 enddo!IR
1797 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1798 value(1:3) = zero
1799 ENDDO !I
1800 !--------
1801 ! LAW78
1802 !--------
1803 ELSEIF (mlw == 78) THEN
1804 IF(id == -1) THEN ! somme of all backstresses
1805 DO i=1,nel
1806 DO ir=1,nptr
1807 DO is=1,npts
1808 DO it=1,nptt
1809 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1810 DO k=1,3
1811 value(k) = value(k) + (lbuf%SIGA(jj(k) + i)+lbuf%SIGB(jj(k) + i))/npg/nptt
1812 ENDDO !k
1813 enddo!IT
1814 ENDDO !IS
1815 ENDDO !IR
1816 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1817 ENDDO ! I
1818 ELSEIF(id ==1 ) THEN
1819 DO i=1,nel
1820 DO ir=1,nptr
1821 DO is=1,npts
1822 DO it=1,nptt
1823 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1824 DO k=1,3
1825 value(k) = value(k) + lbuf%SIGA(jj(k) + i)/npg/nptt
1826 ENDDO !k
1827 enddo!IT
1828 ENDDO !IS
1829 ENDDO !IR
1830 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1831 ENDDO ! I
1832 ELSEIF(id ==2 ) THEN
1833 DO i=1,nel
1834 DO ir=1,nptr
1835 DO is=1,npts
1836 DO it=1,nptt
1837 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1838 DO k=1,3
1839 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg/nptt
1840 ENDDO !k
1841 enddo!IT
1842 ENDDO !IS
1843 ENDDO !IR
1844 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1845 ENDDO ! I
1846 ELSEIF(id ==3 ) THEN
1847 DO i=1,nel
1848 DO ir=1,nptr
1849 DO is=1,npts
1850 DO it=1,nptt
1851 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1852 DO k=1,3
1853 value(k) = value(k) + lbuf%SIGC(jj(k) + i)/npg/nptt
1854 ENDDO !k
1855 enddo!IT
1856 ENDDO !IS
1857 ENDDO !IR
1858 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1859 ENDDO ! I
1860 ENDIF !ID == -1
1861 !--------
1862 ! LAW87
1863 !--------
1864 ELSEIF( mlw == 87 .AND. chard > zero) THEN
1865 IF(id == -1) THEN ! somme of all backstresses
1866 DO i=1,nel
1867 DO ir=1,nptr
1868 DO is=1,npts
1869 DO it=1,nptt
1870 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1871 DO k=1,3
1872 value(k) = value(k) + (lbuf%SIGB(jj(k) + i )
1873 . +lbuf%SIGB(jj(k+3) + i )
1874 . +lbuf%SIGB(jj(k+6) + i )
1875 . +lbuf%SIGB(jj(k+9) + i ))/npg/nptt
1876
1877 ENDDO !k
1878 enddo!IT
1879 ENDDO !IS
1880 ENDDO !IR
1881 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1882 ENDDO ! I
1883 ELSEIF(id ==1 ) THEN
1884 DO i=1,nel
1885 DO ir=1,nptr
1886 DO is=1,npts
1887 DO it=1,nptt
1888 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1889 DO k=1,3
1890 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg/nptt
1891 ENDDO !k
1892 enddo!IT
1893 ENDDO !IS
1894 ENDDO !IR
1895 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1896 ENDDO ! I
1897 ELSEIF(id ==2 ) THEN
1898 DO i=1,nel
1899 DO ir=1,nptr
1900 DO is=1,npts
1901 DO it=1,nptt
1902 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1903 DO k=1,3
1904 value(k) = value(k) + lbuf%SIGB(jj(k+3) + i )/npg/nptt
1905 ENDDO !k
1906 enddo!IT
1907 ENDDO !IS
1908 ENDDO !IR
1909 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1910 ENDDO ! I
1911 ELSEIF(id ==3 ) THEN
1912 DO i=1,nel
1913 DO ir=1,nptr
1914 DO is=1,npts
1915 DO it=1,nptt
1916 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1917 DO k=1,3
1918 value(k) = value(k) + lbuf%SIGB(jj(k+6) + i )/npg/nptt
1919 ENDDO !k
1920 enddo!IT
1921 ENDDO !IS
1922 ENDDO !IR
1923 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1924 ENDDO ! I
1925 ELSEIF(id ==4 ) THEN
1926 DO i=1,nel
1927 DO ir=1,nptr
1928 DO is=1,npts
1929 DO it=1,nptt
1930 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
1931 DO k=1,3
1932 value(k) = value(k) + lbuf%SIGB(jj(k+9) + i )/npg
1933 ENDDO !k
1934 enddo!IT
1935 ENDDO !IS
1936 ENDDO !IR
1937 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1938 ENDDO ! I
1939 ENDIF !ID == -1
1940 ENDIF !(MLW ==
1941 ENDIF !IGTYP ==
1942C---------------------------------------------------------------------------
1943c ILAYER=-1 NPT=IPT IPLY=-1 ID=-1
1944C---------------------------------------------------------------------------
1945 ELSE IF(ilay == -1 .AND. ipt > 0 .AND. ipt<=mpt .AND. iply == -1 ) THEN ! output for each layer/PLY as if LAYER==ALL
1946 DO j=1,nlay
1947 bufly => elbuf_tab(ng)%BUFLY(j)
1948 nptt = bufly%NPTT
1949 IF (ipt <= nptt ) THEN
1950 !--------
1951 ! LAW36
1952 !--------
1953 IF (mlw == 36.AND. (id==-1 .OR. id==1) .AND. chard > zero) THEN
1954 DO i=1,nel
1955 DO ir=1,nptr
1956 DO is=1,npts
1957 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1958 DO k=1,3
1959 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg
1960 ENDDO !k
1961 ENDDO !IS
1962 ENDDO !IR
1963 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1964 ENDDO ! I
1965 !--------
1966 ! LAW78
1967 !--------
1968 ELSEIF (mlw == 78) THEN
1969 IF(id == -1) THEN ! somme of all backstresses
1970 DO i=1,nel
1971 DO ir=1,nptr
1972 DO is=1,npts
1973 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1974 DO k=1,3
1975 value(k) = value(k) + (lbuf%SIGA(jj(k) + i)+lbuf%SIGB(jj(k) + i))/npg
1976 ENDDO !k
1977 ENDDO !IS
1978 ENDDO !IR
1979 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1980 ENDDO ! I
1981 ELSEIF(id ==1 ) THEN
1982 DO i=1,nel
1983 DO ir=1,nptr
1984 DO is=1,npts
1985 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1986 DO k=1,3
1987 value(k) = value(k) + lbuf%SIGA(jj(k) + i)/npg
1988 ENDDO !k
1989 ENDDO !IS
1990 ENDDO !IR
1991 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
1992 ENDDO ! I
1993 ELSEIF(id ==2 ) THEN
1994 DO i=1,nel
1995 DO ir=1,nptr
1996 DO is=1,npts
1997 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
1998 DO k=1,3
1999 value(k) = value(k) + lbuf%SIGB(jj(k) + i)/npg
2000 ENDDO !k
2001 ENDDO !IS
2002 ENDDO !IR
2003 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
2004 ENDDO ! I
2005 ELSEIF(id ==3 ) THEN
2006 DO i=1,nel
2007 DO ir=1,nptr
2008 DO is=1,npts
2009 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
2010 DO k=1,3
2011 value(k) = value(k) + lbuf%SIGC(jj(k) + i)/npg
2012 ENDDO !k
2013 ENDDO !IS
2014 ENDDO !IR
2015 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
2016 ENDDO ! I
2017 ENDIF !ID == -1
2018 !--------
2019 ! LAW87
2020 !--------
2021 ELSEIF( mlw == 87 .AND. chard > zero) THEN
2022 IF(id == -1) THEN ! somme of all backstresses
2023 DO i=1,nel
2024 DO ir=1,nptr
2025 DO is=1,npts
2026 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
2027 DO k=1,3
2028 value(k) = value(k) + (lbuf%SIGB(jj(k) + i )
2029 . +lbuf%SIGB(jj(k+3) + i )
2030 . +lbuf%SIGB(jj(k+6) + i )
2031 . +lbuf%SIGB(jj(k+9) + i ))/npg
2032
2033 ENDDO !k
2034 ENDDO !IS
2035 ENDDO !IR
2036 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
2037 ENDDO ! I
2038 ELSEIF(id ==1 ) THEN
2039 DO i=1,nel
2040 DO ir=1,nptr
2041 DO is=1,npts
2042 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
2043 DO k=1,3
2044 value(k) = value(k) + lbuf%SIGB(jj(k) + i )/npg
2045 ENDDO !k
2046 ENDDO !IS
2047 ENDDO !IR
2048 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
2049 ENDDO ! I
2050 ELSEIF(id ==2 ) THEN
2051 DO i=1,nel
2052 DO ir=1,nptr
2053 DO is=1,npts
2054 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
2055 DO k=1,3
2056 value(k) = value(k) + lbuf%SIGB(jj(k+3) + i )/npg
2057 ENDDO !k
2058 ENDDO !IS
2059 ENDDO !IR
2060 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
2061 ENDDO ! I
2062 ELSEIF(id ==3 ) THEN
2063 DO i=1,nel
2064 DO ir=1,nptr
2065 DO is=1,npts
2066 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
2067 DO k=1,3
2068 value(k) = value(k) + lbuf%SIGB(jj(k+6) + i )/npg
2069 ENDDO !k
2070 ENDDO !IS
2071 ENDDO !IR
2072 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
2073 ENDDO ! I
2074 ELSEIF(id ==4 ) THEN
2075 DO i=1,nel
2076 DO ir=1,nptr
2077 DO is=1,npts
2078 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
2079 DO k=1,3
2080 value(k) = value(k) + lbuf%SIGB(jj(k+9) + i )/npg
2081 ENDDO !k
2082 ENDDO !IS
2083 ENDDO !IR
2084 CALL h3d_write_sh_tensor(iok_part,iselect,is_written_shell,shell_tensor,i,offset,nft,VALUE)
2085 ENDDO ! I
2086 ENDIF !ID == -1
2087 ENDIF !(MLW ==
2088 endif! (IPT <= NPTT )
2089 ENDDO !JLAY
2090
2091 END IF
2092c ........
2093C---------------------------------------------------------------------------
2094c ELSEIF (KEYWORD == 'NEWKEY') THEN ! New Output Example
2095C---------------------------------------------------------------------------
2096c ILAYER=NULL NPT=NULL
2097c IF ( ILAY == -1 .AND. IPT == -1 .AND. IPLY == -1) THEN
2098c DO I=1,NEL
2099c VALUE(I) =
2100c ENDDO
2101c PLY=IPLY NPT=IPT
2102c ELSEIF ( IPLY > 0 .AND. IPT <= MPT .AND. IPT > 0 ) THEN
2103c IF (IGTYP == 17 .OR. IGTYP == 51 .OR. IGTYP == 52) THEN
2104c
2105c ENDIF
2106c
2107c PLY=NULL ILAYER=ILAY NPT=IPT
2108c ELSEIF (IPLY == -1 .AND. ILAY <= NLAY .AND. ILAY > 0 .AND. IPT <= MPT .AND. IPT > 0 ) THEN
2109c IF (IGTYP == 51 .OR. IGTYP == 52) THEN
2110c
2111c ENDIF
2112c PLY=NULL ILAYER=IL NPT=NULL
2113c ELSEIF (IPLY == -1 .AND. ILAY <= NLAY .AND. ILAY > 0 .AND. IPT == -1 ) THEN
2114c IF (IGTYP == 10 .OR. IGTYP == 11 .OR. IGTYP == 16 .OR. IGTYP == 17) THEN
2115c
2116c ELSEIF (IGTYP == 51 .OR. IGTYP == 52) THEN
2117c
2118c ENDIF
2119c PLY=NULL ILAYER=NULL NPT=IPT
2120c ELSEIF ( IPT <= MPT .AND. IPT > 0) THEN
2121c IF (IGTYP == 1 .OR. IGTYP == 9) THEN
2122c
2123c ENDIF
2124c ENDIF
2125 ENDIF
2126
2127C-----------------------------------------------
2128C RNUR
2129C-----------------------------------------------
2130 ELSEIF (ity == 50) THEN
2131c DO I=1,NEL
2132c N = I + NFT
2133c TENS(1,EL2FA(NN9+N)) = ZERO
2134c TENS(2,EL2FA(NN9+N)) = ZERO
2135c TENS(3,EL2FA(NN9+N)) = ZERO
2136c ENDDO
2137C-----------------------------------------------
2138 ELSE
2139 ENDIF ! IF (MLW /= 13)
2140 ENDIF ! IF(ITY == 2)
2141 490 CONTINUE
2142 500 CONTINUE
2143C-----------------------------------------------
2144C
2145 RETURN
2146 END
2147!||====================================================================
2148!|| sh4_tstrain ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2149!||--- called by ------------------------------------------------------
2150!|| h3d_shell_tensor ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2151!||--- calls -----------------------------------------------------
2152!|| c4sysg2l ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.f
2153!||====================================================================
2154 SUBROUTINE sh4_tstrain(XN,YN,ZN,DX,DY,DZ,STRAIN,NEL)
2155C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2156#include "implicit_f.inc"
2157C---------+---------+---+---+--------------------------------------------
2158C VAR | SIZE |TYP| RW| DEFINITION
2159C---------+---------+---+---+--------------------------------------------
2160C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
2161C XN | 4*NEL | R | R | X-coordinate ARRAY (4n quad)
2162C YN | 4*NEL | R | R | Y-coordinate ARRAY (4n quad)
2163C ZN | 4*NEL | R | R | Z-coordinate ARRAY (4n quad)
2164C DX | 4*NEL | R | R | X-Displ ARRAY (4n quad)
2165C DY | 4*NEL | R | R | Y-Displ ARRAY (4n quad)
2166C DZ | 4*NEL | R | R | D-Displ ARRAY (4n quad)
2167C STRAIN | 3*NEL | R | W | STRAIN ARRAY
2168C---------+---------+---+---+--------------------------------------------
2169C-----------------------------------------------
2170C D U M M Y A R G U M E N T S
2171C-----------------------------------------------
2172 INTEGER NEL
2173 my_real
2174 . XN(4,NEL) , YN(4,NEL) , ZN(4,NEL),
2175 . DX(4,NEL) , DY(4,NEL) , DZ(4,NEL),STRAIN(3,*)
2176C-----------------------------------------------
2177C L O C A L V A R I A B L E S
2178C-----------------------------------------------
2179 INTEGER I,NNOD
2180 PARAMETER (NNOD = 4)
2181 my_real
2182 . x0n(nnod,nel) , y0n(nnod,nel) , z0n(nnod,nel),
2183 . lxyz0(3),corel(2,nnod),
2184 . xl2(nel),xl3(nel),xl4(nel),yl2(nel),
2185 . yl3(nel),yl4(nel),zl1(nel),zl(nel),
2186 . x13,x24,y13,y24,ux1,uy1,ux2,ux3,ux4,uy2,uy3,uy4,
2187 . px1(nel),px2(nel),py1(nel),py2(nel),
2188 . ux13(nel),ux24(nel),uy13(nel),uy24(nel),
2189 . x0l2(nel),x0l3(nel),x0l4(nel),y0l2(nel),
2190 . y0l3(nel),y0l4(nel),area(nel),area_i(nel),fxx,fyy,fxy,fyx
2191C----------------------------------------------
2192C------Compute coordinates in elementary local sys: actual configuration first
2193 CALL c4sysg2l(xn,yn,zn,xl2,yl2,xl3,yl3,xl4,yl4,zl1,area,nel)
2194C------initial configuration :
2195 DO i=1,nel
2196 x0n(1:nnod,i) = xn(1:nnod,i)-dx(1:nnod,i)
2197 y0n(1:nnod,i) = yn(1:nnod,i)-dy(1:nnod,i)
2198 z0n(1:nnod,i) = zn(1:nnod,i)-dz(1:nnod,i)
2199 ENDDO
2200 CALL c4sysg2l(x0n,y0n,z0n,x0l2,y0l2,x0l3,y0l3,x0l4,y0l4,zl,area,nel)
2201C-------[B0]---remove the origine to the center--------------
2202 DO i=1,nel
2203C-----------EM20=1.0E-20,FOURTH=0.25,HALF=0.5,ZERO=0.--------------
2204 area_i(i) = one/max(em20,area(i))
2205 lxyz0(1)=fourth*(x0l2(i)+x0l3(i)+x0l4(i))
2206 lxyz0(2)=fourth*(y0l2(i)+y0l3(i)+y0l4(i))
2207 corel(1,1)=-lxyz0(1)
2208 corel(1,2)=x0l2(i)-lxyz0(1)
2209 corel(1,3)=x0l3(i)-lxyz0(1)
2210 corel(1,4)=x0l4(i)-lxyz0(1)
2211 corel(2,1)=-lxyz0(2)
2212 corel(2,2)=y0l2(i)-lxyz0(2)
2213 corel(2,3)=y0l3(i)-lxyz0(2)
2214 corel(2,4)=y0l4(i)-lxyz0(2)
2215C----
2216 x13 =(corel(1,1)-corel(1,3))*half
2217 x24 =(corel(1,2)-corel(1,4))*half
2218 y13 =(corel(2,1)-corel(2,3))*half
2219 y24 =(corel(2,2)-corel(2,4))*half
2220 py2(i) =x13*area_i(i)
2221 py1(i) =-x24*area_i(i)
2222 px2(i) =-y13*area_i(i)
2223 px1(i) =y24*area_i(i)
2224 ENDDO
2225C------ objective disp (or projected disp to free rigide rotation)
2226 ux1= zero
2227 uy1= zero
2228 DO i=1,nel
2229 ux2 = xl2(i)-x0l2(i)
2230 uy2 = yl2(i)-y0l2(i)
2231 ux3 = xl3(i)-x0l3(i)
2232 uy3 = yl3(i)-y0l3(i)
2233 ux4 = xl4(i)-x0l4(i)
2234 uy4 = yl4(i)-y0l4(i)
2235 ux13(i)=ux1-ux3
2236 ux24(i)=ux2-ux4
2237 uy13(i)=uy1-uy3
2238 uy24(i)=uy2-uy4
2239 ENDDO
2240C---------------
2241C MEMBRANE [F]-1 = [B0]{d}, [e] ({d} has been projected to free rigide rotation)
2242C---------------
2243 DO i=1,nel
2244 fxx = px1(i)*ux13(i)+px2(i)*ux24(i)
2245 fyy = py1(i)*uy13(i)+py2(i)*uy24(i)
2246 fyx = px1(i)*uy13(i)+px2(i)*uy24(i)
2247 fxy = py1(i)*ux13(i)+py2(i)*ux24(i)
2248 strain(1,i) = fxx
2249 strain(2,i) = fyy
2250 strain(3,i) = 0.5*(fxy+fyx)
2251 ENDDO
2252C
2253 RETURN
2254 END
2255!||====================================================================
2256!|| c4sysg2l ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2257!||--- called by ------------------------------------------------------
2258!|| sh4_tstrain ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2259!||--- calls -----------------------------------------------------
2260!|| clsys3 ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2261!||====================================================================
2262 SUBROUTINE c4sysg2l(XN,YN,ZN,XL2,YL2,XL3,YL3,XL4,YL4,ZL1,AREA,NEL)
2263C-----------------------------------------------
2264C I m p l i c i t T y p e s
2265C-----------------------------------------------
2266#include "implicit_f.inc"
2267C---------+---------+---+---+--------------------------------------------
2268C VAR | SIZE |TYP| RW| DEFINITION
2269C---------+---------+---+---+--------------------------------------------
2270C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
2271C XN | 4*NEL | R | R | X-coordinate ARRAY (4n quad)
2272C YN | 4*NEL | R | R | Y-coordinate ARRAY (4n quad)
2273C ZN | 4*NEL | R | R | Z-coordinate ARRAY (4n quad)
2274C XL2 | NEL | R | W | Local X-coordinate of N2 (relative to N1)
2275C YL2 | NEL | R | W | Local Y-coordinate of N2 (relative to N1)
2276C ZL1 | NEL | R | W | Local Z-coordinate of N1 (Z1 -Z1 Z1 -Z1)
2277C XL3 | NEL | R | W | Local X-coordinate of N3 (relative to N1)
2278C YL3 | NEL | R | W | Local Y-coordinate of N3 (relative to N1)
2279C XL4 | NEL | R | W | Local X-coordinate of N4 (relative to N1)
2280C YL4 | NEL | R | W | Local Y-coordinate of N4 (relative to N1)
2281C AREA | NEL | R | W | AREA of quad
2282C---------+---------+---+---+--------------------------------------------
2283C-----------------------------------------------
2284C D u m m y A r g u m e n t s
2285C-----------------------------------------------
2286 INTEGER NEL
2287 my_real
2288 . XN(4,*) , YN(4,*) , ZN(4,*),
2289 . XL2(NEL),XL3(NEL),XL4(NEL),YL2(NEL),
2290 . YL3(NEL),YL4(NEL),ZL1(NEL),AREA(NEL)
2291C-----------------------------------------------
2292C L o c a l V a r i a b l e s
2293C-----------------------------------------------
2294 INTEGER I
2295 my_real
2296 . RX(NEL),RY(NEL),RZ(NEL),SX(NEL),SY(NEL),SZ(NEL),
2297 . VQ(3,3,NEL), LXYZ0(3),DETA1(NEL),XX,YY,ZZ
2298C-----------------------------------------------
2299 DO I=1,nel
2300 rx(i)=xn(2,i)+xn(3,i)-xn(1,i)-xn(4,i)
2301 ry(i)=yn(2,i)+yn(3,i)-yn(1,i)-yn(4,i)
2302 rz(i)=zn(2,i)+zn(3,i)-zn(1,i)-zn(4,i)
2303 sx(i)=xn(3,i)+xn(4,i)-xn(1,i)-xn(2,i)
2304 sy(i)=yn(3,i)+yn(4,i)-yn(1,i)-yn(2,i)
2305 sz(i)=zn(3,i)+zn(4,i)-zn(1,i)-zn(2,i)
2306 ENDDO
2307C------Local elem. base:
2308 CALL clsys3(rx, ry, rz, sx, sy, sz,
2309 . vq, deta1,nel)
2310C------ Global -> Local Coordinate FOURTH=0.25 ;
2311 DO i=1,nel
2312 lxyz0(1)=fourth*(xn(1,i)+xn(2,i)+xn(3,i)+xn(4,i))
2313 lxyz0(2)=fourth*(yn(1,i)+yn(2,i)+yn(3,i)+yn(4,i))
2314 lxyz0(3)=fourth*(zn(1,i)+zn(2,i)+zn(3,i)+zn(4,i))
2315 xx=xn(2,i)-xn(1,i)
2316 yy=yn(2,i)-yn(1,i)
2317 zz=zn(2,i)-zn(1,i)
2318 xl2(i)=vq(1,1,i)*xx+vq(2,1,i)*yy+vq(3,1,i)*zz
2319 yl2(i)=vq(1,2,i)*xx+vq(2,2,i)*yy+vq(3,2,i)*zz
2320 xx=xn(2,i)-lxyz0(1)
2321 yy=yn(2,i)-lxyz0(2)
2322 zz=zn(2,i)-lxyz0(3)
2323 zl1(i)=vq(1,3,i)*xx+vq(2,3,i)*yy+vq(3,3,i)*zz
2324C
2325 xx=xn(3,i)-xn(1,i)
2326 yy=yn(3,i)-yn(1,i)
2327 zz=zn(3,i)-zn(1,i)
2328 xl3(i)=vq(1,1,i)*xx+vq(2,1,i)*yy+vq(3,1,i)*zz
2329 yl3(i)=vq(1,2,i)*xx+vq(2,2,i)*yy+vq(3,2,i)*zz
2330C
2331 xx=xn(4,i)-xn(1,i)
2332 yy=yn(4,i)-yn(1,i)
2333 zz=zn(4,i)-zn(1,i)
2334 xl4(i)=vq(1,1,i)*xx+vq(2,1,i)*yy+vq(3,1,i)*zz
2335 yl4(i)=vq(1,2,i)*xx+vq(2,2,i)*yy+vq(3,2,i)*zz
2336 area(i)=fourth*deta1(i)
2337 ENDDO
2338c-----------
2339 RETURN
2340 END
2341!||====================================================================
2342!|| clsys3 ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2343!||--- called by ------------------------------------------------------
2344!|| c3sysg2l ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2345!|| c4sysg2l ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2346!|| sdlensh ../engine/source/elements/thickshell/solidec/sdlensh.F
2347!|| sdlensh8 ../engine/source/elements/thickshell/solide8c/sdlensh8.F
2348!||====================================================================
2349 SUBROUTINE clsys3(RX, RY, RZ, SX, SY, SZ, VQ, DET,NEL)
2350C-----------------------------------------------
2351C I m p l i c i t T y p e s
2352C-----------------------------------------------
2353#include "implicit_f.inc"
2354C-----------------------------------------------
2355C D u m m y A r g u m e n t s
2356C-----------------------------------------------
2357 INTEGER NEL
2358 my_real
2359 . RX(*) , RY(*) , RZ(*),
2360 . SX(*) , SY(*) , SZ(*),
2361 . DET(*),VQ(3,3,*)
2362C---------+---------+---+---+--------------------------------------------
2363C VAR | SIZE |TYP| RW| DEFINITION
2364C---------+---------+---+---+--------------------------------------------
2365C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
2366C RX | NEL | R | R | X-of covariant vecter g1
2367C RY | NEL | R | R | Y-of covariant vecter g1
2368C RZ | NEL | R | R | Z-of covariant vecter g1
2369C SX | NEL | R | R | X-of covariant vecter g2
2370C SY | NEL | R | R | Y-of covariant vecter g2
2371C SZ | NEL | R | R | Z-of covariant vecter g2
2372C VQ |3*3*NEL | R | W | Local elem sys bases
2373C DET | NEL | R | W | det of g1 ^ g2
2374C---------+---------+---+---+--------------------------------------------
2375C-----------------------------------------------
2376C L o c a l V a r i a b l e s
2377C-----------------------------------------------
2378 INTEGER I
2379C REAL
2380 my_real
2381 . E1X(NEL), E1Y(NEL), E1Z(NEL),
2382 . E2X(NEL), E2Y(NEL), E2Z(NEL),
2383 . E3X(NEL), E3Y(NEL), E3Z(NEL),
2384 . c1,c1c1,c2c2,c1_1,c2_1
2385C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2386 DO i=1,nel
2387C---------E3------------
2388 e3x(i) = ry(i) * sz(i) - rz(i) * sy(i)
2389 e3y(i) = rz(i) * sx(i) - rx(i) * sz(i)
2390 e3z(i) = rx(i) * sy(i) - ry(i) * sx(i)
2391 det(i) = sqrt(e3x(i)*e3x(i) + e3y(i)*e3y(i) + e3z(i)*e3z(i))
2392C ----- EM20=1.0E-20
2393 det(i) = max(em20,det(i))
2394 e3x(i) = e3x(i) / det(i)
2395 e3y(i) = e3y(i) / det(i)
2396 e3z(i) = e3z(i) / det(i)
2397 ENDDO
2398C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2399 DO i=1,nel
2400 c1c1 = rx(i)*rx(i) + ry(i)*ry(i) + rz(i)*rz(i)
2401 c2c2 = sx(i)*sx(i) + sy(i)*sy(i) + sz(i)*sz(i)
2402C ----- ZERO=0., ONE=1.0
2403 IF(c1c1 /= zero) THEN
2404 c2_1 = sqrt(c2c2/max(em20,c1c1))
2405 c1_1 = one
2406 ELSEIF(c2c2 /= zero)THEN
2407 c2_1 = one
2408 c1_1 = sqrt(c1c1/max(em20,c2c2))
2409 END IF
2410 e1x(i) = rx(i)*c2_1+(sy(i)*e3z(i)-sz(i)*e3y(i))*c1_1
2411 e1y(i) = ry(i)*c2_1+(sz(i)*e3x(i)-sx(i)*e3z(i))*c1_1
2412 e1z(i) = rz(i)*c2_1+(sx(i)*e3y(i)-sy(i)*e3x(i))*c1_1
2413 ENDDO
2414C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2415 DO i=1,nel
2416 c1 = sqrt(e1x(i)*e1x(i) + e1y(i)*e1y(i) + e1z(i)*e1z(i))
2417 IF ( c1 /= zero) c1 = one / max(em20,c1)
2418 e1x(i) = e1x(i)*c1
2419 e1y(i) = e1y(i)*c1
2420 e1z(i) = e1z(i)*c1
2421 e2x(i) = e3y(i) * e1z(i) - e3z(i) * e1y(i)
2422 e2y(i) = e3z(i) * e1x(i) - e3x(i) * e1z(i)
2423 e2z(i) = e3x(i) * e1y(i) - e3y(i) * e1x(i)
2424 ENDDO
2425 DO i=1,nel
2426 vq(1,1,i)=e1x(i)
2427 vq(2,1,i)=e1y(i)
2428 vq(3,1,i)=e1z(i)
2429 vq(1,2,i)=e2x(i)
2430 vq(2,2,i)=e2y(i)
2431 vq(3,2,i)=e2z(i)
2432 vq(1,3,i)=e3x(i)
2433 vq(2,3,i)=e3y(i)
2434 vq(3,3,i)=e3z(i)
2435 ENDDO
2436c-----------
2437 RETURN
2438 END
2439!||====================================================================
2440!|| sh3_tstrain ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2441!||--- called by ------------------------------------------------------
2442!|| h3d_shell_tensor ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2443!||--- calls -----------------------------------------------------
2444!|| c3sysg2l ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2445!|| u_from_f2 ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2446!||====================================================================
2447 SUBROUTINE sh3_tstrain(XN,YN,ZN,DX,DY,DZ,STRAIN,NEL,IHBE)
2448C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2449#include "implicit_f.inc"
2450C---------+---------+---+---+--------------------------------------------
2451C VAR | SIZE |TYP| RW| DEFINITION
2452C---------+---------+---+---+--------------------------------------------
2453C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
2454C XN | 3*NEL | R | R | X-coordinate ARRAY (3n tria)
2455C YN | 3*NEL | R | R | Y-coordinate ARRAY (3n tria)
2456C ZN | 3*NEL | R | R | Z-coordinate ARRAY (3n tria)
2457C DX | 3*NEL | R | R | X-Displ ARRAY (3n tria)
2458C DY | 3*NEL | R | R | Y-Displ ARRAY (3n tria)
2459C DZ | 3*NEL | R | R | D-Displ ARRAY (3n tria)
2460C STRAIN | 3*NEL | R | W | STRAIN ARRAY
2461C---------+---------+---+---+--------------------------------------------
2462C-----------------------------------------------
2463C D U M M Y A R G U M E N T S
2464C-----------------------------------------------
2465 INTEGER NEL,IHBE
2466 my_real
2467 . XN(3,NEL) , YN(3,NEL) , ZN(3,NEL),
2468 . DX(3,NEL) , DY(3,NEL) , DZ(3,NEL),STRAIN(3,*)
2469C-----------------------------------------------
2470C L O C A L V A R I A B L E S
2471C-----------------------------------------------
2472 INTEGER I,NNOD
2473 PARAMETER (NNOD = 3)
2474 my_real
2475 . x0n(nnod,nel) , y0n(nnod,nel) , z0n(nnod,nel),
2476 .
2477 . xl2(nel),xl3(nel),yl2(nel),yl3(nel),
2478 . ux1,uy1,ux2,ux3,uy2,uy3,
2479 . px2(nel),px3(nel),py2(nel),py3(nel),
2480 . ux21(nel),ux31(nel),uy21(nel),uy31(nel),
2481 . x0l2(nel),x0l3(nel),y0l2(nel),
2482 . y0l3(nel),area(nel),area_i(nel),f(2,2,nel)
2483C----------------------------------------------
2484C------Compute coordinates in elementary local sys: actual configuration first
2485 CALL c3sysg2l(xn,yn,zn,xl2,yl2,xl3,yl3,area,nel)
2486C------initial configuration :
2487 DO i=1,nel
2488 x0n(1:nnod,i) = xn(1:nnod,i)-dx(1:nnod,i)
2489 y0n(1:nnod,i) = yn(1:nnod,i)-dy(1:nnod,i)
2490 z0n(1:nnod,i) = zn(1:nnod,i)-dz(1:nnod,i)
2491 ENDDO
2492 CALL c3sysg2l(x0n,y0n,z0n,x0l2,y0l2,x0l3,y0l3,area,nel)
2493C-----------[B0]-----------------
2494 DO i=1,nel
2495 area_i(i)=0.5/area(i)
2496 px2(i)= y0l3(i)*area_i(i)
2497 py2(i)=-x0l3(i)*area_i(i)
2498 px3(i)=-y0l2(i)*area_i(i)
2499 py3(i)= x0l2(i)*area_i(i)
2500 ENDDO
2501C------ objective disp (free rigide rotation)
2502 ux1= zero
2503 uy1= zero
2504 DO i=1,nel
2505 ux2 = xl2(i)-x0l2(i)
2506 uy2 = yl2(i)-y0l2(i)
2507 ux3 = xl3(i)-x0l3(i)
2508 uy3 = yl3(i)-y0l3(i)
2509 ux21(i)=ux2-ux1
2510 ux31(i)=ux3-ux1
2511 uy21(i)=uy2-uy1
2512 uy31(i)=uy3-uy1
2513 ENDDO
2514C---------------
2515C MEMBRANE [F]-1 = [B0]{d}, [e] ({d} has been projected to free rigide rotation)
2516 DO i=1,nel
2517 f(1,1,i) = px2(i)*ux21(i)+px3(i)*ux31(i)
2518 f(2,2,i) = py2(i)*uy21(i)+py3(i)*uy31(i)
2519 f(2,1,i) = px2(i)*uy21(i)+px3(i)*uy31(i)
2520 f(1,2,i) = py2(i)*ux21(i)+py3(i)*ux31(i)
2521 strain(1,i) = f(1,1,i)
2522 strain(2,i) = f(2,2,i)
2523 strain(3,i) = 0.5*(f(1,2,i) + f(2,1,i))
2524 ENDDO
2525C--- local sys for 3N isn't really material ----
2526 IF (ihbe==3) CALL u_from_f2(f,strain,nel )
2527C
2528 RETURN
2529 END
2530!||====================================================================
2531!|| c3sysg2l ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2532!||--- called by ------------------------------------------------------
2533!|| sh3_tstrain ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2534!||--- calls -----------------------------------------------------
2535!|| clsys3 ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2536!||====================================================================
2537 SUBROUTINE c3sysg2l(XN,YN,ZN,XL2,YL2,XL3,YL3,AREA,NEL)
2538C-----------------------------------------------
2539C I m p l i c i t T y p e s
2540C-----------------------------------------------
2541#include "implicit_f.inc"
2542C---------+---------+---+---+--------------------------------------------
2543C VAR | SIZE |TYP| RW| DEFINITION
2544C---------+---------+---+---+--------------------------------------------
2545C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
2546C XN | 3*NEL | R | R | X-coordinate ARRAY (3n tria)
2547C YN | 3*NEL | R | R | Y-coordinate ARRAY (3n tria)
2548C ZN | 3*NEL | R | R | Z-coordinate ARRAY (3n tria)
2549C XL2 | NEL | R | W | Local X-coordinate of N2 (relative to N1)
2550C YL2 | NEL | R | W | Local Y-coordinate of N2 (relative to N1)
2551C XL3 | NEL | R | W | Local X-coordinate of N3 (relative to N1)
2552C YL3 | NEL | R | W | Local Y-coordinate of N3 (relative to N1)
2553C AREA | NEL | R | W | AREA of tria
2554C---------+---------+---+---+--------------------------------------------
2555C-----------------------------------------------
2556C D u m m y A r g u m e n t s
2557C-----------------------------------------------
2558 INTEGER NEL
2559 my_real
2560 . XN(3,*) , YN(3,*) , ZN(3,*),
2561 . XL2(NEL),XL3(NEL),YL2(NEL),YL3(NEL),AREA(NEL)
2562C-----------------------------------------------
2563C L o c a l V a r i a b l e s
2564C-----------------------------------------------
2565 INTEGER I
2566 my_real
2567 . RX(NEL),RY(NEL),RZ(NEL),SX(NEL),SY(NEL),SZ(NEL),
2568 . VQ(3,3,NEL), DETA1(NEL),XX,YY,ZZ
2569 DO I=1,nel
2570 rx(i)=xn(2,i)-xn(1,i)
2571 ry(i)=yn(2,i)-yn(1,i)
2572 rz(i)=zn(2,i)-zn(1,i)
2573 sx(i)=xn(3,i)-xn(1,i)
2574 sy(i)=yn(3,i)-yn(1,i)
2575 sz(i)=zn(3,i)-zn(1,i)
2576 ENDDO
2577 CALL clsys3(rx, ry, rz, sx, sy, sz,
2578 . vq, deta1,nel)
2579C------ Global -> Local Coordinate
2580 DO i=1,nel
2581 xx=rx(i)
2582 yy=ry(i)
2583 zz=rz(i)
2584 xl2(i)=vq(1,1,i)*xx+vq(2,1,i)*yy+vq(3,1,i)*zz
2585 yl2(i)=vq(1,2,i)*xx+vq(2,2,i)*yy+vq(3,2,i)*zz
2586C
2587 xx=sx(i)
2588 yy=sy(i)
2589 zz=sz(i)
2590 xl3(i)=vq(1,1,i)*xx+vq(2,1,i)*yy+vq(3,1,i)*zz
2591 yl3(i)=vq(1,2,i)*xx+vq(2,2,i)*yy+vq(3,2,i)*zz
2592 area(i)=0.5*deta1(i)
2593 ENDDO
2594c-----------
2595 RETURN
2596 END
2597!||====================================================================
2598!|| u_from_f2 ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2599!||--- called by ------------------------------------------------------
2600!|| sh3_tstrain ../engine/source/output/h3d/h3d_results/h3d_shell_tensor.F
2601!||====================================================================
2602 SUBROUTINE u_from_f2(F,STRAIN,NEL )
2603C-----------------------------------------------
2604C I m p l i c i t T y p e s
2605C-----------------------------------------------
2606#include "implicit_f.inc"
2607C-----------------------------------------------
2608C D u m m y A r g u m e n t s
2609C-----------------------------------------------
2610 INTEGER NEL
2611 my_real
2612 . F(2,2,NEL), STRAIN(3,*)
2613C-----------------------------------------------
2614C L o c a l V a r i a b l e s
2615C-----------------------------------------------
2616 INTEGER I
2617 DOUBLE PRECISION
2618 . FMAT(2,2),UM(2,2),IC,
2619 . C11,C12,C22,
2620 . CC11,CC12,CC22,
2621 . B,ALPHA,C_A,S_A,C_A2,S_A2
2622C-----------------------------------------------
2623 DO i=1,nel
2624 fmat(1,1)= f(1,1,i)+1.0
2625 fmat(2,2)= f(2,2,i)+1.0
2626 fmat(1,2)= f(1,2,i)
2627 fmat(2,1)= f(2,1,i)
2628 c11 = fmat(1,1)*fmat(1,1)+fmat(2,1)*fmat(2,1)
2629 c12 = fmat(1,1)*fmat(1,2)+fmat(2,1)*fmat(2,2)
2630 c22 = fmat(1,2)*fmat(1,2)+fmat(2,2)*fmat(2,2)
2631 cc12 = 0.5*(c11-c22)
2632 ic = 0.5*(c11+c22)
2633 b = sqrt(cc12*cc12+c12*c12)
2634 cc11 = sqrt(ic + b)
2635 cc22 = sqrt(ic - b)
2636 IF (abs(cc12)<em20) THEN
2637 c_a = 0.0
2638 s_a = 1.0
2639 ELSE
2640 alpha =0.5*(atan(c12/cc12))
2641 c_a = cos(alpha)
2642 s_a = sin(alpha)
2643 END IF
2644 c_a2 = c_a*c_a
2645 s_a2 = s_a*s_a
2646 um(1,1) = cc11*c_a2+cc22*s_a2
2647 um(2,2) = cc11*s_a2+cc22*c_a2
2648 um(1,2) = (cc11-cc22)*s_a*c_a
2649 strain(1,i) = um(1,1)-one
2650 strain(2,i) = um(2,2)-one
2651 strain(3,i) = um(1,2)
2652 END DO
2653C
2654 RETURN
2655 END
2656
my_real
#define my_real
Definition cppsort.cpp:32
u_from_f2
subroutine u_from_f2(f, strain, nel)
Definition h3d_shell_tensor.F:2603
c3sysg2l
subroutine c3sysg2l(xn, yn, zn, xl2, yl2, xl3, yl3, area, nel)
Definition h3d_shell_tensor.F:2538
sh3_tstrain
subroutine sh3_tstrain(xn, yn, zn, dx, dy, dz, strain, nel, ihbe)
Definition h3d_shell_tensor.F:2448
clsys3
subroutine clsys3(rx, ry, rz, sx, sy, sz, vq, det, nel)
Definition h3d_shell_tensor.F:2350
sh4_tstrain
subroutine sh4_tstrain(xn, yn, zn, dx, dy, dz, strain, nel)
Definition h3d_shell_tensor.F:2155
h3d_shell_tensor
subroutine h3d_shell_tensor(elbuf_tab, shell_tensor, iparg, itens, invert, nelcut, el2fa, nbf, tens, epsdot, iadp, nbf_l, nbpart, iadg, x, ixc, igeo, ixtg, ipm, stack, id_elem, ity_elem, info1, info2, is_written_shell, ipartc, iparttg, layer_input, ipt_input, ply_input, gauss_input, iuvar_input, h3d_part, keyword, d, id, bufmat, mat_param, geo, drape_sh4n, drape_sh3n, drapeg)
Definition h3d_shell_tensor.F:52
c4sysg2l
subroutine c4sysg2l(xn, yn, zn, xl2, yl2, xl3, yl3, xl4, yl4, zl1, area, nel)
Definition h3d_shell_tensor.F:2263
h3d_write_sh_tensor
subroutine h3d_write_sh_tensor(iok_part, iselect, is_written, tensor, i, offset, nft, value)
Definition h3d_write_sh_tensor.F:31
h3d_write_sh_tensor_array
subroutine h3d_write_sh_tensor_array(iok_part, iselect, nel, offset, nft, is_written, tensor, value)
Definition h3d_write_sh_tensor_array.F:31
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
layini
subroutine layini(elbuf_str, jft, jlt, geo, igeo, mat, pid, thkly, matly, posly, igtyp, ixfem, ixlay, nlay, npt, isubstack, stack, drape, nft, thk, nel, ratio_thkly, indx_drape, sedrape, numel_drape)
Definition layini.F:47
max
#define max(a, b)
Definition macros.h:21
drape_mod
Definition drape_mod.F:54
drape_mod::scdrape
integer scdrape
Definition drape_mod.F:92
drape_mod::stdrape
integer stdrape
Definition drape_mod.F:92
names_and_titles_mod
Definition names_and_titles_mod.F:1002
names_and_titles_mod::ncharline100
integer, parameter ncharline100
Definition names_and_titles_mod.F:1007
stack_mod
Definition stack_mod.F:75
stack_mod::ply_info
integer, dimension(:,:), allocatable ply_info
Definition stack_mod.F:133
roto_tens2d
subroutine roto_tens2d(nel, sig, dir)
Definition roto_tens2d.F:29
roto_tens2d_aniso
subroutine roto_tens2d_aniso(nel, tens, dir_a, dir_b)
Definition roto_tens2d_aniso.F:29
uroto_tens2d
subroutine uroto_tens2d(nel, sig, dir)
Definition uroto_tens2d.F:30
uroto_tens2d_aniso
subroutine uroto_tens2d_aniso(nel, tens, dir_a, dir_b)
Definition uroto_tens2d_aniso.F:30