OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
h3d_shell_tensor.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "param_c.inc"
#include "tabsiz_c.inc"

Go to the source code of this file.

Functions/Subroutines

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)
subroutine sh4_tstrain (xn, yn, zn, dx, dy, dz, strain, nel)
subroutine c4sysg2l (xn, yn, zn, xl2, yl2, xl3, yl3, xl4, yl4, zl1, area, nel)
subroutine clsys3 (rx, ry, rz, sx, sy, sz, vq, det, nel)
subroutine sh3_tstrain (xn, yn, zn, dx, dy, dz, strain, nel, ihbe)
subroutine c3sysg2l (xn, yn, zn, xl2, yl2, xl3, yl3, area, nel)
subroutine u_from_f2 (f, strain, nel)

Function/Subroutine Documentation

◆ c3sysg2l()

subroutine c3sysg2l ( xn,
yn,
zn,
xl2,
yl2,
xl3,
yl3,
area,
integer nel )

Definition at line 2537 of file h3d_shell_tensor.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
clsys3
subroutine clsys3(rx, ry, rz, sx, sy, sz, vq, det, nel)
Definition h3d_shell_tensor.F:2350
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567

◆ c4sysg2l()

subroutine c4sysg2l ( xn,
yn,
zn,
xl2,
yl2,
xl3,
yl3,
xl4,
yl4,
zl1,
area,
integer nel )

Definition at line 2262 of file h3d_shell_tensor.F.

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

◆ clsys3()

subroutine clsys3 ( rx,
ry,
rz,
sx,
sy,
sz,
vq,
det,
integer nel )

Definition at line 2349 of file h3d_shell_tensor.F.

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
max
#define max(a, b)
Definition macros.h:21

◆ h3d_shell_tensor()

subroutine h3d_shell_tensor ( type (elbuf_struct_), dimension(ngroup), target elbuf_tab,
shell_tensor,
integer, dimension(nparg,*) iparg,
integer itens,
integer, dimension(*) invert,
integer nelcut,
integer, dimension(*) el2fa,
integer nbf,
tens,
epsdot,
integer, dimension(*) iadp,
integer nbf_l,
integer nbpart,
integer, dimension(nspmd,*) iadg,
x,
integer, dimension(nixc,*) ixc,
integer, dimension(npropgi,*) igeo,
integer, dimension(nixtg,*) ixtg,
integer, dimension(npropmi,*) ipm,
type (stack_ply) stack,
integer, dimension(*) id_elem,
integer, dimension(*) ity_elem,
integer info1,
integer info2,
integer, dimension(*) is_written_shell,
integer, dimension(*) ipartc,
integer, dimension(*) iparttg,
integer layer_input,
integer ipt_input,
integer ply_input,
integer gauss_input,
integer iuvar_input,
integer, dimension(*) h3d_part,
character(len=ncharline100) keyword,
d,
integer id,
dimension(sbufmat), target bufmat,
type (matparam_struct_), dimension(nummat), intent(in) mat_param,
dimension(npropg,numgeo), intent(in) geo,
type (drape_), dimension(numelc_drape), intent(in) drape_sh4n,
type (drape_), dimension(numeltg_drape), intent(in) drape_sh3n,
type (drapeg_), intent(in) drapeg )

Definition at line 45 of file h3d_shell_tensor.F.

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
sh3_tstrain
subroutine sh3_tstrain(xn, yn, zn, dx, dy, dz, strain, nel, ihbe)
Definition h3d_shell_tensor.F:2448
sh4_tstrain
subroutine sh4_tstrain(xn, yn, zn, dx, dy, dz, strain, nel)
Definition h3d_shell_tensor.F:2155
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
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
id
initmumps id
Definition multiplerhs_example.m:4
drape_mod
Definition drape_mod.F:54
drape_mod::numeltg_drape
integer numeltg_drape
Definition drape_mod.F:92
drape_mod::scdrape
integer scdrape
Definition drape_mod.F:92
drape_mod::stdrape
integer stdrape
Definition drape_mod.F:92
drape_mod::numelc_drape
integer numelc_drape
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
nbfunct
subroutine nbfunct(nfunct, ntable, npts, lsubmodel)
Definition nbfunc.F:37
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

◆ sh3_tstrain()

subroutine sh3_tstrain ( xn,
yn,
zn,
dx,
dy,
dz,
strain,
integer nel,
integer ihbe )

Definition at line 2447 of file h3d_shell_tensor.F.

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
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

◆ sh4_tstrain()

subroutine sh4_tstrain ( xn,
yn,
zn,
dx,
dy,
dz,
strain,
integer nel )

Definition at line 2154 of file h3d_shell_tensor.F.

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
c4sysg2l
subroutine c4sysg2l(xn, yn, zn, xl2, yl2, xl3, yl3, xl4, yl4, zl1, area, nel)
Definition h3d_shell_tensor.F:2263

◆ u_from_f2()

subroutine u_from_f2 ( f,
strain,
integer nel )

Definition at line 2602 of file h3d_shell_tensor.F.

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
alpha
#define alpha
Definition eval.h:35