OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dfuncc.F File Reference
#include "implicit_f.inc"
#include "vect01_c.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "scr25_c.inc"
#include "param_c.inc"
#include "task_c.inc"
#include "spmd_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine dfuncc (elbuf_tab, func, ifunc, iparg, geo, ixq, ixc, ixtg, mass, pm, el2fa, nbf, iadp, itherm, nbf_l, ehour, anim, nbpart, iadg, ipm, igeo, thke, err_thk_sh4, err_thk_sh3, invert, x, v, w, ale_connectivity, nv46, nercvois, nesdvois, lercvois, lesdvois, stack, bufmat, multi_fvm, mat_param)
subroutine clsconv3 (rx, ry, rz, sx, sy, sz, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)

Function/Subroutine Documentation

◆ clsconv3()

subroutine clsconv3 ( rx,
ry,
rz,
sx,
sy,
sz,
e1x,
e1y,
e1z,
e2x,
e2y,
e2z,
e3x,
e3y,
e3z )

Definition at line 5191 of file dfuncc.F.

5195C-----------------------------------------------
5196C I m p l i c i t T y p e s
5197C-----------------------------------------------
5198#include "implicit_f.inc"
5199C-----------------------------------------------
5200C D u m m y A r g u m e n t s
5201C-----------------------------------------------
5202
5203 my_real
5204 . rx , ry , rz, sx , sy , sz,
5205 . e1x, e1y, e1z,e2x, e2y, e2z,e3x, e3y, e3z
5206C-----------------------------------------------
5207C L o c a l V a r i a b l e s
5208C-----------------------------------------------
5209
5210 my_real
5211 . c1,c2,cc,det
5212C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5213C---------E3------------
5214 e3x = ry * sz - rz * sy
5215 e3y = rz * sx - rx * sz
5216 e3z = rx * sy - ry * sx
5217 det = sqrt(e3x*e3x + e3y*e3y + e3z*e3z)
5218 det =max(em20,det)
5219 cc = one / det
5220 e3x = e3x * cc
5221 e3y = e3y * cc
5222 e3z = e3z * cc
5223C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5224 c1 = sqrt(rx*rx + ry*ry + rz*rz)
5225 c2 = sqrt(sx*sx + sy*sy + sz*sz)
5226 e1x = rx*c2+(sy*e3z-sz*e3y)*c1
5227 e1y = ry*c2+(sz*e3x-sx*e3z)*c1
5228 e1z = rz*c2+(sx*e3y-sy*e3x)*c1
5229C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5230 c1 = sqrt(e1x*e1x + e1y*e1y + e1z*e1z)
5231 IF ( c1/=zero) c1 = one / c1
5232 e1x = e1x*c1
5233 e1y = e1y*c1
5234 e1z = e1z*c1
5235 e2x = e3y * e1z - e3z * e1y
5236 e2y = e3z * e1x - e3x * e1z
5237 e2z = e3x * e1y - e3y * e1x
5238C
5239 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21

◆ dfuncc()

subroutine dfuncc ( type (elbuf_struct_), dimension(ngroup), target elbuf_tab,
func,
integer ifunc,
integer, dimension(nparg,ngroup) iparg,
geo,
integer, dimension(nixq,numelq) ixq,
integer, dimension(nixc,numelc) ixc,
integer, dimension(nixtg,numeltg) ixtg,
mass,
pm,
integer, dimension(*) el2fa,
integer nbf,
integer, dimension(*) iadp,
integer, intent(in) itherm,
integer nbf_l,
ehour,
anim,
integer nbpart,
integer, dimension(nspmd,*) iadg,
integer, dimension(npropmi,nummat) ipm,
integer, dimension(npropgi,numgeo) igeo,
thke,
err_thk_sh4,
err_thk_sh3,
integer, dimension(*) invert,
x,
v,
w,
type(t_ale_connectivity), intent(in) ale_connectivity,
integer nv46,
integer, dimension(*) nercvois,
integer, dimension(*) nesdvois,
integer, dimension(*) lercvois,
integer, dimension(*) lesdvois,
type (stack_ply) stack,
target bufmat,
type(multi_fvm_struct), intent(in) multi_fvm,
type (matparam_struct_), dimension(nummat), intent(in) mat_param )

Definition at line 47 of file dfuncc.F.

55C-----------------------------------------------
56C M o d u l e s
57C-----------------------------------------------
58 USE initbuf_mod
59 USE elbufdef_mod
60 USE schlieren_mod
61 USE stack_mod
62 USE multi_fvm_mod
63 USE ale_connectivity_mod
64 USE alefvm_mod , only:alefvm_param
65 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas
66 USE matparam_def_mod , ONLY : matparam_struct_
67 USE my_alloc_mod
68 use element_mod , only : nixq,nixc,nixtg
69C-----------------------------------------------
70C I m p l i c i t T y p e s
71C-----------------------------------------------
72#include "implicit_f.inc"
73C-----------------------------------------------
74C C o m m o n B l o c k s
75C-----------------------------------------------
76#include "vect01_c.inc"
77#include "mvsiz_p.inc"
78#include "com01_c.inc"
79#include "com04_c.inc"
80#include "scr25_c.inc"
81#include "param_c.inc"
82#include "task_c.inc"
83#include "spmd_c.inc"
84C-----------------------------------------------
85C D u m m y A r g u m e n t s
86C-----------------------------------------------
87 my_real
88 . func(*),mass(*),x(3,numnod),v(3,numnod),w(3,numnod),thke(*),ehour(*),geo(npropg,numgeo),
89 . anim(*),pm(npropm,nummat),err_thk_sh4(*), err_thk_sh3(*),bufmat(*)
90 INTEGER IPARG(NPARG,NGROUP),IXC(NIXC,NUMELC),IXTG(NIXTG,NUMELTG),EL2FA(*),
91 . IXQ(NIXQ,NUMELQ),IFUNC,NBF,
92 . IADP(*),NBF_L, NBPART,IADG(NSPMD,*),IPM(NPROPMI,NUMMAT),
93 . IGEO(NPROPGI,NUMGEO),INVERT(*), NV46
94 INTEGER, INTENT(IN) :: ITHERM
95 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP), TARGET :: ELBUF_TAB
96 TYPE (STACK_PLY) :: STACK
97 TYPE(BUF_MAT_) ,POINTER :: MBUF
98 TYPE(MULTI_FVM_STRUCT), INTENT(IN) :: MULTI_FVM
99 TYPE (MATPARAM_STRUCT_) ,DIMENSION(NUMMAT) ,INTENT(IN) :: MAT_PARAM
100C-----------------------------------------------
101C L o c a l V a r i a b l e s
102C-----------------------------------------------
103 my_real
104 . evar(mvsiz),dam1(mvsiz),dam2(mvsiz),
105 . wpla(mvsiz),dmax(mvsiz),wpmax(mvsiz),fail(mvsiz),
106 . epst1(mvsiz),epst2(mvsiz),epsf1(mvsiz),epsf2(mvsiz),
107 . sig1(mvsiz),sig2(mvsiz),sig3(mvsiz),
108 . a002(mvsiz),values(mvsiz)
109 my_real
110 . off, p,vonm2,s1,s2,s12,s3,VALUE,value1,value2,dmgmx,fac,
111 . dir1_1,dir1_2,dir2_1,dir2_2,aa,bb,v1,v2,v3,x21,x32,x34,
112 . x41,y21,y32,y34,y41,z21,z32,z34,z41,suma,vr,vs,x31,y31,
113 . z31,e11,e12,e13,e21,e22,e23,sum_,area,x2l,var,rindx,
114 . e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,rx,ry,rz,sx,sy,sz,
115 . vg(5),vly(5),ve(5),dmgmx_ly,evar_tmp,a01,a02,a03,a12,a,
116 . m,ff,gg,nn,
117 . vel(0:3),vfrac(mvsiz,21),phi,err
118 INTEGER I,IDX,I1,II,J,NG,NEL,NPTR,NPTS,NPTT,NLAY,L,IFAIL,ILAY,
119 . IR,IS,IT,IL,MLW, NUVAR,IUS,PTF,PTM,PTS,NFAIL,
120 . N,K,K1,K2,JTURB,MT,IMID,IALEL,IPID,ISH3N,NNI,
121 . NN1,NN2,NN3,NN4,NN5,NN6,NN9,NF,BUF,NVARF,
122 . OFFSET,IHBE,NPTM,NPG, MPT,IPT,IADD,IADR,IPMAT,IFAILT,
123 . IIGEO,IADI,ISUBSTACK,ITHK,NERCVOIS(*),NESDVOIS(*),
124 . LERCVOIS(*),LESDVOIS(*),ID_PLY,NB_PLYOFF,IFRAM_OLD,
125 . JJ(6),NPGT,IADBUF,NUPARAM,IMAT,NS,NRATE,EXPA,IVISC,IU(4),NFRAC
126 INTEGER PID(MVSIZ),MAT(MVSIZ),MATLY(MVSIZ*100),FAILG(MVSIZ),
127 . PTE(4),PTP(4),PTMAT(4),PTVAR(4),LENCOM,NPT_ALL,IPLY,ITRIMAT,IPOS,
128 . ISUBMAT, ISH_EINT, IS_ALE, IS_EULER,IPG,IPINCH,
129 . IMAT_TILLOTSON, NTILLOTSON,NVAREOS,IEOS,IDRAPE,IVAR
130 REAL,DIMENSION(:),ALLOCATABLE:: WAL
131 REAL R4
132 TYPE(G_BUFEL_) ,POINTER :: GBUF
133 TYPE(L_BUFEL_) ,POINTER :: LBUF
134 TYPE(BUF_LAY_) ,POINTER :: BUFLY
135 TYPE(BUF_FAIL_) ,POINTER :: FBUF
136 TYPE(BUF_EOS_) ,POINTER :: EBUF
137 TYPE(L_BUFEL_DIR_) ,POINTER :: LBUF_DIR
138
139 my_real,DIMENSION(:), POINTER :: uvar,offl
140 TYPE(L_BUFEL_) ,POINTER :: LBUF1,LBUF2
141 my_real, DIMENSION(:) ,POINTER :: uparam
142 TARGET :: bufmat
143 my_real tmp(3,4)
144 DATA ns/10/
145 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECTIVITY
146 INTEGER MID
147 my_real :: vi(21) !< submaterial volumes at reference densities (max submat : 21)
148 my_real :: v0i(21) !< submaterial volumes at reference densities (max submat : 21)
149 my_real :: v0g !< global volume at reference density (mixture)
150 my_real :: rho0i(21) !< submaterial initial mass densities (max submat : 21)
151 my_real :: rhoi(21) !< submaterial mass densities (max submat : 21)
152 my_real :: rho0g !< global initial mass density (mixture)
153 LOGICAL detected
154 INTEGER :: IDX0,IDX1,IDX2
155C-----------------------------------------------
156 CALL my_alloc(wal,nbf_l)
157 nn1 = 1
158 nn3 = 1
159 nn4 = nn3 + numelq
160 nn5 = nn4 + numelc
161 nn6 = nn5 + numeltg
162 nn9 = nn6
163 ish_eint = 13242 + 4*mx_ply_anim + 2
164 idx0 = 15921 + 4*mx_ply_anim ! VEL-Y 0:10
165 idx1 = idx0+10+1 ! vel-z 0:10
166 idx2 = idx1+10+1 ! VEL 0:10
167
168 !-------------------------------------------------------!
169 ! INITIALIZATION IF SCHLIEREN DEFINED !
170 !-------------------------------------------------------!
171 IF(ifunc==10672)THEN
172 CALL schlieren_buffer_gathering(nercvois ,nesdvois ,lercvois ,lesdvois, iparg, elbuf_tab, multi_fvm,
173 . itherm)
174 endif!(IFUNC==10672)
175
176C-----------------------------------------------
177 DO ng=1,ngroup
178 CALL initbuf(iparg ,ng ,
179 2 mlw ,nel ,nft ,iad ,ity ,
180 3 npt ,jale ,ismstr ,jeul ,jturb ,
181 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
182 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
183 6 irep ,iint ,igtyp ,israt ,isrot ,
184 7 icsen ,isorth ,isorthg ,ifailure,jsms )
185 IF(mlw /= 13) THEN
186 DO offset = 0,nel-1,nvsiz
187 nft =iparg(3,ng) + offset
188 iad =iparg(4,ng)
189 lft=1
190 llt=min(nvsiz,nel-offset)
191 isubstack = iparg(71,ng)
192 ivisc = iparg(61,ng)
193 is_ale=iparg(7,ng)
194 is_euler=iparg(11,ng)
195 idrape = elbuf_tab(ng)%IDRAPE
196!
197 DO i=1,6
198 jj(i) = nel*(i-1)
199 ENDDO
200!
201c
202 DO i=lft,llt
203! EVAR(I) = ZERO
204 sig1(i) = zero
205 sig2(i) = zero
206 sig3(i) = zero
207 a002(i) = zero
208 ENDDO
209C-----------------------------------------------
210C QUAD OR TRIA
211C-----------------------------------------------
212 IF (ity == 2 .OR.(ity == 7.AND.n2d/=0) ) THEN
213 gbuf => elbuf_tab(ng)%GBUF
214 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
215 uvar => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)%VAR
216 jale=(iparg(7,ng)+iparg(11,ng))
217 jturb=iparg(12,ng)*jale
218 nptr = elbuf_tab(ng)%NPTR
219 npts = elbuf_tab(ng)%NPTS
220 nptt = elbuf_tab(ng)%NPTT
221 nlay = elbuf_tab(ng)%NLAY
222 nuvar = elbuf_tab(ng)%BUFLY(1)%NVAR_MAT
223C
224 DO i=lft,llt
225 func(el2fa(nn3+nft+i))= zero
226 ENDDO
227C
228 IF (ifunc == 1)THEN ! EPSP
229 IF (mlw == 10 .OR. mlw == 21) THEN
230 DO i=lft,llt
231 func(el2fa(nn3+nft+i)) = lbuf%EPSQ(i)
232 ENDDO
233 ELSEIF (mlw == 24) THEN ! and others to add
234 DO i=lft,llt
235 func(el2fa(nn3+nft+i)) = lbuf%VK(i)
236 ENDDO
237 ELSEIF (mlw == 6 .OR. mlw == 17 .OR. mlw == 11) THEN ! and others to add
238 DO i=lft,llt
239 func(el2fa(nn3+nft+i)) = lbuf%RK(i)
240 ENDDO
241 ELSEIF (mlw >=28 .AND. mlw /= 49 .and. nuvar > 0) THEN
242 DO i=lft,llt
243 func(el2fa(nn3+nft+i)) = uvar(i)
244 ENDDO
245 ELSE
246 IF (gbuf%G_PLA > 0) THEN
247 DO i=lft,llt
248 func(el2fa(nn3+nft+i)) = gbuf%PLA(i)
249 ENDDO
250 ENDIF
251 ENDIF
252 ELSEIF(ifunc == 2)THEN
253 DO i=lft,llt
254 func(el2fa(nn3+nft+i)) = gbuf%RHO(i)
255 ENDDO
256 ELSEIF(ifunc == 3)THEN
257 DO i=lft,llt
258 n = i + nft
259 ialel=iparg(7,ng)+iparg(11,ng)
260 IF(ialel == 0)THEN
261 mt=ixq(1,n)
262 VALUE = gbuf%EINT(i)/max(em30,pm(1,mt))
263 ELSE
264 VALUE = gbuf%EINT(i)/max(em30,gbuf%RHO(i))
265 ENDIF
266 func(el2fa(nn3+n)) = VALUE
267 ENDDO
268 ELSEIF(ifunc == 4)THEN
269 IF(gbuf%G_TEMP > 0)THEN
270 DO i=lft,llt
271 func(el2fa(nn3+nft+i)) = gbuf%TEMP(i)
272 ENDDO
273 ELSE
274 DO i=lft,llt
275 func(el2fa(nn3+nft+i)) = zero
276 ENDDO
277 ENDIF
278 ELSEIF(ifunc == 6)THEN
279 DO i=lft,llt
280 p = - (gbuf%SIG(jj(1) + i)
281 . + gbuf%SIG(jj(2) + i)
282 . + gbuf%SIG(jj(3) + i))*third
283 func(el2fa(nn3+nft+i)) = p
284 ENDDO
285 ELSEIF(ifunc == 7)THEN
286 DO i=lft,llt
287 p = - (gbuf%SIG(jj(1) + i)
288 . + gbuf%SIG(jj(2) + i)
289 . + gbuf%SIG(jj(3) + i) )*third
290 s1 = gbuf%SIG(jj(1) + i) + p
291 s2 = gbuf%SIG(jj(2) + i) + p
292 s3 = gbuf%SIG(jj(3) + i) + p
293 vonm2 = three*(gbuf%SIG(jj(4) + i)**2
294 . + half*(s1**2+s2**2+s3**2) )
295 func(el2fa(nn3+nft+i)) = sqrt(vonm2)
296 ENDDO
297 ELSEIF(ifunc == 8 .AND. jturb/=0)THEN
298 DO i=lft,llt
299 func(el2fa(nn3+nft+i)) = gbuf%RK(i)
300 ENDDO
301 ELSEIF(ifunc == 9 )THEN
302 IF (mlw == 6 .OR. mlw == 17.AND.jturb/=0) THEN
303 DO i=lft,llt
304 n = i + nft
305 mt=ixq(1,n)
306 func(el2fa(nn3+n))=pm(81,mt)*gbuf%RK(i)**2/
307 . max(em15,gbuf%RE(i))
308 ENDDO
309 ELSEIF(mlw == 46 .OR. mlw == 47)THEN
310 DO i=lft,llt
311 func(el2fa(nn3+nft+i))= uvar(i)
312 ENDDO
313 ENDIF
314 ELSEIF(ifunc == 10 )THEN ! VORTX
315 IF (mlw == 6 .OR. mlw == 17) THEN
316 DO i=lft,llt
317 func(el2fa(nn3+nft+i)) = lbuf%VK(i)
318 ENDDO
319 ELSEIF(mlw == 46 .OR. mlw == 47)THEN
320 DO i=lft,llt
321 func(el2fa(nn3+nft+i)) = uvar(nel+i)
322 ENDDO
323 ENDIF
324 ELSEIF((ifunc == 11.OR.ifunc == 12.OR.ifunc == 13)
325 . .AND.mlw == 24)THEN ! DAM(3)
326 DO i=lft,llt
327 func(el2fa(nn3+nft+i)) = lbuf%DAM(jj(ifunc-10) + i)
328 ENDDO
329 ELSEIF(ifunc == 14)THEN ! Sigx
330 DO i=lft,llt
331 func(el2fa(nn3+nft+i)) = gbuf%SIG(jj(3) + i)
332 ENDDO
333 ELSEIF(ifunc == 15)THEN ! Sigy
334 DO i=lft,llt
335 func(el2fa(nn3+nft+i)) = gbuf%SIG(jj(1) + i)
336 ENDDO
337 ELSEIF(ifunc == 16)THEN ! Sigz
338 DO i=lft,llt
339 func(el2fa(nn3+nft+i)) = gbuf%SIG(jj(2) + i)
340 ENDDO
341 ELSEIF(ifunc == 17.OR.ifunc == 18)THEN ! Sigxy
342 DO i=lft,llt
343 func(el2fa(nn3+nft+i)) = gbuf%SIG(jj(4) + i)
344 ENDDO
345c----
346 ELSEIF(ifunc>=20.AND.ifunc<=24.AND.
347 . (mlw == 28.OR.mlw == 29.OR.mlw == 30.OR.
348 . mlw == 31.OR.mlw == 52.OR.mlw == 79))THEN
349c---- USER VARIABLES from 1 to 5
350
351 ius = ifunc - 20
352 DO i=lft,llt
353 n = i + nft
354 mt = ixq(1,n)
355 nuvar = ipm(8,mt)
356 IF (nuvar > ius) func(el2fa(nn3+n)) = uvar(ius*nel + i)
357 ENDDO
358 ELSEIF(ifunc == 25)THEN
359 DO i=lft,llt
360 n = i + nft
361 func(el2fa(nn3+nft+i)) = ehour(n)
362 ENDDO
363c---
364 ELSEIF (ifunc == 26) THEN
365 IF (gbuf%G_EPSD > 0) THEN
366 DO i = lft,llt
367 func(el2fa(nn3+nft+i)) = gbuf%EPSD(i)
368 ENDDO
369 ENDIF
370c----
371 ELSEIF (ifunc>=27 .AND. ifunc<=39 .AND.
372 . (mlw == 28.OR.mlw == 29.OR.mlw == 30.OR.mlw == 31.OR.
373 . mlw == 52))THEN
374c---- USER VARIABLES from 6 to 60
375 ius = ifunc - 22 ! IUS = (n-1)'th user variable in UVAR
376 DO i=lft,llt
377 n = i + nft
378 mt=ixq(1,n)
379 nuvar = ipm(8,mt)
380 IF (nuvar>ius) func(el2fa(nn3+n)) = uvar(ius*nel + i)
381 ENDDO
382
383 !/ANIM/ELEM/VFRAC law 20
384 ELSEIF(mlw == 20 .AND. (ifunc == 10248.OR.ifunc == 10249))THEN
385 DO i=lft,llt
386 func(el2fa(nn3+nft+i)) =
387 . elbuf_tab(ng)%BUFLY(ifunc-10248+1)%LBUF(1,1,1)%VOL(i)
388 . / elbuf_tab(ng)%GBUF%VOL(i)
389 ENDDO
390
391 !/ANIM/ELEM/VFRAC law 37
392 ELSEIF(mlw == 37 .AND. (ifunc == 10248.OR.ifunc == 10249))THEN
393 ius=3+(ifunc-10248) !law37 user4 and user5
394 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
395 DO i=lft,llt
396 func(el2fa(nn3+nft+i)) = mbuf%VAR(i+ius*nel)
397 ENDDO
398
399 !/ANIM/ELEM/VFRAC law 51
400 ELSEIF(mlw == 51 .AND. (ifunc >= 10248.AND.ifunc <= 10251))THEN
401 imat = ixq(1,nft+1)
402 iadbuf = ipm(7,imat)
403 nuparam= ipm(9,imat)
404 uparam => bufmat(iadbuf:iadbuf+nuparam-1)
405 isubmat = (ifunc-10247)
406 isubmat = uparam(276+isubmat) !bijective order
407 ius=m51_n0phas+(isubmat-1)*m51_nvphas
408 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
409 DO i=lft,llt
410 func(el2fa(nn3+nft+i)) = mbuf%VAR(i+ius*nel)
411 ENDDO
412
413 ELSEIF(mlw == 151 .AND. (ifunc >= 10248.AND.ifunc <= 10250))THEN
414 ius= ifunc - 10248 + 1
415 IF(ius<=nlay)THEN
416 DO i=lft,llt
417 func(el2fa(nn3+nft+i)) = elbuf_tab(ng)%BUFLY(ius)%LBUF(1,1,1)%VOL(i) / gbuf%VOL(i)
418 ENDDO
419 ELSE
420 DO i=lft,llt
421 func(el2fa(nn3+nft+i)) = zero
422 ENDDO
423 ENDIF
424
425 !Burn Fraction
426 ELSEIF(ifunc == 10252)THEN
427 IF(elbuf_tab(ng)%GBUF%G_BFRAC > 0 .AND. n2d > 0)THEN
428 DO i=lft,llt
429 func(el2fa(nn3+nft+i)) = elbuf_tab(ng)%GBUF%BFRAC(i)
430 ENDDO
431 ELSE
432 DO i=lft,llt
433 func(el2fa(nn3+nft+i)) = zero
434 ENDDO
435 ENDIF
436
437 ELSEIF(ifunc == 10671)THEN
438 IF (mlw == 151) THEN
439 DO i = 1, nel
440 func(el2fa(nn3+nft+i)) = multi_fvm%SOUND_SPEED(i + nft)
441 ENDDO
442 ELSE
443 l = elbuf_tab(ng)%BUFLY(1)%L_SSP
444 IF(elbuf_tab(ng)%BUFLY(1)%L_SSP /= 0)THEN
445 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
446 DO i=1,nel
447 func(el2fa(nn3+nft+i)) = lbuf%SSP(i)
448 ENDDO
449 ENDIF
450 ENDIF
451
452 ELSEIF(ifunc == 10672) THEN
453 ialel=iparg(7,ng)+iparg(11,ng)
454 IF(ialel == 0)THEN
455 DO i=lft,llt
456 func(el2fa(nn3+nft+i)) = zero
457 ENDDO
458 ELSE
459 IF(ity == 2)THEN
460 CALL output_schlieren(
461 1 evar , ixq , x ,
462 2 iparg , wa_l , elbuf_tab , ale_connectivity , gbuf%VOL,
463 3 ng , nixq , ity)
464 ELSEIF(ity == 7 .AND. n2d /= 0)THEN
465 CALL output_schlieren(
466 1 evar , ixtg , x ,
467 2 iparg , wa_l , elbuf_tab , ale_connectivity , gbuf%VOL,
468 3 ng , nixtg , ity)
469 ENDIF
470 DO i=lft,llt
471 func(el2fa(nn3+nft+i)) = evar(i)
472 ENDDO
473 ENDIF
474!---
475 ELSEIF (ifunc == 10677) THEN ! /ANIM/ELEM/SIGEQ
476 ! equivalent stress - other than VON MISES
477 IF (gbuf%G_SEQ > 0) THEN
478C------------------
479 ! Total number of integration points
480 npgt = 0
481 DO il=1,nlay
482 bufly => elbuf_tab(ng)%BUFLY(il)
483 npgt = npgt + bufly%NPTT*nptr*npts
484 ENDDO
485 ! Average equivalent stress on integration points
486 DO i=lft,llt
487 evar_tmp = zero
488 DO il=1,nlay
489 bufly => elbuf_tab(ng)%BUFLY(il)
490 DO it=1,bufly%NPTT
491 DO ir=1,nptr
492 DO is=1,npts
493 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
494 evar_tmp = evar_tmp + lbuf%SEQ(i)/npgt
495 ENDDO
496 ENDDO
497 ENDDO
498 ENDDO
499 func(el2fa(nn3+nft+i)) = evar_tmp
500 ENDDO
501C------------------
502 ELSE ! VON MISES
503 DO i=lft,llt
504 p = - (gbuf%SIG(jj(1) + i)
505 . + gbuf%SIG(jj(2) + i)
506 . + gbuf%SIG(jj(3) + i))*third
507 s1 = gbuf%SIG(jj(1) + i) + p
508 s2 = gbuf%SIG(jj(2) + i) + p
509 s3 = gbuf%SIG(jj(3) + i) + p
510 vonm2 = three*(gbuf%SIG(jj(4) + i)**2
511 . + half*(s1**2+s2**2+s3**2))
512 func(el2fa(nn3+nft+i)) = sqrt(vonm2)
513 ENDDO
514 ENDIF ! IF (GBUF%G_SEQ > 0)
515!---
516 !=== ARTIFICIAL VISCOSITY ===!
517 ELSEIF(ifunc == 11888)THEN
518 IF (gbuf%G_QVIS > 0) THEN
519 DO i=lft,llt
520 func(el2fa(nn3+nft+i)) = gbuf%QVIS(i)
521 ENDDO
522 ELSE
523 DO i=lft,llt
524 func(el2fa(nn3+nft+i)) = zero
525 ENDDO
526 ENDIF
527 !=== DETONATION TIME ===!
528 ELSEIF (ifunc == 11889) THEN
529 IF (gbuf%G_TB > 0) THEN
530 DO i=lft,llt
531 func(el2fa(nn3+nft+i)) = -gbuf%TB(i)
532 ENDDO
533 ELSE
534 DO i=lft,llt
535 func(el2fa(nn3+nft+i)) = zero
536 ENDDO
537 ENDIF
538 ! === LAW 20 === !
539 !DENS BIMAT LAW20 - or future 2d law51
540 ELSEIF(ifunc>=11890 .AND. ifunc<=11893)THEN
541 IF (mlw == 51) THEN
542 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
543 DO i = lft, llt
544 n = i + nft
545 VALUE = zero
546 itrimat = ifunc - 11890 + 1
547 ipos = 12
548 k = llt * (m51_n0phas + (itrimat - 1) * m51_nvphas + ipos - 1)
549 VALUE = mbuf%VAR(k + i)
550 func(el2fa(nn3+n)) = VALUE
551 ENDDO
552 ELSE
553 DO i=lft,llt
554 VALUE = zero
555 n = i + nft
556 ialel = iparg(7,ng)+iparg(11,ng)
557 IF(ialel /= 0 .AND. mlw == 20)THEN
558 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
559 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
560 value1 = lbuf1%RHO(i)
561 value2 = lbuf2%RHO(i)
562 VALUE = zero
563 IF(ifunc == 11890)VALUE=value1
564 IF(ifunc == 11891)VALUE=value2
565 ENDIF
566 func(el2fa(nn3+n)) = VALUE
567 ENDDO
568 ENDIF
569 !ENER BY VOLUME BIMAT LAW20 - or future 2d law51
570 ELSEIF(ifunc>=11894 .AND. ifunc<=11897)THEN
571 IF (mlw == 51) THEN
572 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
573 DO i = lft, llt
574 n = i + nft
575 VALUE = zero
576 itrimat = ifunc - 11894 + 1
577 ipos = 8
578 k = llt * (m51_n0phas + (itrimat - 1) * m51_nvphas + ipos - 1)
579 k2 = llt * ((m51_n0phas + (itrimat-1)*m51_nvphas )+12-1)
580 VALUE = mbuf%VAR(k + i) / mbuf%VAR(k2+i)
581 func(el2fa(nn3+n)) = VALUE
582 ENDDO
583 ELSE
584 DO i=lft,llt
585 VALUE = zero
586 n = i + nft
587 ialel = iparg(7,ng)+iparg(11,ng)
588 IF(ialel /= 0 .AND. mlw == 20)THEN
589 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
590 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
591 value1 = lbuf1%EINT(i)/max(em30,lbuf1%RHO(i))
592 value2 = lbuf2%EINT(i)/max(em30,lbuf2%RHO(i))
593 VALUE = zero
594 IF(ifunc == 11894)VALUE=value1
595 IF(ifunc == 11895)VALUE=value2
596 ENDIF
597 func(el2fa(nn3+n)) = VALUE
598 ENDDO
599 ENDIF
600 !TEMP BIMAT LAW20 - or future 2d law51
601 ELSEIF(ifunc>=11898 .AND. ifunc<=11901)THEN
602 IF (mlw == 51) THEN
603 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
604 DO i = lft, llt
605 n = i + nft
606 VALUE = zero
607 itrimat = ifunc - 11898 + 1
608 ipos = 16
609 k = llt * (m51_n0phas + (itrimat - 1) * m51_nvphas + ipos - 1)
610 VALUE = mbuf%VAR(k + i)
611 func(el2fa(nn3+n)) = VALUE
612 ENDDO
613 ELSE
614 DO i=lft,llt
615 VALUE = zero
616 n = i + nft
617 ialel = iparg(7,ng)+iparg(11,ng)
618 IF(ialel /= 0 .AND. mlw == 20)THEN
619 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
620 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
621 IF(elbuf_tab(ng)%BUFLY(1)%L_TEMP>0)value1 = lbuf1%TEMP(i)
622 IF(elbuf_tab(ng)%BUFLY(2)%L_TEMP>0)value2 = lbuf2%TEMP(i)
623 VALUE = zero
624 IF(ifunc == 11898)VALUE=value1
625 IF(ifunc == 11899)VALUE=value2
626 ENDIF
627 func(el2fa(nn3+n)) = VALUE
628 ENDDO
629 ENDIF
630 !PRES BIMAT LAW20 - or future 2d law51
631 ELSEIF(ifunc>=11902 .AND. ifunc<=11905)THEN
632 IF (mlw == 51) THEN
633 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
634 DO i = lft, llt
635 n = i + nft
636 VALUE = zero
637 itrimat = ifunc - 11902 + 1
638 ipos = 18
639 k = llt * (m51_n0phas + (itrimat - 1) * m51_nvphas + ipos - 1)
640 VALUE = mbuf%VAR(k + i)
641 func(el2fa(nn3+n)) = VALUE
642 ENDDO
643 ELSE
644 DO i=lft,llt
645 VALUE = zero
646 n = i + nft
647 ialel = iparg(7,ng)+iparg(11,ng)
648 IF(ialel /= 0 .AND. mlw == 20)THEN
649 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
650 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
651 value1 = - (lbuf1%SIG(jj(1) + i) +
652 . lbuf1%SIG(jj(2) + i) +
653 . lbuf1%SIG(jj(3) + i))*third
654 value2 = - (lbuf2%SIG(jj(1) + i) +
655 . lbuf2%SIG(jj(2) + i) +
656 . lbuf2%SIG(jj(3) + i))*third
657 VALUE = zero
658 IF(ifunc == 11902)VALUE=value1
659 IF(ifunc == 11903)VALUE=value2
660 ENDIF
661 func(el2fa(nn3+n)) = VALUE
662 ENDDO
663 ENDIF
664 !PLASTIC STRAIN BIMAT LAW20 - or future 2d law51
665 ELSEIF(ifunc>=11906 .AND. ifunc<=11909)THEN
666 DO i=lft,llt
667 VALUE = zero
668 value1 = zero
669 value2 = zero
670 n = i + nft
671 ialel = iparg(7,ng)+iparg(11,ng)
672 IF(ialel /= 0 .AND. mlw == 20)THEN
673 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
674 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
675 IF(elbuf_tab(ng)%BUFLY(1)%L_PLA>0)value1 = lbuf1%PLA(i)
676 IF(elbuf_tab(ng)%BUFLY(2)%L_PLA>0)value2 = lbuf2%PLA(i)
677 VALUE = zero
678 IF(ifunc == 11906)VALUE=value1
679 IF(ifunc == 11907)VALUE=value2
680 ENDIF
681 func(el2fa(nn3+n)) = VALUE
682 ENDDO
683 !SOUND SPEED BIMAT LAW20 - or future 2d law51
684 ELSEIF(ifunc>=11910 .AND. ifunc<=11913)THEN
685 IF (mlw == 51) THEN
686 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
687 DO i = lft, llt
688 n = i + nft
689 VALUE = zero
690 itrimat = ifunc - 11910 + 1
691 ipos = 14
692 k = llt * (m51_n0phas + (itrimat - 1) * m51_nvphas + ipos - 1)
693 VALUE = mbuf%VAR(k + i)
694 func(el2fa(nn3+n)) = VALUE
695 ENDDO
696 ELSE
697 DO i=lft,llt
698 VALUE = zero
699 n = i + nft
700 ialel = iparg(7,ng)+iparg(11,ng)
701 IF(ialel /= 0 .AND. mlw == 20)THEN
702 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
703 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
704 value1 = lbuf1%SSP(i)
705 value2 = lbuf2%SSP(i)
706 VALUE = zero
707 IF(ifunc == 11910)VALUE=value1
708 IF(ifunc == 11911)VALUE=value2
709 ENDIF
710 func(el2fa(nn3+n)) = VALUE
711 ENDDO
712 ENDIF
713!VOLUME SPEED BIMAT LAW20 - or future 2d law51
714 ELSEIF(ifunc>=11914 .AND. ifunc<=11917)THEN
715 IF (mlw == 51) THEN
716 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
717 DO i = lft, llt
718 n = i + nft
719 VALUE = zero
720 itrimat = ifunc - 11914 + 1
721 ipos = 11
722 k = llt * (m51_n0phas + (itrimat - 1) * m51_nvphas + ipos - 1)
723 VALUE = mbuf%VAR(k + i)
724 func(el2fa(nn3+n)) = VALUE
725 ENDDO
726 ELSE
727 DO i=lft,llt
728 VALUE = zero
729 n = i + nft
730 ialel = iparg(7,ng)+iparg(11,ng)
731 IF(ialel /= 0 .AND. mlw == 20)THEN
732 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
733 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
734 value1 = lbuf1%VOL(i)
735 value2 = lbuf2%VOL(i)
736 VALUE = zero
737 IF(ifunc == 11914)VALUE=value1
738 IF(ifunc == 11915)VALUE=value2
739 ENDIF
740 func(el2fa(nn3+n)) = VALUE
741 ENDDO
742 ENDIF
743 !MASS BIMAT LAW20 - or future 2d law51
744 ELSEIF(ifunc>=11918 .AND. ifunc<=11921)THEN
745 IF (mlw == 51) THEN
746 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
747 DO i = lft, llt
748 n = i + nft
749 VALUE = zero
750 itrimat = ifunc - 11918 + 1
751 ipos = 0
752 k = llt * (m51_n0phas + (itrimat - 1) * m51_nvphas + ipos - 1)
753 VALUE = mbuf%VAR(k + i)
754 func(el2fa(nn3+n)) = VALUE
755 ENDDO
756 ELSE
757 DO i=lft,llt
758 VALUE = zero
759 n = i + nft
760 ialel = iparg(7,ng)+iparg(11,ng)
761 IF(ialel /= 0 .AND. mlw == 20)THEN
762 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
763 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
764 value1 = lbuf1%VOL(i) * lbuf1%RHO(i)
765 value2 = lbuf2%VOL(i) * lbuf2%RHO(i)
766 VALUE = zero
767 IF(ifunc == 11918)VALUE=value1
768 IF(ifunc == 11919)VALUE=value2
769 ENDIF
770 func(el2fa(nn3+n)) = VALUE
771 ENDDO
772 ENDIF
773 !ARTIFICIAL VISCOSITY BIMAT LAW20 - or future 2d law51
774 ELSEIF(ifunc>=11922 .AND. ifunc<=11925)THEN
775 DO i=lft,llt
776 VALUE = zero
777 n = i + nft
778 ialel = iparg(7,ng)+iparg(11,ng)
779 IF(ialel /= 0 .AND. mlw == 20)THEN
780 lbuf1 => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
781 lbuf2 => elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)
782 value1 = lbuf1%QVIS(i)
783 value2 = lbuf2%QVIS(i)
784 VALUE = zero
785 IF(ifunc == 11922)VALUE=value1
786 IF(ifunc == 11923)VALUE=value2
787 ENDIF
788 func(el2fa(nn3+n)) = VALUE
789 ENDDO
790 ELSEIF(ifunc == 13242 + 4*mx_ply_anim )THEN
791 IF(gbuf%G_DT>0)THEN
792 DO i=lft,llt
793 func(el2fa(nn3+nft+i)) = gbuf%DT(i)
794 ENDDO
795 ENDIF
796c------------------------------------ Mach Number
797 ELSEIF(ifunc == 13547 + 4*mx_ply_anim + 1000 + 2)THEN
798 IF (mlw == 151) THEN
799 DO i = 1, nel
800 vel(1) = multi_fvm%VEL(1, i + nft)
801 vel(2) = multi_fvm%VEL(2, i + nft)
802 vel(3) = multi_fvm%VEL(3, i + nft)
803 vel(0) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))
804 func(el2fa(nn3+nft+i)) = vel(0)/multi_fvm%SOUND_SPEED(i + nft)
805 ENDDO
806 ELSEIF(alefvm_param%ISOLVER>1)THEN
807 l = elbuf_tab(ng)%BUFLY(1)%L_SSP
808 IF(elbuf_tab(ng)%BUFLY(1)%L_SSP /= 0)THEN
809 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
810 DO i=1,nel
811 vel(1) = gbuf%MOM(jj(1) + i) / gbuf%RHO(i)
812 vel(2) = gbuf%MOM(jj(2) + i) / gbuf%RHO(i)
813 vel(3) = gbuf%MOM(jj(3) + i) / gbuf%RHO(i)
814 vel(0) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))
815 func(el2fa(nn3+nft+i)) = vel(0)/lbuf%SSP(i)
816 ENDDO
817 ENDIF
818 ELSE
819 l = elbuf_tab(ng)%BUFLY(1)%L_SSP
820 IF(elbuf_tab(ng)%BUFLY(1)%L_SSP /= 0)THEN
821 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
822 IF(is_ale /= 0)THEN
823 !ale
824 DO i=1,nel
825 tmp(1,1:4)=v(1,ixq(2:5,i+nft))-w(1,ixq(2:5,i+nft))
826 tmp(2,1:4)=v(2,ixq(2:5,i+nft))-w(2,ixq(2:5,i+nft))
827 tmp(3,1:4)=v(3,ixq(2:5,i+nft))-w(3,ixq(2:5,i+nft))
828 vel(1) = sum(tmp(1,1:4))*fourth
829 vel(2) = sum(tmp(2,1:4))*fourth
830 vel(3) = sum(tmp(3,1:4))*fourth
831 func(el2fa(nn3+nft+i)) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))/lbuf%SSP(i)
832 ENDDO
833 ELSE
834 !euler and lagrange
835 DO i=1,nel
836 tmp(1,1:4)=v(1,ixq(2:5,i+nft))
837 tmp(2,1:4)=v(2,ixq(2:5,i+nft))
838 tmp(3,1:4)=v(3,ixq(2:5,i+nft))
839 vel(1) = sum(tmp(1,1:4))*fourth
840 vel(2) = sum(tmp(2,1:4))*fourth
841 vel(3) = sum(tmp(3,1:4))*fourth
842 func(el2fa(nn3+nft+i)) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))/lbuf%SSP(i)
843 ENDDO
844 ENDIF
845 ENDIF
846 ENDIF
847c------------------------------------ Color Function
848 ELSEIF(ifunc == 13547 + 4*mx_ply_anim + 1000 + 3)THEN
849 gbuf => elbuf_tab(ng)%GBUF
850 IF (mlw == 151) THEN
851 nfrac=nlay
852 DO imat=1,nlay
853 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1,1,1)
854 DO i=1,nel
855 vfrac(i,imat) = lbuf%VOL(i) / gbuf%VOL(i)
856 ENDDO
857 ENDDO
858 ELSEIF(mlw == 20)THEN
859 nfrac=2
860 DO i=1,nel
861 vfrac(i,1) = elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)%VOL(i) / gbuf%VOL(i)
862 vfrac(i,2) = elbuf_tab(ng)%BUFLY(2)%LBUF(1,1,1)%VOL(i) / gbuf%VOL(i)
863 ENDDO
864 ELSEIF(mlw == 37)THEN
865 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
866 nfrac=2
867 DO i=1,nel
868 vfrac(i,1) = mbuf%VAR(i+3*nel)
869 vfrac(i,2) = mbuf%VAR(i+4*nel)
870 ENDDO
871 ELSEIF(mlw == 51)THEN
872 !get UPARAM
873 imat = ixq(1,nft+1)
874 iadbuf = ipm(7,imat)
875 nuparam= ipm(9,imat)
876 uparam => bufmat(iadbuf:iadbuf+nuparam-1)
877 !bijective order !indexes
878 isubmat = uparam(276+1); iu(1)=m51_n0phas+(isubmat-1)*m51_nvphas
879 isubmat = uparam(276+2); iu(2)=m51_n0phas+(isubmat-1)*m51_nvphas
880 isubmat = uparam(276+3); iu(3)=m51_n0phas+(isubmat-1)*m51_nvphas
881 isubmat = uparam(276+4); iu(4)=m51_n0phas+(isubmat-1)*m51_nvphas
882 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
883 nfrac=4
884 DO i=1,nel
885 vfrac(i,1) = mbuf%VAR(i+iu(1)*nel)
886 vfrac(i,2) = mbuf%VAR(i+iu(2)*nel)
887 vfrac(i,3) = mbuf%VAR(i+iu(3)*nel)
888 vfrac(i,4) = mbuf%VAR(i+iu(4)*nel)
889 ENDDO
890 ELSE
891 nfrac=0
892 vfrac(1:nel,1:21)=zero
893 ENDIF
894 IF(nfrac>0)THEN
895 DO i=1,nel
896 values(i)=zero
897 DO imat=1,nfrac
898 values(i) = values(i) + vfrac(i,imat)*imat
899 ENDDO
900 func(el2fa(nn3+nft+i))=values(i)
901 ENDDO
902 ELSE
903 DO i=1,nel
904 func(el2fa(nn3+nft+i))=zero
905 ENDDO
906 ENDIF
907 ELSEIF(ifunc == 4*mx_ply_anim + 14566)THEN
908 IF(n2d == 1)THEN
909 fac = two*3.141592653589793238 !axi
910 ELSE
911 fac = one !plane stress
912 ENDIF
913 DO i = lft, llt
914 n = i + nft
915 func(el2fa(nn3+n)) = fac*gbuf%VOL(i)
916 ENDDO
917 !ELSEIF(IFUNC == 4*MX_PLY_ANIM + 14567)THEN
918 ! ...
919 !ELSEIF(IFUNC == 4*MX_PLY_ANIM + 14568)THEN
920 ! ...
921C-------------------------------------
922 ELSE IF (ifunc == 10676) THEN
923C-------------------------------------
924 !SPMD DOMAIN
925 DO i=lft,llt
926 func(el2fa(nn3+nft+i)) = ispmd
927 ENDDO
928C-------------------------------------
929 ELSEIF (ifunc == 14595 + 4*mx_ply_anim .AND. (gbuf%G_TSAIWU > 0)) THEN
930C-------------------------------------
931 ! /ANIM/ELEM/TSAIWU
932 bufly => elbuf_tab(ng)%BUFLY(1)
933 DO i=lft,llt
934 DO ir=1,nptr
935 DO is=1,npts
936 DO it=1,nptt
937 func(el2fa(nn3+nft+i)) =
938 . func(el2fa(nn3+nft+i))
939 . + bufly%LBUF(ir,is,it)%TSAIWU(i)/(nptt*nptr*npts)
940 ENDDO
941 ENDDO
942 ENDDO
943 ENDDO
944c------------------------------------ Tillotson Region id
945 ! /ANIM/ELEM/TILLOTSON
946 ELSEIF( ifunc == 15898 + 4*mx_ply_anim ) THEN
947 DO i=lft,llt
948 func(el2fa(nn3+nft+i)) = zero
949 ENDDO
950 mt = ixq(1,nft+1)
951 IF (mlw == 151) THEN
952 nlay = elbuf_tab(ng)%NLAY
953 !count number of submaterial based on /EOS/TILLOTSON (IEOS=3)
954 ntillotson = 0
955 DO imat=1,nlay
956 ieos = ipm(4, mat_param(mt)%MULTIMAT%MID(imat) )
957 IF(ieos == 3)THEN
958 ntillotson = ntillotson + 1
959 imat_tillotson = imat
960 ENDIF
961 ENDDO
962 !several Tillotson EoS Value= sum ( Region_i*10**(i-1), i=1,imat)
963 IF(ntillotson > 1)THEN
964 fac=one
965 DO imat=1,nlay
966 ieos = ipm(4, mat_param(mt)%MULTIMAT%MID(imat) )
967 IF(ieos == 3)THEN
968 ebuf => elbuf_tab(ng)%BUFLY(imat)%EOS(1,1,1)
969 nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
970 DO i=1,nel
971 func(el2fa(nn3+nft+i)) = func(el2fa(nn3+nft+i)) + ebuf%VAR(i) * fac
972 ENDDO
973 ENDIF
974 fac=fac*ten
975 ENDDO
976 !single Tillotson EoS Value= Region_i
977 ELSEIF(ntillotson == 1)THEN
978 ebuf => elbuf_tab(ng)%BUFLY(imat_tillotson)%EOS(1,1,1)
979 nvareos = elbuf_tab(ng)%BUFLY(imat_tillotson)%NVAR_EOS
980 DO i=1,nel
981 func(el2fa(nn3+nft+i)) = ebuf%VAR(i)
982 ENDDO
983 ENDIF
984 ELSE
985 !monomaterial law
986 ieos = ipm(4,mt)
987 IF(ieos == 3)THEN
988 ebuf => elbuf_tab(ng)%BUFLY(1)%EOS(1,1,1)
989 nvareos = elbuf_tab(ng)%BUFLY(1)%NVAR_EOS
990 DO i=1,nel
991 func(el2fa(nn3+nft+i)) = ebuf%VAR(i)
992 ENDDO
993 ENDIF
994 ENDIF
995
996c------------------------------------ Volumetric Strain (GLOBAL VSTRAIN)
997 elseif(ifunc == 15899 + 4*mx_ply_anim .and. n2d > 0) then
998!--------------------------------------------------
999 DO i=lft,llt
1000 func(el2fa(nn3+nft+i)) = zero
1001 ENDDO
1002
1003 if(ity == 2)then
1004 mt = ixq(1,nft+1)
1005 elseif(ity == 7 .and. n2d > 0)then
1006 mt = ixtg(1,nft+1)
1007 endif
1008 imat = mt
1009
1010 do i=1,nel
1011
1012 if(mlw == 151)then
1013 !multimaterial 151 (collocated scheme)
1014 do ilay=1,multi_fvm%nbmat
1015 mid = mat_param(mt)%multimat%mid(ilay)
1016 rho0i(ilay) = pm(89,mid)
1017 vi(ilay) = multi_fvm%phase_alpha(ilay,i+nft) * gbuf%vol(i)
1018 v0i(ilay) = multi_fvm%phase_rho(ilay,i+nft) * vi(ilay) / rho0i(ilay) !rho0.V0 = rho.V
1019 enddo
1020 v0g = sum(v0i)
1021 rho0g = zero
1022 do ilay=1,multi_fvm%nbmat
1023 rho0g = rho0g + rho0i(ilay)*v0i(ilay)
1024 end do
1025 rho0g = rho0g / v0g
1026 func(el2fa(nn3+nft+i)) = multi_fvm%rho(i+nft) / rho0g - one
1027
1028 elseif(mlw == 51)then
1029 !multimaterial 51 (staggered scheme)
1030 iadbuf = ipm(7,imat)
1031 nuparam= ipm(9,imat)
1032 uparam => bufmat(iadbuf:iadbuf+nuparam-1)
1033 mbuf => elbuf_tab(ng)%bufly(1)%mat(1,1,1)
1034 ipos = 1
1035 !bijective order !indexes
1036 isubmat = nint(uparam(276+1)); iu(1)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1037 isubmat = nint(uparam(276+2)); iu(2)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1038 isubmat = nint(uparam(276+3)); iu(3)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1039 isubmat = nint(uparam(276+4)); iu(4)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1040 vfrac(i,1) = mbuf%var(i+iu(1)*nel)
1041 vfrac(i,2) = mbuf%var(i+iu(2)*nel)
1042 vfrac(i,3) = mbuf%var(i+iu(3)*nel)
1043 vfrac(i,4) = mbuf%var(i+iu(4)*nel)
1044 ipos = 12
1045 !bijective order !indexes
1046 isubmat = nint(uparam(276+1)); iu(1)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1047 isubmat = nint(uparam(276+2)); iu(2)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1048 isubmat = nint(uparam(276+3)); iu(3)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1049 isubmat = nint(uparam(276+4)); iu(4)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1050 rhoi(1) = mbuf%var(i+iu(1)*nel)
1051 rhoi(2) = mbuf%var(i+iu(2)*nel)
1052 rhoi(3) = mbuf%var(i+iu(3)*nel)
1053 rhoi(4) = mbuf%var(i+iu(4)*nel)
1054 do ilay=1,4
1055 mid = mat_param(mt)%multimat%mid(ilay)
1056 rho0i(ilay) = pm(89,mid)
1057 vi(ilay) = vfrac(i,ilay) * gbuf%vol(i)
1058 ipos = 12
1059 v0i(ilay) = rhoi(ilay) * vi(ilay) / rho0i(ilay) !rho0.V0 = rho.V
1060 enddo
1061 v0g = sum(v0i)
1062 rho0g = zero
1063 do ilay=1,4
1064 rho0g = rho0g + rho0i(ilay)*v0i(ilay)
1065 end do
1066 rho0g = rho0g / v0g
1067 func(el2fa(nn3+nft+i))= gbuf%rho(i) / rho0g - one
1068
1069 elseif(mlw == 37)then
1070 !multimaterial 37 (staggered scheme)
1071 iadbuf = ipm(7,imat)
1072 nuparam= ipm(9,imat)
1073 uparam => bufmat(iadbuf:iadbuf+nuparam-1)
1074 mbuf => elbuf_tab(ng)%bufly(1)%mat(1,1,1)
1075 rho0i(1) = uparam(11)
1076 rho0i(2) = uparam(12)
1077 vi(1) = mbuf%var(i+3*nel) * gbuf%vol(i) !UVAR(I,4) = VFRAC1
1078 vi(2) = mbuf%var(i+4*nel) * gbuf%vol(i) !UVAR(I,5) = VFRAC2
1079 rhoi(1) = mbuf%var(i+2*nel) !UVAR(I,3) = RHO1
1080 rhoi(2) = mbuf%var(i+1*nel) !UVAR(I,2) = RHO2
1081 v0i(1) = rhoi(1) * vi(1) / rho0i(1) !rho0.V0 = rho.V
1082 v0i(2) = rhoi(2) * vi(2) / rho0i(2) !rho0.V0 = rho.V
1083 v0g = sum(v0i)
1084 rho0g = zero
1085 do ilay=1,2
1086 rho0g = rho0g + rho0i(ilay)*v0i(ilay)
1087 end do
1088 rho0g = rho0g / v0g
1089 func(el2fa(nn3+nft+i)) = gbuf%rho(i) / rho0g - one
1090
1091 elseif(mlw == 20)then
1092 !multimaterial 20 (staggered scheme)
1093 lbuf1 => elbuf_tab(ng)%bufly(1)%lbuf(1,1,1)
1094 lbuf2 => elbuf_tab(ng)%bufly(2)%lbuf(1,1,1)
1095 mid = mat_param(mt)%multimat%mid(1)
1096 rho0i(1) = pm(89,mid)
1097 mid = mat_param(mt)%multimat%mid(2)
1098 rho0i(2) = pm(89,mid)
1099 vi(1) = lbuf1%vol(i)
1100 vi(2) = lbuf2%vol(i)
1101 rhoi(1) = lbuf1%rho(i)
1102 rhoi(2) = lbuf2%rho(i)
1103 v0i(1) = rhoi(1) * vi(1) / rho0i(1) !rho0.V0 = rho.V
1104 v0i(2) = rhoi(2) * vi(2) / rho0i(2) !rho0.V0 = rho.V
1105 v0g = sum(v0i)
1106 rho0g = zero
1107 do ilay=1,2
1108 rho0g = rho0g + rho0i(ilay)*v0i(ilay)
1109 end do
1110 rho0g = rho0g / v0g
1111 func(el2fa(nn3+nft+i)) = gbuf%rho(i) / rho0g - one
1112
1113 else
1114 !general case (monomaterial law)
1115 if(pm(89,mt) > zero)then
1116 func(el2fa(nn3+nft+i))= gbuf%rho(i) / pm(89,mt) - one
1117 end if
1118 end if
1119
1120 enddo
1121c------------------------------------ Volumetric Strain (submaterial VSTRAIN)
1122 elseif( ifunc >= 15899 + 4*mx_ply_anim +1
1123 . .and. ifunc <= 15899 + 4*mx_ply_anim +10
1124 . .and. n2d > 0) then
1125!--------------------------------------------------
1126 detected = .false.
1127 ilay = ifunc - (15899 + 4*mx_ply_anim)
1128 if(mlw == 151 .and. ilay <= min(10,multi_fvm%nbmat))detected = .true.
1129 if(mlw == 51 .and. ilay <= 4 )detected = .true.
1130 if(mlw == 37 .and. ilay <= 2 )detected = .true.
1131 if(mlw == 20 .and. ilay <= 2 )detected = .true.
1132
1133 if(detected)then
1134
1135 if(ity == 2)then
1136 mt = ixq(1,nft+1)
1137 elseif(ity == 7 .and. n2d > 0)then
1138 mt = ixtg(1,nft+1)
1139 endif
1140 imat = mt
1141
1142 do i=1,nel
1143
1144 if(mlw == 151)then
1145 !multimaterial 151 (collocated scheme)
1146 mid = mat_param(mt)%multimat%mid(ilay)
1147 rho0i(ilay) = pm(89,mid)
1148 vi(ilay) = multi_fvm%phase_alpha(ilay,i+nft) * gbuf%vol(i)
1149 v0i(ilay) = multi_fvm%phase_rho(ilay,i+nft) * vi(ilay) / rho0i(ilay) !rho0.V0 = rho.V
1150 func(el2fa(nn3+nft+i)) = multi_fvm%phase_rho(ilay,i+nft) / rho0i(ilay) - one
1151
1152 elseif(mlw == 51)then
1153 !multimaterial 51 (staggered scheme)
1154 iadbuf = ipm(7,imat)
1155 nuparam= ipm(9,imat)
1156 uparam => bufmat(iadbuf:iadbuf+nuparam-1)
1157 mbuf => elbuf_tab(ng)%bufly(1)%mat(1,1,1)
1158 mid = mat_param(mt)%multimat%mid(ilay)
1159 rho0i(ilay) = pm(89,mid)
1160 ipos = 1
1161 !bijective order !indexes
1162 isubmat = nint(uparam(276+ilay)); iu(1)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1163 vfrac(i,ilay) = mbuf%var(i+iu(ilay)*nel)
1164 vi(ilay) = vfrac(i,ilay) * gbuf%vol(i)
1165 ipos = 12
1166 !bijective order !indexes
1167 isubmat = nint(uparam(276+ilay)); iu(ilay)=m51_n0phas+(isubmat-1)*m51_nvphas + ipos-1
1168 rhoi(ilay) = mbuf%var(i+iu(ilay)*nel)
1169 v0i(ilay) = rhoi(ilay) * vi(ilay) / rho0i(ilay) !rho0.V0 = rho.V
1170 func(el2fa(nn3+nft+i)) = rhoi(ilay) / rho0i(ilay) - one
1171
1172 elseif(mlw == 37)then
1173 !multimaterial 37 (staggered scheme)
1174 iadbuf = ipm(7,imat)
1175 nuparam= ipm(9,imat)
1176 uparam => bufmat(iadbuf:iadbuf+nuparam-1)
1177 mbuf => elbuf_tab(ng)%bufly(1)%mat(1,1,1)
1178 rho0i(ilay) = uparam(10+ilay)
1179 vi(ilay) = mbuf%var(i+(ilay+2)*nel) * gbuf%vol(i)
1180 rhoi(ilay) = mbuf%var(i+(3-ilay)*nel) !UVAR(I,3) = RHO1
1181 v0i(ilay) = rhoi(ilay) * vi(ilay) / rho0i(ilay)
1182 func(el2fa(nn3+nft+i)) = rhoi(ilay) / rho0i(ilay) - one
1183
1184 elseif(mlw == 20)then
1185 !multimaterial 20 (staggered scheme)
1186 lbuf => elbuf_tab(ng)%bufly(ilay)%lbuf(1,1,1)
1187 mid = mat_param(mt)%multimat%mid(ilay)
1188 rho0i(ilay) = pm(89,mid)
1189 vi(ilay) = lbuf%vol(i)
1190 rhoi(ilay) = lbuf%rho(i)
1191 v0i(ilay) = rhoi(ilay) * vi(ilay) / rho0i(ilay) !rho0.V0 = rho.V
1192 func(el2fa(nn3+nft+i)) = rhoi(ilay) / rho0i(ilay) - one
1193
1194 else
1195 !general case (monomaterial law)
1196 func(el2fa(nn3+nft+i)) = zero
1197 end if
1198 enddo
1199 end if
1200
1201
1202c------------------------------------ 2D velocity VEL-GLOBAL
1203 elseif( ifunc >= idx0 .AND. ifunc <= idx0+10) then
1204 ilay = ifunc - idx0
1205 IF(mlw == 151 .AND. ilay == 0)THEN
1206 DO i=lft,llt
1207 vel(1) = multi_fvm%VEL(1, i + nft)
1208 vel(2) = multi_fvm%VEL(2, i + nft)
1209 vel(3) = multi_fvm%VEL(3, i + nft)
1210 func(el2fa(nn3+nft+i)) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))
1211 ENDDO
1212 ELSE
1213 DO i=lft,llt
1214 func(el2fa(nn3+nft+i)) = zero
1215 ENDDO
1216 ENDIF
1217c------------------------------------ 2D velocity VEL-Y
1218 elseif( ifunc >= idx1 .AND. ifunc <= idx1+10) then
1219 ilay = ifunc - idx1
1220 IF(mlw == 151 .AND. ilay == 0)THEN
1221 DO i=lft,llt
1222 vel(2) = multi_fvm%VEL(2, i + nft)
1223 func(el2fa(nn3+nft+i)) = vel(2)
1224 ENDDO
1225 ELSE
1226 DO i=lft,llt
1227 func(el2fa(nn3+nft+i)) = zero
1228 ENDDO
1229 ENDIF
1230c------------------------------------ 2D velocity VEL-Z
1231 elseif( ifunc >= idx2 .AND. ifunc <= idx2+10) then
1232 ilay = ifunc - idx2
1233 IF(mlw == 151 .AND. ilay == 0)THEN
1234 DO i=lft,llt
1235 vel(3) = multi_fvm%VEL(3, i + nft)
1236 func(el2fa(nn3+nft+i)) = vel(3)
1237 ENDDO
1238 ELSE
1239 DO i=lft,llt
1240 func(el2fa(nn3+nft+i)) = zero
1241 ENDDO
1242 ENDIF
1243c------------------------------------
1244 ELSE
1245 DO i=lft,llt
1246 func(el2fa(nn3+nft+i)) = zero
1247 ENDDO
1248 ENDIF
1249C-----------------------------------------------
1250C COQUES 3 N 4 N
1251C-----------------------------------------------
1252 ELSEIF (ity == 3.OR.(ity == 7.AND.n2d==0)) THEN
1253 IF(ity == 3)THEN
1254 nni = nn4
1255 ELSE
1256 nni = nn5
1257 ENDIF
1258 gbuf => elbuf_tab(ng)%GBUF
1259 npt = iparg(6,ng)
1260 ihbe = iparg(23,ng)
1261 irep = iparg(35,ng)
1262 igtyp = iparg(38,ng)
1263 ithk = iparg(28,ng)
1264 mpt = iabs(npt)
1265 nptr = elbuf_tab(ng)%NPTR
1266 npts = elbuf_tab(ng)%NPTS
1267 nptt = elbuf_tab(ng)%NPTT
1268 nlay = elbuf_tab(ng)%NLAY
1269 npg = nptr*npts
1270 nuvar = 0
1271 ipinch= iparg(90,ng)
1272 IF (ihbe==3.AND.ish3nfram==0) THEN
1273 ifram_old =0
1274 ELSE
1275 ifram_old =1
1276 END IF
1277C
1278 IF (igtyp == 51 .OR. igtyp == 52) THEN
1279 npt_all = 0
1280 DO ipt=1,nlay
1281 npt_all = npt_all + elbuf_tab(ng)%BUFLY(ipt)%NPTT
1282 ENDDO
1283 IF (nlay == 1) mpt = max(1,npt_all)
1284 ENDIF
1285C---------------------
1286 DO i=lft,llt
1287 evar(i) = zero ! Default = zero in all cases !
1288 ENDDO
1289C---------------------
1290C
1291 IF (mlw == 0 .OR. mlw == 13) THEN
1292c---
1293 ELSEIF (ifunc == 1 .AND. (mlw /= 15 .AND. mlw /= 25)) THEN ! plastic strain (mid layer : npt/2 + 1)
1294! law15 and law25 it's the plastic work instead plastic strain. it can be outp using /ANIM/ELEM/WPLA ( ifunc=4*MX_PLY_ANIM +13245)
1295 IF (gbuf%G_PLA > 0) THEN
1296 ilay = 1
1297 IF (nlay > 1) ilay = iabs(nlay)/2 + 1
1298 bufly => elbuf_tab(ng)%BUFLY(ilay)
1299 IF (bufly%L_PLA > 0) THEN
1300 IF (npg > 1) THEN
1301 IF(ity == 3) THEN
1302 IF(igtyp == 51 .OR. igtyp == 52) THEN
1303 nptt = bufly%NPTT
1304 DO is = 1,npts
1305 DO ir = 1,nptr
1306 DO it = 1, nptt
1307 DO i=1,nel
1308 evar(i) = evar(i) + fourth*bufly%LBUF(ir,is,it)%PLA(i)/nptt
1309 ENDDO
1310 ENDDO
1311 ENDDO
1312 ENDDO
1313 ELSE
1314 DO i=1,nel
1315 evar(i) = fourth*(bufly%LBUF(1,1,1)%PLA(i) + bufly%LBUF(2,1,1)%PLA(i) +
1316 . bufly%LBUF(1,2,1)%PLA(i) + bufly%LBUF(2,2,1)%PLA(i))
1317 ENDDO
1318 ENDIF ! igtyp
1319 ELSE ! ITY == 7
1320 IF(igtyp == 51 .OR. igtyp == 52) THEN
1321 nptt = bufly%NPTT
1322 DO it = 1,nptt
1323 DO ir =1,npg
1324 DO i=1,nel
1325 evar(i) = evar(i) + third*bufly%LBUF(ir,1,it)%PLA(i)/nptt
1326 ENDDO
1327 ENDDO
1328 ENDDO
1329 ELSE
1330 DO i=1,nel
1331 evar(i) = third*(bufly%LBUF(1,1,1)%PLA(i) + bufly%LBUF(1,1,1)%PLA(i) +
1332 . bufly%LBUF(1,1,1)%PLA(i))
1333 ENDDO
1334 ENDIF ! igtyp
1335 ENDIF ! ity
1336 ELSE ! NPG == 1
1337 IF(igtyp == 51 .OR. igtyp == 52) THEN
1338 nptt = bufly%NPTT ! needed for PID51 (new shell prop)
1339 DO it = 1,nptt
1340 DO i=1,nel
1341 evar(i) = evar(i) + abs(bufly%LBUF(1,1,it)%PLA(i))/nptt
1342 ENDDO
1343 ENDDO
1344 ELSE
1345 nptt = bufly%NPTT !
1346 ipt = iabs(nptt)/2 + 1
1347 DO i=1,nel
1348 evar(i) = abs(bufly%LBUF(1,1,ipt)%PLA(i))
1349 ENDDO
1350 ENDIF ! igtyp
1351 ENDIF ! npg
1352 ENDIF ! BUFLY%L_PLA
1353 ENDIF ! GBUF
1354
1355c---
1356 ELSEIF (ifunc == 2) THEN
1357 IF (mlw == 151) THEN
1358 DO i=lft,llt
1359 evar(i) = gbuf%RHO(i)
1360 ENDDO
1361 ELSE
1362 IF (ity == 3) THEN
1363 DO i=lft,llt
1364 evar(i) = pm(1,ixc(1,nft+i))
1365 ENDDO
1366 ELSEIF (ity == 7) THEN
1367 DO i=lft,llt
1368 evar(i) = pm(1,ixtg(1,nft+i))
1369 ENDDO
1370 ENDIF
1371 ENDIF
1372c---
1373 ELSEIF (ifunc == 3 .AND. mlw == 151) THEN
1374 DO i=lft,llt
1375 evar(i) = gbuf%EINT(i) / gbuf%RHO(i)
1376 ENDDO
1377c---
1378 ELSEIF (ifunc == 3 .OR. ifunc == ish_eint) THEN ! EINT
1379 DO i=lft,llt
1380cc K1 = 2 * I
1381cc K2 = K1 - 1
1382!! K1 = 2*(I-1)+1
1383!! K2 = 2*(I-1)+2
1384 evar(i) = gbuf%EINT(i) + gbuf%EINT(i+llt)
1385 ENDDO
1386c---
1387 ELSEIF (ifunc == 4) THEN ! ELEMENT TEMPERATURE
1388 IF (jthe /= 0) THEN
1389 evar(1:nel) = gbuf%TEMP(1:nel)
1390 ELSE
1391 evar(1:nel) = zero
1392 nptt = 0
1393 DO il=1,nlay
1394 IF (elbuf_tab(ng)%BUFLY(il)%L_TEMP > 0) THEN
1395 nptt = nptt + elbuf_tab(ng)%BUFLY(il)%NPTT
1396 ENDIF
1397 END DO
1398 npg = nptr*npts*nptt
1399 DO il=1,nlay
1400 IF (elbuf_tab(ng)%BUFLY(il)%L_TEMP > 0) THEN
1401 DO it=1,elbuf_tab(ng)%BUFLY(il)%NPTT
1402 DO is=1,npts
1403 DO ir=1,nptr
1404 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
1405 evar(1:nel) = evar(1:nel) + lbuf%TEMP(1:nel)/npg
1406 ENDDO
1407 ENDDO
1408 ENDDO
1409 ENDIF
1410 ENDDO
1411 END IF
1412c---
1413 ELSEIF (ifunc == 5) THEN ! THK
1414 IF (ithk >0) THEN
1415 DO i=lft,llt
1416 evar(i) = gbuf%THK(i)
1417 ENDDO
1418 ELSE
1419 IF (ity == 3) THEN ! SHELL
1420 DO i=lft,llt
1421 evar(i) = thke(nft+i)
1422 ENDDO
1423 ELSEIF (ity == 7) THEN ! SH3N
1424 DO i=lft,llt
1425 evar(i) = thke(nft+i+numelc)
1426 ENDDO
1427 ENDIF ! IF (ITY == 3)
1428 END IF
1429c---
1430 ELSEIF (ifunc == 6 .AND. mlw == 151) THEN
1431 DO i=lft,llt
1432 evar(i) = - third * (gbuf%SIG(i) + gbuf%SIG(i + nel) + gbuf%SIG(i + 2 * nel))
1433 ENDDO
1434c---
1435 ELSEIF (ifunc == 7) THEN ! Von Mises
1436 DO i=lft,llt
1437 s1 = gbuf%FOR(jj(1)+i)
1438 s2 = gbuf%FOR(jj(2)+i)
1439 s12= gbuf%FOR(jj(3)+i)
1440 vonm2= s1*s1 + s2*s2 - s1*s2 + three*s12*s12
1441 evar(i) = sqrt(vonm2)
1442 ENDDO
1443c---
1444 ELSEIF (ifunc == 11) THEN ! DAM1
1445c IF (MLW == 25 .OR. MLW == 15 )THEN
1446c NB13 = NB12 + 2*NEL*MAX(1,NPT)
1447c NB15 = NB12 + 2*NEL*MAX(1,NPT)
1448c NB12 = NB12 + 2*OFFSET*MAX(1,NPT)
1449c NB13 = NB13 + 4*OFFSET*MAX(1,NPT)
1450c DO J=1,NPT
1451c N = (IPT-1)*NEL
1452c DO I=LFT,LLT
1453c K1 = NB13 + N+I
1454c K2 = NB15 + 2*N+2*I-1
1455c EVAR(I)=EVAR(I)+BUFEL(K1)*BUFEL(K2)
1456c ENDDO
1457c ENDDO
1458c ENDIF
1459c---
1460 ELSEIF(ifunc == 12)THEN ! DAM2
1461c IF(MLW == 25.OR.MLW == 15)THEN
1462c NB13 = NB12 + 2*NEL*MAX(1,NPT)
1463c NB15 = NB12 + 2*NEL*MAX(1,NPT)
1464c NB12 = NB12 + 2*OFFSET*MAX(1,NPT)
1465c NB13 = NB13 + 4*OFFSET*MAX(1,NPT)
1466c DO J=1,NPT
1467c N = (IPT-1)*NEL
1468c DO I=LFT,LLT
1469c K1 = NB13 + N+I
1470c K2 = NB15 + 2*N+2*I
1471c
1472c ENDDO
1473c ENDDO
1474c ENDIF
1475c---
1476 ELSEIF(ifunc == 13)THEN ! DAM3
1477c IF(MLW == 25.OR.MLW == 15)THEN
1478c IADD = NB12 + 6*NEL*MAX(1,NPT)
1479c IADD = IADD + OFFSET
1480c DO I=LFT,LLT
1481c EVAR(I) = BUFEL(IADD+I)
1482c ENDDO
1483c ENDIF
1484c---
1485 ELSEIF (ifunc >= 14 .AND. ifunc <= 15) THEN ! Sigx, Sigy - global
1486 ius = ifunc-13
1487 DO i=lft,llt
1488 evar(i) = gbuf%FOR(jj(ius)+i)
1489 ENDDO
1490c---
1491 ELSEIF (ifunc == 16 .AND. ihbe == 11 .AND. ipinch == 1) THEN ! Sigz for pinching shell
1492 DO i=lft,llt
1493 evar(i) = zero
1494 DO ipg=1,4
1495 evar(i) = evar(i) + fourth*gbuf%FORPGPINCH(nel*(ipg-1)+i)
1496 ENDDO
1497 ENDDO
1498c---
1499 ELSEIF (ifunc >= 17 .AND. ifunc <= 19) THEN ! Sigyx, Sigyz, Sigzx - global
1500 ius = ifunc-14
1501 DO i=lft,llt
1502 evar(i) = gbuf%FOR(jj(ius)+i)
1503 ENDDO
1504c---
1505 ELSEIF (ifunc == 26) THEN
1506 evar(lft:llt) = gbuf%EPSD(lft:llt)
1507c---
1508 ELSEIF(ifunc == 2155)THEN
1509 DO i=lft,llt
1510 evar(i) = hundred *(gbuf%THK_I(i)-gbuf%THK(i))/gbuf%THK_I(i)
1511 ENDDO
1512c------------------------------------
1513 ELSEIF (ifunc>=20 .AND. ifunc<=24) THEN
1514C USER VARIABLES from 1 to 5
1515c------------------------------------
1516 IF (mlw==29 .OR. mlw==30 .OR. mlw==31 .OR. mlw>=33) THEN
1517c
1518 ius = ifunc - 20
1519 i1 = ius*nel
1520 IF (nlay > 1) THEN
1521 il = iabs(nlay)/2 + 1
1522 ipt = 1
1523 ELSE
1524 il = 1
1525 ipt = iabs(npt)/2 + 1
1526 ENDIF
1527 nuvar = elbuf_tab(ng)%BUFLY(il)%NVAR_MAT
1528c
1529 IF (mlw == 58 .or. mlw == 158) THEN
1530 DO i=lft,llt
1531 DO ir = 1, nptr
1532 DO is = 1, npts
1533 uvar=>elbuf_tab(ng)%BUFLY(il)%MAT(ir,is,ipt)%VAR
1534 IF (ius==4 .OR. ius==5) THEN ! convert true to eng strain
1535 evar(i) = evar(i) + exp(uvar(i1+i) - one) / npg
1536 ELSE
1537 evar(i) = evar(i) + uvar(i1 + i) / npg
1538 ENDIF
1539 ENDDO
1540 ENDDO
1541 ENDDO
1542 ELSE
1543 DO i=lft,llt
1544 IF (nuvar > ius) THEN
1545 DO ir = 1, nptr
1546 DO is = 1, npts
1547 uvar=>elbuf_tab(ng)%BUFLY(il)%MAT(ir,is,ipt)%VAR
1548 evar(i) = evar(i) + uvar(i1 + i)/npg
1549 ENDDO
1550 ENDDO
1551 ENDIF
1552 ENDDO
1553 ENDIF !MLW=58
1554 ENDIF
1555c------------------------------------
1556 ELSEIF(ifunc >= 27 .AND. ifunc < 40) THEN
1557c USER VARIABLES from 6 to 60 (mid-layer)
1558c------------------------------------
1559 IF (mlw == 29.OR.mlw == 30.OR.mlw == 31.OR.mlw>=33) THEN
1560 ius = ifunc - 22
1561 IF (nlay > 1) THEN
1562 il = iabs(nlay)/2 + 1
1563 ipt = 1
1564 ELSE
1565 il = 1
1566 ipt = iabs(npt)/2 + 1
1567 ENDIF
1568 nuvar = elbuf_tab(ng)%BUFLY(il)%NVAR_MAT
1569 IF (nuvar > ius .and. npt >= ipt*il) THEN
1570 i1 = ius*nel
1571 DO i=lft,llt
1572 DO ir = 1, nptr
1573 DO is = 1, npts
1574 uvar=>elbuf_tab(ng)%BUFLY(il)%MAT(ir,is,ipt)%VAR
1575 evar(i) = evar(i) + uvar(i1 + i)/npg
1576 ENDDO
1577 ENDDO
1578 ENDDO
1579 ENDIF
1580 ENDIF
1581c------------------------------------
1582 ELSEIF((ifunc > 39 .AND. ifunc < 2040) .OR.
1583 . (ifunc > 2239 .AND. ifunc < 10140)) THEN ! /usri/layer
1584c------------------------------------
1585 IF (ifunc > 39 .and. ifunc < 2040) THEN
1586 ius = (ifunc - 39)/100
1587 ipt = mod((ifunc - 39), 100)
1588 ELSEIF (ifunc > 2239 .AND. ifunc < 10140) THEN
1589 ius = ((ifunc - 2239)/100) + 20
1590 ipt = mod((ifunc - 2239), 100)
1591 ENDIF
1592 IF (ipt == 0) THEN
1593 ipt = 100
1594 ius = ius -1
1595 ENDIF
1596 IF (nlay > 1) THEN
1597 il = ipt
1598 ipt = 1
1599 ELSE
1600 il = 1
1601 ENDIF
1602 nuvar = elbuf_tab(ng)%BUFLY(il)%NVAR_MAT
1603 IF (nuvar > ius .and. (npt >= ipt*il)) THEN
1604 i1 = ius*nel
1605 DO i=lft,llt
1606 DO ir = 1, nptr
1607 DO is = 1, npts
1608 uvar=>elbuf_tab(ng)%BUFLY(il)%MAT(ir,is,ipt)%VAR
1609 evar(i) = evar(i) + uvar(i1 + i)/npg
1610 ENDDO
1611 ENDDO
1612 ENDDO
1613 ENDIF ! USER LAW
1614c------------------------------------
1615 ELSEIF( (ifunc>=10140.AND.ifunc<=10239)
1616 . .OR. ifunc == 10673.OR. ifunc == 10674
1617 . .OR. ifunc == 10675 ) THEN
1618 IF (ifunc == 10673) THEN ! phi Output : Mid Layer
1619 il = iabs(nlay)/2 + 1
1620 ELSEIF (ifunc == 10674) THEN ! phi Output : Upper Layer
1621 il = nlay
1622 ELSEIF (ifunc == 10675) THEN ! phi Output : Lower Layer
1623 il = 1
1624 ELSE
1625 il = ifunc - 10139 ! phi Output : layer 1->99
1626 ENDIF
1627c------------------------------------
1628c
1629 IF (il <= nlay) THEN
1630 bufly => elbuf_tab(ng)%BUFLY(il)
1631 IF (ity == 3) THEN
1632 IF (igtyp == 9 .OR. igtyp == 10 .OR.igtyp == 11 .OR.
1633 . igtyp == 16 .OR. igtyp == 17 .OR. igtyp == 51 .OR.
1634 . igtyp == 52 ) THEN
1635 IF (mlw /= 0 .AND. mlw /= 13) THEN
1636
1637 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
1638 lbuf_dir => elbuf_tab(ng)%BUFLY(il)%LBUF_DIR(1)
1639 DO i=lft,llt
1640 n = i + nft
1641 x21 = x(1,ixc(3,n))-x(1,ixc(2,n))
1642 x32 = x(1,ixc(4,n))-x(1,ixc(3,n))
1643 x34 = x(1,ixc(4,n))-x(1,ixc(5,n))
1644 x41 = x(1,ixc(5,n))-x(1,ixc(2,n))
1645
1646 y21 = x(2,ixc(3,n))-x(2,ixc(2,n))
1647 y32 = x(2,ixc(4,n))-x(2,ixc(3,n))
1648 y34 = x(2,ixc(4,n))-x(2,ixc(5,n))
1649 y41 = x(2,ixc(5,n))-x(2,ixc(2,n))
1650
1651 z21 = x(3,ixc(3,n))-x(3,ixc(2,n))
1652 z32 = x(3,ixc(4,n))-x(3,ixc(3,n))
1653 z34 = x(3,ixc(4,n))-x(3,ixc(5,n))
1654 z41 = x(3,ixc(5,n))-x(3,ixc(2,n))
1655
1656 e1x = (x21+x34)
1657 e1y = (y21+y34)
1658 e1z = (z21+z34)
1659
1660 e2x = (x32+x41)
1661 e2y = (y32+y41)
1662 e2z = (z32+z41)
1663
1664 e3x = e1y*e2z-e1z*e2y
1665 e3y = e1z*e2x-e1x*e2z
1666 e3z = e1x*e2y-e1y*e2x
1667 IF (irep > 0) THEN
1668 rx = e1x
1669 ry = e1y
1670 rz = e1z
1671 sx = e2x
1672 sy = e2y
1673 sz = e2z
1674 ENDIF
1675 IF (ishfram == 0 ) THEN
1676C----- Symmetrical convected user-version 5 (DEFAULT)
1677 suma = e3x*e3x+e3y*e3y+e3z*e3z
1678 suma = one / max(sqrt(suma),em20)
1679 e3x = e3x * suma
1680 e3y = e3y * suma
1681 e3z = e3z * suma
1682C
1683 s1 = e1x*e1x+e1y*e1y+e1z*e1z
1684 s2 = e2x*e2x+e2y*e2y+e2z*e2z
1685 suma = sqrt(s1/s2)
1686 e1x = e1x + (e2y*e3z-e2z*e3y)*suma
1687 e1y = e1y + (e2z*e3x-e2x*e3z)*suma
1688 e1z = e1z + (e2x*e3y-e2y*e3x)*suma
1689C
1690 suma = e1x*e1x+e1y*e1y+e1z*e1z
1691 suma = one / max(sqrt(suma),em20)
1692 e1x = e1x * suma
1693 e1y = e1y * suma
1694 e1z = e1z * suma
1695C
1696 e2x = e3y * e1z - e3z * e1y
1697 e2y = e3z * e1x - e3x * e1z
1698 e2z = e3x * e1y - e3y * e1x
1699 ELSEIF (ishfram == 2) THEN
1700C----- Nononmetrical convected marker-version 4
1701 suma = e2x*e2x+e2y*e2y+e2z*e2z
1702 e1x = e1x*suma + e2y*e3z-e2z*e3y
1703 e1y = e1y*suma + e2z*e3x-e2x*e3z
1704 e1z = e1z*suma + e2x*e3y-e2y*e3x
1705 suma = e1x*e1x+e1y*e1y+e1z*e1z
1706 suma = one/max(sqrt(suma),em20)
1707 e1x = e1x*suma
1708 e1y = e1y*suma
1709 e1z = e1z*suma
1710C
1711 suma = e3x*e3x+e3y*e3y+e3z*e3z
1712 suma = one / max(sqrt(suma),em20)
1713 e3x = e3x * suma
1714 e3y = e3y * suma
1715 e3z = e3z * suma
1716C
1717 e2x = e3y*e1z-e3z*e1y
1718 e2y = e3z*e1x-e3x*e1z
1719 e2z = e3x*e1y-e3y*e1x
1720 suma = e2x*e2x+e2y*e2y+e2z*e2z
1721 suma = one/max(sqrt(suma),em20)
1722 e2x = e2x*suma
1723 e2y = e2y*suma
1724 e2z = e2z*suma
1725 ENDIF
1726 IF (irep >= 1) THEN
1727 aa = lbuf_dir%DIRA(i)
1728 bb = lbuf_dir%DIRA(i+nel)
1729 v1 = aa*rx + bb*sx
1730 v2 = aa*ry + bb*sy
1731 v3 = aa*rz + bb*sz
1732 vr = v1*e1x+ v2*e1y + v3*e1z
1733 vs = v1*e2x+ v2*e2y + v3*e2z
1734 suma=sqrt(vr*vr + vs*vs)
1735 dir1_1 = vr/suma
1736 dir1_2 = vs/suma
1737 ELSE
1738 dir1_1 = lbuf_dir%DIRA(i)
1739 dir1_2 = lbuf_dir%DIRA(i+nel)
1740 ENDIF
1741C
1742 phi =(hundred80/pi)*atan2(dir1_2,dir1_1)
1743 err = (abs(phi) - ninety)/ninety
1744 evar(i) = phi
1745 IF(abs(err) < em02) evar(i) = sign(ninety,phi)
1746 IF(abs(evar(i)) < one) evar(i) = zero
1747 !!EVAR(I) = (HUNDRED80/PI)*ATAN2(DIR1_2,DIR1_1)
1748 ENDDO
1749 ELSE
1750 DO i=lft,llt
1751 n = i + nft
1752 x21 = x(1,ixc(3,n))-x(1,ixc(2,n))
1753 x32 = x(1,ixc(4,n))-x(1,ixc(3,n))
1754 x34 = x(1,ixc(4,n))-x(1,ixc(5,n))
1755 x41 = x(1,ixc(5,n))-x(1,ixc(2,n))
1756
1757 y21 = x(2,ixc(3,n))-x(2,ixc(2,n))
1758 y32 = x(2,ixc(4,n))-x(2,ixc(3,n))
1759 y34 = x(2,ixc(4,n))-x(2,ixc(5,n))
1760 y41 = x(2,ixc(5,n))-x(2,ixc(2,n))
1761
1762 z21 = x(3,ixc(3,n))-x(3,ixc(2,n))
1763 z32 = x(3,ixc(4,n))-x(3,ixc(3,n))
1764 z34 = x(3,ixc(4,n))-x(3,ixc(5,n))
1765 z41 = x(3,ixc(5,n))-x(3,ixc(2,n))
1766
1767 e1x = (x21+x34)
1768 e1y = (y21+y34)
1769 e1z = (z21+z34)
1770
1771 e2x = (x32+x41)
1772 e2y = (y32+y41)
1773 e2z = (z32+z41)
1774
1775 e3x = e1y*e2z-e1z*e2y
1776 e3y = e1z*e2x-e1x*e2z
1777 e3z = e1x*e2y-e1y*e2x
1778 IF (irep > 0) THEN
1779 rx = e1x
1780 ry = e1y
1781 rz = e1z
1782 sx = e2x
1783 sy = e2y
1784 sz = e2z
1785 ENDIF
1786 IF (ishfram == 0 .OR. igtyp == 16 ) THEN
1787C----- Symmetrical convected user-version 5 (DEFAULT)
1788 suma = e3x*e3x+e3y*e3y+e3z*e3z
1789 suma = one / max(sqrt(suma),em20)
1790 e3x = e3x * suma
1791 e3y = e3y * suma
1792 e3z = e3z * suma
1793C
1794 s1 = e1x*e1x+e1y*e1y+e1z*e1z
1795 s2 = e2x*e2x+e2y*e2y+e2z*e2z
1796 suma = sqrt(s1/s2)
1797 e1x = e1x + (e2y*e3z-e2z*e3y)*suma
1798 e1y = e1y + (e2z*e3x-e2x*e3z)*suma
1799 e1z = e1z + (e2x*e3y-e2y*e3x)*suma
1800C
1801 suma = e1x*e1x+e1y*e1y+e1z*e1z
1802 suma = one / max(sqrt(suma),em20)
1803 e1x = e1x * suma
1804 e1y = e1y * suma
1805 e1z = e1z * suma
1806C
1807 e2x = e3y * e1z - e3z * e1y
1808 e2y = e3z * e1x - e3x * e1z
1809 e2z = e3x * e1y - e3y * e1x
1810 ELSEIF (ishfram == 2) THEN
1811C----- Nononmetrical convected marker-version 4
1812 suma = e2x*e2x+e2y*e2y+e2z*e2z
1813 e1x = e1x*suma + e2y*e3z-e2z*e3y
1814 e1y = e1y*suma + e2z*e3x-e2x*e3z
1815 e1z = e1z*suma + e2x*e3y-e2y*e3x
1816 suma = e1x*e1x+e1y*e1y+e1z*e1z
1817 suma = one/max(sqrt(suma),em20)
1818 e1x = e1x*suma
1819 e1y = e1y*suma
1820 e1z = e1z*suma
1821C
1822 suma = e3x*e3x+e3y*e3y+e3z*e3z
1823 suma = one / max(sqrt(suma),em20)
1824 e3x = e3x * suma
1825 e3y = e3y * suma
1826 e3z = e3z * suma
1827C
1828 e2x = e3y*e1z-e3z*e1y
1829 e2y = e3z*e1x-e3x*e1z
1830 e2z = e3x*e1y-e3y*e1x
1831 suma = e2x*e2x+e2y*e2y+e2z*e2z
1832 suma = one/max(sqrt(suma),em20)
1833 e2x = e2x*suma
1834 e2y = e2y*suma
1835 e2z = e2z*suma
1836 ENDIF
1837 IF (irep >= 1) THEN
1838 aa = bufly%DIRA(i)
1839 bb = bufly%DIRA(i+nel)
1840 v1 = aa*rx + bb*sx
1841 v2 = aa*ry + bb*sy
1842 v3 = aa*rz + bb*sz
1843 vr = v1*e1x+ v2*e1y + v3*e1z
1844 vs = v1*e2x+ v2*e2y + v3*e2z
1845 suma=sqrt(vr*vr + vs*vs)
1846 dir1_1 = vr/suma
1847 dir1_2 = vs/suma
1848 ELSE
1849 dir1_1 = bufly%DIRA(i)
1850 dir1_2 = bufly%DIRA(i+nel)
1851 ENDIF
1852C
1853 phi =(hundred80/pi)*atan2(dir1_2,dir1_1)
1854 err = (abs(phi) - ninety)/ninety
1855 evar(i) = phi
1856 IF(abs(err) < em02) evar(i) = sign(ninety,phi)
1857 IF(abs(evar(i)) < one) evar(i) = zero
1858 !!EVAR(I) = (HUNDRED80/PI)*ATAN2(DIR1_2,DIR1_1)
1859 ENDDO
1860 ENDIF ! IDRAPE
1861 ENDIF ! MLW
1862 ENDIF ! IGTYP
1863
1864 ELSEIF (ity == 7) THEN
1865 npg = iparg(48,ng)
1866 IF (igtyp == 9 .OR. igtyp == 10 .OR. igtyp == 11 .OR.
1867 . igtyp == 16 .OR. igtyp == 17 .OR. igtyp == 51 .OR.
1868 . igtyp == 52 ) THEN
1869 IF (mlw /= 0 .AND. mlw /= 13) THEN
1870 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
1871 lbuf_dir => elbuf_tab(ng)%BUFLY(il)%LBUF_DIR(1)
1872 DO i=lft,llt
1873 n = i + nft
1874 x21 = x(1,ixtg(3,n))-x(1,ixtg(2,n))
1875 x31 = x(1,ixtg(4,n))-x(1,ixtg(2,n))
1876 x32 = x(1,ixtg(4,n))-x(1,ixtg(3,n))
1877
1878 y21 = x(2,ixtg(3,n))-x(2,ixtg(2,n))
1879 y31 = x(2,ixtg(4,n))-x(2,ixtg(2,n))
1880 y32 = x(2,ixtg(4,n))-x(2,ixtg(3,n))
1881
1882 z21 = x(3,ixtg(3,n))-x(3,ixtg(2,n))
1883 z31 = x(3,ixtg(4,n))-x(3,ixtg(2,n))
1884 z32 = x(3,ixtg(4,n))-x(3,ixtg(3,n))
1885 IF (irep > 0) THEN
1886 e11 = x21
1887 e12 = y21
1888 e13 = z21
1889 e21 = x31
1890 e22 = y31
1891 e23 = z31
1892 ENDIF
1893 IF(ifram_old ==0 ) THEN
1894 CALL clsconv3(x21,y21,z21,x31,y31,z31,
1895 + e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z)
1896 ELSE
1897 e1x= x21
1898 e1y= y21
1899 e1z= z21
1900 x2l = sqrt(e1x*e1x+e1y*e1y+e1z*e1z)
1901 e1x=e1x/x2l
1902 e1y=e1y/x2l
1903 e1z=e1z/x2l
1904C
1905 e3x=y31*z32-z31*y32
1906 e3y=z31*x32-x31*z32
1907 e3z=x31*y32-y31*x32
1908 sum_ = sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
1909 e3x=e3x/sum_
1910 e3y=e3y/sum_
1911 e3z=e3z/sum _
1912 area = half * sum_
1913 e2x=e3y*e1z-e3z*e1y
1914 e2y=e3z*e1x-e3x*e1z
1915 e2z=e3x*e1y-e3y*e1x
1916 sum_ = sqrt(e2x*e2x+e2y*e2y+e2z*e2z)
1917 e2x=e2x/sum_
1918 e2y=e2y/sum_
1919 e2z=e2z/sum_
1920 END if!(ISH3NFRAM ==0 ) THEN
1921 IF (irep >= 1) THEN
1922 aa = lbuf_dir%DIRA(i)
1923 bb = lbuf_dir%DIRA(i+nel)
1924 v1 = aa*e11 + bb*e21
1925 v2 = aa*e12 + bb*e22
1926 v3 = aa*e13 + bb*e23
1927 vr = v1*e1x + v2*e1y + v3*e1z
1928 vs = v1*e2x + v2*e2y + v3*e2z
1929 suma=sqrt(vr*vr + vs*vs)
1930 dir1_1 = vr/suma
1931 dir1_2 = vs/suma
1932 ELSE
1933 dir1_1 = lbuf_dir%DIRA(i)
1934 dir1_2 = lbuf_dir%DIRA(i+nel)
1935 ENDIF
1936 phi =(hundred80/pi)*atan2(dir1_2,dir1_1)
1937 err = (abs(phi) - ninety)/ninety
1938 evar(i) = phi
1939 IF(abs(err) < em02) evar(i) = sign(ninety,phi)
1940 IF(abs(evar(i)) < one) evar(i) = zero
1941 ENDDO
1942 ELSE ! IDRAPE
1943 DO i=lft,llt
1944 n = i + nft
1945 x21 = x(1,ixtg(3,n))-x(1,ixtg(2,n))
1946 x31 = x(1,ixtg(4,n))-x(1,ixtg(2,n))
1947 x32 = x(1,ixtg(4,n))-x(1,ixtg(3,n))
1948
1949 y21 = x(2,ixtg(3,n))-x(2,ixtg(2,n))
1950 y31 = x(2,ixtg(4,n))-x(2,ixtg(2,n))
1951 y32 = x(2,ixtg(4,n))-x(2,ixtg(3,n))
1952
1953 z21 = x(3,ixtg(3,n))-x(3,ixtg(2,n))
1954 z31 = x(3,ixtg(4,n))-x(3,ixtg(2,n))
1955 z32 = x(3,ixtg(4,n))-x(3,ixtg(3,n))
1956 IF (irep > 0) THEN
1957 e11 = x21
1958 e12 = y21
1959 e13 = z21
1960 e21 = x31
1961 e22 = y31
1962 e23 = z31
1963 ENDIF
1964 IF(ifram_old ==0 ) THEN
1965 CALL clsconv3(x21,y21,z21,x31,y31,z31,
1966 + e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z)
1967 ELSE
1968 e1x= x21
1969 e1y= y21
1970 e1z= z21
1971 x2l = sqrt(e1x*e1x+e1y*e1y+e1z*e1z)
1972 e1x=e1x/x2l
1973 e1y=e1y/x2l
1974 e1z=e1z/x2l
1975C
1976 e3x=y31*z32-z31*y32
1977 e3y=z31*x32-x31*z32
1978 e3z=x31*y32-y31*x32
1979 sum_ = sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
1980 e3x=e3x/sum_
1981 e3y=e3y/sum_
1982 e3z=e3z/sum _
1983 area = half * sum_
1984 e2x=e3y*e1z-e3z*e1y
1985 e2y=e3z*e1x-e3x*e1z
1986 e2z=e3x*e1y-e3y*e1x
1987 sum_ = sqrt(e2x*e2x+e2y*e2y+e2z*e2z)
1988 e2x=e2x/sum_
1989 e2y=e2y/sum_
1990 e2z=e2z/sum_
1991 END if!(ISH3NFRAM ==0 ) THEN
1992 IF (irep >= 1) THEN
1993 aa = bufly%DIRA(i)
1994 bb = bufly%DIRA(i+nel)
1995 v1 = aa*e11 + bb*e21
1996 v2 = aa*e12 + bb*e22
1997 v3 = aa*e13 + bb*e23
1998 vr = v1*e1x + v2*e1y + v3*e1z
1999 vs = v1*e2x + v2*e2y + v3*e2z
2000 suma=sqrt(vr*vr + vs*vs)
2001 dir1_1 = vr/suma
2002 dir1_2 = vs/suma
2003 ELSE
2004 dir1_1 = bufly%DIRA(i)
2005 dir1_2 = bufly%DIRA(i+nel)
2006 ENDIF
2007 phi =(hundred80/pi)*atan2(dir1_2,dir1_1)
2008 err = (abs(phi) - ninety)/ninety
2009 evar(i) = phi
2010 IF(abs(err) < em02) evar(i) = sign(ninety,phi)
2011 IF(abs(evar(i)) < one) evar(i) = zero
2012 ENDDO
2013 ENDIF ! IDRAPE
2014 ENDIF
2015 ENDIF
2016 ENDIF
2017 ELSE
2018 DO i=lft,llt
2019 evar(i) = zero
2020 ENDDO
2021 ENDIF ! IL
2022c------------------------------------
2023 ELSEIF (ifunc == 2040 .AND. mlw /= 15 .AND. mlw /= 25) THEN ! EPSP/UPPER
2024c------------------------------------
2025 IF (nlay > 1) THEN
2026 il = max(1,npt)
2027 ipt = 1
2028 ELSE
2029 il = 1
2030 ipt = max(1,npt)
2031 ENDIF
2032 bufly => elbuf_tab(ng)%BUFLY(il)
2033 IF (bufly%L_PLA > 0) THEN
2034 IF (npg > 1) THEN
2035 IF (igtyp == 51 .OR. igtyp == 52) ipt = bufly%NPTT
2036 DO i=lft,llt
2037 DO ir=1,nptr
2038 DO is=1,npts
2039 lbuf => bufly%LBUF(ir,is,ipt)
2040 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
2041 ENDDO
2042 ENDDO
2043 ENDDO
2044 ELSE
2045 IF (igtyp == 51 .OR. igtyp == 52) ipt = bufly%NPTT
2046 DO i=lft,llt
2047 evar(i) = abs(bufly%LBUF(1,1,ipt)%PLA(i))
2048 ENDDO
2049 ENDIF
2050 ENDIF
2051c------------------------------------
2052 ELSEIF (ifunc == 2041 .AND. mlw /= 15 .AND. mlw /= 25) THEN ! EPSP/LOWER
2053c------------------------------------
2054 bufly => elbuf_tab(ng)%BUFLY(1)
2055 IF (bufly%L_PLA > 0) THEN
2056 IF (npg > 1) THEN
2057 DO i=lft,llt
2058 DO ir=1,nptr
2059 DO is=1,npts
2060 lbuf => bufly%LBUF(ir,is,1)
2061 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
2062 ENDDO
2063 ENDDO
2064 ENDDO
2065 ELSE
2066 DO i=lft,llt
2067 evar(i) = abs(bufly%LBUF(1,1,1)%PLA(i))
2068 ENDDO
2069 ENDIF
2070 ENDIF
2071c------------------------------------
2072 ELSEIF (ifunc > 2041 .AND. ifunc < 2142 .AND. mlw /= 15 .AND. mlw /= 25) THEN ! anim/shell/epsp/layer
2073c------------------------------------
2074 ilay = mod((ifunc - 2041), 100)
2075 IF (ilay == 0) ilay = 100
2076 IF ((ilay <= nlay .or. ilay <= mpt) .and. gbuf%G_PLA > 0) THEN
2077 IF (npt == 0) THEN
2078 il = 1
2079 ipt = 1
2080 ELSEIF (nlay > 1) THEN
2081 il = ilay
2082 ipt = 1
2083 ELSE
2084 il = 1
2085 ipt = ilay
2086 ENDIF
2087 bufly => elbuf_tab(ng)%BUFLY(il)
2088 IF (bufly%L_PLA > 0) THEN
2089 IF (npg > 1) THEN
2090 IF (igtyp == 51 .OR. igtyp == 52) THEN
2091C PROP/51 : more than 1 integration point by layer: average value per layer
2092 nptt = bufly%NPTT
2093 npgt = npg*nptt
2094 DO i=lft,llt
2095 DO it=1,nptt
2096 DO ir=1,nptr
2097 DO is=1,npts
2098 lbuf => bufly%LBUF(ir,is,it)
2099 evar(i) = evar(i) + abs(lbuf%PLA(i))/npgt
2100 ENDDO
2101 ENDDO
2102 ENDDO
2103 ENDDO
2104 ELSE
2105 DO i=lft,llt
2106 DO ir=1,nptr
2107 DO is=1,npts
2108 lbuf => bufly%LBUF(ir,is,ipt)
2109 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
2110 ENDDO
2111 ENDDO
2112 ENDDO
2113 ENDIF
2114 ELSE
2115 IF (igtyp == 51 .OR. igtyp == 52) THEN
2116 nptt = bufly%NPTT
2117 DO i=lft,llt
2118 DO it=1,nptt
2119 evar(i) = evar(i) + abs(bufly%LBUF(1,1,it)%PLA(i))/nptt
2120 ENDDO
2121 ENDDO
2122 ELSE
2123 DO i=lft,llt
2124 evar(i) = abs(bufly%LBUF(1,1,ipt)%PLA(i))
2125 ENDDO
2126 ENDIF
2127 ENDIF
2128 ENDIF
2129 ENDIF
2130c------------------------------------
2131 ELSEIF (ifunc == 10253.OR.ifunc == 10254.OR.ifunc == 10255)THEN
2132C NXT failure model output
2133C Max values: Damage factor + normalized stress
2134c------------------------------------
2135 IF (ifunc == 10253) THEN ! output Damage
2136 DO il=1,nlay
2137 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2138 DO is=1,npts
2139 DO it=1,nptt
2140 DO ir=1,nptr
2141 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,it)
2142 DO ifail=1,nfail
2143 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2144 DO i=lft,llt
2145 evar(i) = max(evar(i),fbuf%FLOC(ifail)%DAM(i))
2146 ENDDO
2147 ENDIF
2148 ENDDO
2149 ENDDO
2150 ENDDO
2151 ENDDO
2152 ENDDO
2153 ELSEIF (ifunc == 10254) THEN ! output Sig1
2154 DO il=1,nlay
2155 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2156 DO is=1,npts
2157 DO it=1,nptt
2158 DO ir=1,nptr
2159 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,it)
2160 DO ifail=1,nfail
2161 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2162 nvarf = fbuf%FLOC(ifail)%NVAR
2163 DO i=lft,llt
2164 var = fbuf%FLOC(ifail)%VAR(nvarf*(i-1)+1) ! Sig1
2165 evar(i) = max(evar(i), var)
2166 ENDDO
2167 ENDIF
2168 ENDDO
2169 ENDDO
2170 ENDDO
2171 ENDDO
2172 ENDDO
2173 ELSEIF (ifunc == 10255) THEN ! output Sig2
2174 DO il=1,nlay
2175 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2176 DO is=1,npts
2177 DO it=1,nptt
2178 DO ir=1,nptr
2179 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,it)
2180 DO ifail=1,nfail
2181 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2182 nvarf = fbuf%FLOC(ifail)%NVAR
2183 DO i=lft,llt
2184 var = fbuf%FLOC(ifail)%VAR(nvarf*(i-1)+2) ! Sig2
2185 evar(i) = max(evar(i), var)
2186 ENDDO
2187 ENDIF
2188 ENDDO
2189 ENDDO
2190 ENDDO
2191 ENDDO
2192 ENDDO
2193 ENDIF
2194C--------------------------------------------------
2195 ELSE IF (ifunc >= 10360 .and. ifunc <= 10668) THEN
2196C--------------------------------------------------
2197C NXT failure model output par layer : upper/lower/membrane
2198C Damage factor + normalized stress
2199c---
2200 IF (ifunc == 10360) THEN ! Damage factor - upper
2201 IF (nlay > 1) THEN
2202 il = nlay
2203 ipt = 1
2204 ELSE
2205 il = 1
2206 ipt = nptt
2207 ENDIF
2208cc DO IL=1,NLAY ! loop unneeded for the upper layer only
2209 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2210 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
2211 DO is=1,npts
2212 DO ir=1,nptr
2213 DO it=1,nptt
2214 ipt = it
2215 IF (nlay == 1) ipt = nptt
2216 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2217 DO ifail=1,nfail
2218 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2219 DO i=lft,llt
2220 evar(i) = max(evar(i),fbuf%FLOC(ifail)%DAM(i))
2221 ENDDO
2222 ENDIF
2223 ENDDO
2224 ENDDO
2225 ENDDO
2226 ENDDO
2227cc ENDDO
2228 ELSEIF (ifunc == 10361) THEN ! Damage factor - lower
2229 ipt = 1
2230 il = 1
2231cc DO IL=1,NLAY ! loop unneeded for the lower layer only
2232 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2233 DO is=1,npts
2234 DO ir=1,nptr
2235 DO it=1,nptt
2236 ipt = it
2237 IF (nlay == 1) ipt = 1
2238 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2239 DO ifail=1,nfail
2240 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2241 DO i=lft,llt
2242 evar(i) = max(evar(i),fbuf%FLOC(ifail)%DAM(i))
2243 ENDDO
2244 ENDIF
2245 ENDDO
2246 ENDDO
2247 ENDDO
2248 ENDDO
2249cc ENDDO
2250 ELSEIF (ifunc == 10362) THEN ! Damage factor - membrane
2251 IF (nlay > 1) THEN
2252 il = iabs(nlay) / 2
2253 ipt = 1
2254 ELSE
2255 il = 1
2256 ipt = iabs(nptt) / 2
2257 ENDIF
2258cc DO IL=1,NLAY ! loop unneeded for the membrane layer only
2259 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2260 DO is=1,npts
2261 DO ir=1,nptr
2262 DO it=1,nptt
2263 ipt = it
2264 IF (nlay == 1) ipt = iabs(nptt) / 2
2265 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2266 DO ifail=1,nfail
2267 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2268 DO i=lft,llt
2269 evar(i) = max(evar(i),fbuf%FLOC(ifail)%DAM(i))
2270 ENDDO
2271 ENDIF
2272 ENDDO
2273 ENDDO
2274 ENDDO
2275 ENDDO
2276cc ENDDO
2277 ELSEIF (ifunc == 10363) THEN ! Sig1 - upper
2278 IF (nlay > 1) THEN
2279 il = nlay
2280 ipt = 1
2281 ELSE
2282 il = 1
2283 ipt = nptt
2284 ENDIF
2285 DO il=1,nlay
2286 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2287 DO is=1,npts
2288 DO ir=1,nptr
2289 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2290 DO ifail=1,nfail
2291 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2292 nvarf = fbuf%FLOC(ifail)%NVAR
2293 DO i=lft,llt
2294 var = fbuf%FLOC(ifail)%VAR(nvarf*(i-1)+1) ! Sig1
2295 evar(i) = max(evar(i), var)
2296 ENDDO
2297 ENDIF
2298 ENDDO
2299 ENDDO
2300 ENDDO
2301 ENDDO
2302 ELSEIF (ifunc == 10364) THEN ! Sig1 - lower
2303 ipt = 1
2304 il = 1
2305 DO il=1,nlay
2306 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2307 DO is=1,npts
2308 DO ir=1,nptr
2309 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2310 DO ifail=1,nfail
2311 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2312 nvarf = fbuf%FLOC(ifail)%NVAR
2313 DO i=lft,llt
2314 var = fbuf%FLOC(ifail)%VAR(nvarf*(i-1)+1) ! Sig1
2315 evar(i) = max(evar(i), var)
2316 ENDDO
2317 ENDIF
2318 ENDDO
2319 ENDDO
2320 ENDDO
2321 ENDDO
2322 ELSEIF (ifunc == 10365) THEN ! Sig1 - membrane
2323 IF (nlay > 1) THEN
2324 il = nlay / 2
2325 ipt = 1
2326 ELSE
2327 il = 1
2328 ipt = nptt / 2
2329 ENDIF
2330 DO il=1,nlay
2331 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2332 DO is=1,npts
2333 DO ir=1,nptr
2334 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2335 DO ifail=1,nfail
2336 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2337 nvarf = fbuf%FLOC(ifail)%NVAR
2338 DO i=lft,llt
2339 var = fbuf%FLOC(ifail)%VAR(nvarf*(i-1)+1) ! Sig1
2340 evar(i) = max(evar(i), var)
2341 ENDDO
2342 ENDIF
2343 ENDDO
2344 ENDDO
2345 ENDDO
2346 ENDDO
2347 ELSEIF (ifunc == 10366) THEN ! Sig2 - upper
2348 IF (nlay > 1) THEN
2349 il = nlay
2350 ipt = 1
2351 ELSE
2352 il = 1
2353 ipt = nptt
2354 ENDIF
2355 DO il=1,nlay
2356 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2357 DO is=1,npts
2358 DO ir=1,nptr
2359 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2360 DO ifail=1,nfail
2361 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check nxt model
2362 nvarf = fbuf%FLOC(ifail)%NVAR
2363 DO i=lft,llt
2364 var = fbuf%FLOC(ifail)%VAR(nvarf*(i-1)+2) ! Sig2
2365 evar(i) = max(evar(i), var)
2366 ENDDO
2367 ENDIF
2368 ENDDO
2369 ENDDO
2370 ENDDO
2371 ENDDO
2372 ELSEIF (ifunc == 10367) THEN ! Sig2 - lower
2373 ipt = 1
2374 il = 1
2375 DO il=1,nlay
2376 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2377 DO is=1,npts
2378 DO ir=1,nptr
2379 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2380 DO ifail=1,nfail
2381 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2382 nvarf = fbuf%FLOC(ifail)%NVAR
2383 DO i=lft,llt
2384 var = fbuf%FLOC(ifail)%VAR(nvarf*(i-1)+2) ! Sig2
2385 evar(i) = max(evar(i), var)
2386 ENDDO
2387 ENDIF
2388 ENDDO
2389 ENDDO
2390 ENDDO
2391 ENDDO
2392 ELSEIF (ifunc == 10368) THEN ! Sig2 - membrane
2393 IF (nlay > 1) THEN
2394 il = nlay / 2
2395 ipt = 1
2396 ELSE
2397 il = 1
2398 ipt = nptt / 2
2399 ENDIF
2400 DO il=1,nlay
2401 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2402 DO is=1,npts
2403 DO ir=1,nptr
2404 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
2405 DO ifail=1,nfail
2406 IF (fbuf%FLOC(ifail)%ILAWF == 25) THEN ! check NXT model
2407 nvarf = fbuf%FLOC(ifail)%NVAR
2408 DO i=lft,llt
2409 var = fbuf%FLOC(ifail)%VAR(nvarf*(i-1)+2) ! Sig2
2410 evar(i) = max(evar(i), var)
2411 ENDDO
2412 ENDIF
2413 ENDDO
2414 ENDDO
2415 ENDDO
2416 ENDDO
2417 ENDIF
2418C--------------------------------------------------
2419 ELSE IF (ifunc == 2142) THEN ! FAIL
2420c------------------------------------
2421 IF (igtyp == 10.OR.igtyp == 11.OR.igtyp == 17.OR. igtyp == 51
2422 . .OR. igtyp == 52) THEN
2423 failg = 0
2424 DO i=lft,llt
2425 dam1(i)=zero
2426 dam2(i)=zero
2427 wpla(i)=zero
2428 fail(i)=zero
2429 END DO
2430 IF(ity == 3)THEN
2431 DO i=lft,llt
2432 mat(i)=ixc(1,nft+i)
2433 pid(i)=ixc(6,nft+i)
2434 END DO
2435 ELSE
2436 DO i=lft,llt
2437 mat(i)=ixtg(1,nft+i)
2438 pid(i)=ixtg(5,nft+i)
2439 END DO
2440 END IF
2441 IF (igtyp == 11) THEN
2442 ipmat = 100
2443 DO n=1,npt
2444 iadr = (n-1)*nel
2445 DO i=lft,llt
2446 j = iadr+i
2447 matly(j) = igeo(ipmat+n,pid(i))
2448 END DO
2449 END DO
2450 ELSEIF (igtyp == 10) THEN
2451 DO n=1,npt
2452 iadr = (n-1)*nel
2453 DO i=lft,llt
2454 j = iadr+i
2455 matly(j)=mat(i)
2456 END DO
2457 END DO
2458 ELSEIF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
2459 ipmat = 2 + nlay
2460! old stack organisation IPMAT = 300
2461 DO n=1,nlay
2462 iadr = (n-1)*nel
2463 DO i=lft,llt
2464 j = iadr+i
2465 matly(j) = stack%IGEO(ipmat+n,isubstack)
2466 END DO
2467 END DO
2468 END IF
2469c
2470 IF (ihbe == 11) THEN ! Batoz shell
2471 DO il=1,nlay
2472 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
2473 bufly => elbuf_tab(ng)%BUFLY(il)
2474 iadr = (il-1)*nel
2475 DO it=1,nptt
2476 DO i=lft,llt
2477 wpla(i) = zero
2478 ENDDO
2479 failg = 0
2480 DO ir=1,nptr
2481 DO is=1,npts
2482 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
2483 offl => lbuf%OFF
2484 IF (bufly%L_DAM > 0 .OR. bufly%L_OFF > 0 ) THEN
2485 DO i=lft,llt
2486 j = iadr + i
2487 IF(ipm(2,matly(j)) == 15) THEN
2488 dam1(i)=lbuf%DAM(jj(1)+i)
2489 dam2(i)=lbuf%DAM(jj(2)+i)
2490 wpla(i) = wpla(i) + abs(lbuf%PLA(i))/npg
2491 dmax(i) = pm(64,matly(j))
2492 wpmax(i)= pm(41,matly(j))
2493 IF (dam1(i) >= dmax(i).OR.dam2(i) >= dmax(i)
2494 . .OR.wpla(i) < zero.OR.wpla(i) >= wpmax(i)
2495 . .OR.offl(i) < one) failg(i) = failg(i) + 1
2496 IF (failg(i) == npg ) THEN
2497 fail(i) = fail(i) + one
2498 ENDIF
2499 ELSEIF (ipm(2,matly(j)) == 25) THEN
2500 dam1(i)=lbuf%DMG(jj(2)+i)
2501 dam2(i)=lbuf%DMG(jj(3)+i)
2502 wpla(i) = wpla(i) + abs(lbuf%PLA(i))/npg
2503 dmax(i) = pm(64,matly(j))
2504 wpmax(i)= pm(41,matly(j))
2505 IF (dam1(i) >= dmax(i).OR.dam2(i) >= dmax(i)
2506 . .OR.wpla(i) < zero.OR.wpla(i) >= wpmax(i)
2507 . .OR.offl(i) < one) failg(i) = failg(i) + 1
2508 IF (failg(i) == npg ) THEN
2509 fail(i) = fail(i) + one
2510 ENDIF
2511 ELSE
2512 IF (offl(i) < one) failg(i)= failg(i) + 1
2513 IF (failg(i) == npg ) THEN
2514 fail(i) = fail(i) + one
2515 ENDIF
2516 ENDIF ! law25
2517 ENDDO
2518 ENDIF
2519 ENDDO
2520 ENDDO
2521 ENDDO
2522 ENDDO ! DO IL=1,NLAY
2523 DO i=lft,llt
2524 evar(i) = fail(i)
2525 ENDDO
2526 ELSE ! all shells
2527 DO il=1,nlay
2528 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
2529 bufly => elbuf_tab(ng)%BUFLY(il)
2530
2531 iadr = (il-1)*nel
2532 DO it=1,nptt
2533 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(1,1,it)
2534 offl => lbuf%OFF
2535 IF (bufly%L_DAM > 0 .OR.bufly%L_OFF > 0 ) THEN
2536 DO i=lft,llt
2537 j = iadr + i
2538 IF (ipm(2,matly(j)) == 15) THEN
2539 dam1(i) = lbuf%DAM(jj(1)+i)
2540 dam2(i) = lbuf%DAM(jj(2)+i)
2541 wpla(i) = abs(lbuf%PLA(i))
2542 dmax(i) = pm(64,matly(j))
2543 wpmax(i)= pm(41,matly(j))
2544 IF (dam1(i) >= dmax(i).OR.dam2(i) >= dmax(i).OR.
2545 . wpla(i) < zero.OR.wpla(i) >= wpmax(i) .OR.
2546 . offl(i) < one ) fail(i) = fail(i) + one
2547 ELSEIF (ipm(2,matly(j)) == 25) THEN
2548 dam1(i) = lbuf%DMG(jj(2)+i)
2549 dam2(i) = lbuf%DMG(jj(3)+i)
2550 wpla(i) = abs(lbuf%PLA(i))
2551 dmax(i) = pm(64,matly(j))
2552 wpmax(i)= pm(41,matly(j))
2553 IF (dam1(i) >= dmax(i).OR.dam2(i) >= dmax(i).OR.
2554 . wpla(i) < zero.OR.wpla(i) >= wpmax(i) .OR.
2555 . offl(i) < one ) fail(i) = fail(i) + one
2556 ELSE
2557 IF (offl(i) < one ) fail(i) = fail(i) + one
2558 ENDIF
2559 ENDDO
2560 ENDIF
2561 ENDDO
2562 ENDDO ! DO IL=1,NLAY
2563 DO i=lft,llt
2564 evar(i) = fail(i)
2565 ENDDO
2566!! ELSE
2567!! IF (IFAILA == 1) THEN
2568!! DO I=LFT,LLT
2569!! OFF = GBUF%OFF(I)
2570!! IF (OFF<ZERO) THEN
2571!! FAIL(I)= -ONE
2572! ELSEIF (OFF>ZERO) THEN
2573!! FAIL(I)= ZERO
2574!! ELSE
2575!! FAIL(I)= ONE
2576!! END IF
2577!! EVAR(I)=FAIL(I)
2578!! END DO
2579!! ELSE
2580!! DO I=LFT,LLT
2581
2582!! ENDDO
2583!! ENDIF
2584 ENDIF ! IHBE
2585c
2586 ELSE !
2587!! Should be clarified for others Pid.
2588!! IF (ITY == 3) THEN
2589!! DO I=LFT,LLT
2590!! MAT(I)=IXC(1,NFT+I)
2591!! PID(I)=IXC(6,NFT+I)
2592!! END DO
2593!! ELSE
2594!! DO I=LFT,LLT
2595!! MAT(I)=IXTG(1,NFT+I)
2596!! PID(I)=IXTG(5,NFT+I)
2597!! END DO
2598!! END IF
2599!! IF (IGTYP == 11 .OR. IGTYP == 16) THEN
2600!! IPMAT = 100
2601!! DO N=1,NPT
2602!! IADR = (N-1)*NEL
2603!! DO I=LFT,LLT
2604!! J = IADR+I
2605!! MATLY(J) = IGEO(IPMAT+N,PID(I))
2606!! END DO
2607!! END DO
2608!! ENDIF
2609!! IF (IFAILURE==0 .OR.(IFAILURE /=0 .AND.IFAILA==1)) THEN
2610!! DO I=LFT,LLT
2611!! OFF = GBUF%OFF(I)
2612!! IF (OFF<ZERO)THEN
2613!! FAIL(I)= -ONE
2614!! ELSEIF(OFF>ZERO)THEN
2615!! FAIL(I)= ZERO
2616!! ELSE
2617!! FAIL(I)= ONE
2618!! END IF
2619!! EVAR(I)=FAIL(I)
2620!! END DO
2621!! ELSE
2622!! DO I=LFT,LLT
2623
2624!! ENDDO
2625!! ENDIF
2626 ENDIF ! MLW , IGTYP
2627c---------------------
2628 ELSE IF (ifunc >= 10256 .and. ifunc <= 10359) THEN
2629C---------------------Damage Output : ALL Layers
2630 IF (ifunc == 10257) THEN ! Damage Output : Upper Layer
2631 ipt = npt
2632 ELSEIF (ifunc == 10258) THEN ! Damage Output : Lower Layer
2633 ipt = 1
2634 ELSEIF (ifunc == 10259) THEN ! Damage Output : Membrane Layer
2635 ipt = iabs(npt)/2 + 1
2636 ELSEIF (ifunc >= 10260 .AND. ifunc <= 10359) THEN ! Damage Output : layer IPT
2637 ipt = mod((ifunc - 10259), 100)
2638 IF (ipt == 0) ipt = 100
2639 ENDIF
2640c
2641 DO i = lft, llt
2642 evar(i) = zero
2643 ENDDO
2644c
2645 IF(ifailure > 0) THEN
2646C-------
2647 IF (ifunc == 10256) THEN
2648! Element Damage Output --> Average value over all layers :
2649! for each layer --> get mean value over it's NPTT (for each NPTT --> max over all pts Gauss, for all failure models)
2650C-------
2651 IF (nlay > 1) THEN
2652 DO i = lft,llt
2653 DO n = 1,nlay
2654 nptt = elbuf_tab(ng)%BUFLY(n)%NPTT
2655 dmgmx_ly = zero
2656 DO it = 1,nptt
2657 dmgmx = zero
2658 DO ir = 1,nptr
2659 DO is = 1,npts
2660 fbuf => elbuf_tab(ng)%BUFLY(n)%FAIL(ir,is,it)
2661 DO ifail = 1,elbuf_tab(ng)%BUFLY(n)%NFAIL
2662 dmgmx = max(dmgmx,fbuf%FLOC(ifail)%DAMMX(i))
2663 ENDDO
2664 ENDDO
2665 ENDDO
2666 dmgmx_ly = dmgmx_ly + dmgmx / nptt
2667 ENDDO ! DO IT = 1,NPTT
2668 evar(i) = evar(i) + dmgmx_ly
2669 ENDDO ! N = 1,NLAY
2670 evar(i) = evar(i) / nlay
2671 ENDDO ! DO I = LFT,LLT
2672 ELSEIF (mpt > 0) THEN ! NLAY = 1
2673 nptt = elbuf_tab(ng)%BUFLY(1)%NPTT
2674 DO i = lft,llt
2675 DO it = 1,nptt
2676 dmgmx = zero
2677 DO ir = 1,nptr
2678 DO is = 1,npts
2679 fbuf => elbuf_tab(ng)%BUFLY(1)%FAIL(ir,is,it)
2680 DO ifail = 1,elbuf_tab(ng)%BUFLY(1)%NFAIL
2681 dmgmx = max(dmgmx, fbuf%FLOC(ifail)%DAMMX(i))
2682 ENDDO
2683 ENDDO
2684 ENDDO
2685 evar(i) = evar(i) + dmgmx
2686 ENDDO ! N = 1,NPTT
2687 evar(i) = evar(i) / nptt
2688 ENDDO ! I = LFT,LLT
2689 ENDIF ! IF (NLAY > 1)
2690C-------
2691 ELSEIF (npt >= ipt) THEN
2692! Layer Damage Output (UPPER, LOWER, MEMB, ILAY) :
2693! for layer (ILAY) --> get mean value over it's NPTT (for each NPTT --> max over all pts Gauss, for all failure models)
2694C-------
2695 IF (nlay > 1 .AND. ipt <= nlay) THEN
2696 nptt = elbuf_tab(ng)%BUFLY(ipt)%NPTT
2697 DO i = lft,llt
2698 DO it = 1,nptt
2699 dmgmx = zero
2700 DO ir = 1,nptr
2701 DO is = 1,npts
2702 fbuf => elbuf_tab(ng)%BUFLY(ipt)%FAIL(ir,is,it)
2703 DO ifail = 1,elbuf_tab(ng)%BUFLY(ipt)%NFAIL
2704 dmgmx = max(dmgmx,fbuf%FLOC(ifail)%DAMMX(i))
2705 ENDDO
2706 ENDDO
2707 ENDDO
2708 evar(i) = evar(i) + dmgmx
2709 ENDDO ! DO IT = 1,NPTT
2710 evar(i) = evar(i) / nptt
2711 ENDDO ! I = LFT,LLT
2712 ELSEIF (mpt > 0) THEN ! NLAY = 1
2713 DO i = lft,llt
2714 DO ir = 1, nptr
2715 DO is = 1, npts
2716 fbuf => elbuf_tab(ng)%BUFLY(1)%FAIL(ir,is,ipt)
2717 DO ifail = 1, elbuf_tab(ng)%BUFLY(1)%NFAIL
2718 evar(i) = max(evar(i), fbuf%FLOC(ifail)%DAMMX(i))
2719 ENDDO
2720 ENDDO
2721 ENDDO
2722 ENDDO ! I = LFT,LLT
2723 ENDIF ! IF (NLAY > 1 .AND. IPT <= NLAY)
2724 ENDIF ! Damage IF (IFUNC == 10256)
2725 ENDIF ! IFAILURE
2726C
2727C for outp of dam inside law25
2728C
2729 IF(mlw == 25 .AND. (igtyp == 10 .OR. igtyp == 11 .OR.
2730 . igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52 )) THEN
2731 IF(ity == 3)THEN
2732 DO i=lft,llt
2733 mat(i)=ixc(1,nft+i)
2734 pid(i)=ixc(6,nft+i)
2735 END DO
2736 ELSE
2737 DO i=lft,llt
2738 mat(i)=ixtg(1,nft+i)
2739 pid(i)=ixtg(5,nft+i)
2740 END DO
2741 END IF
2742 IF (igtyp == 11) THEN
2743 ipmat = 100
2744 DO n=1,npt
2745 iadr = (n-1)*nel
2746 DO i=lft,llt
2747 j = iadr+i
2748 matly(j) = igeo(ipmat+n,pid(i))
2749 END DO
2750 END DO
2751 ELSEIF (igtyp == 10) THEN
2752 DO n=1,npt
2753 iadr = (n-1)*nel
2754 DO i=lft,llt
2755 j = iadr+i
2756 matly(j)=mat(i)
2757 END DO
2758 END DO
2759 ELSEIF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
2760 ipmat = 2 + nlay
2761 DO n=1,nlay
2762 iadr = (n-1)*nel
2763 DO i=lft,llt
2764 j = iadr+i
2765 matly(j) = stack%IGEO(ipmat+n,isubstack)
2766 END DO
2767 END DO
2768 END IF
2769C
2770CC IF (IHBE == 11) THEN Batoz NPG = 4
2771 IF (ifunc == 10256) THEN
2772 DO i=lft,llt
2773 ve(1:5) = zero
2774 DO il=1,nlay
2775 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
2776 bufly => elbuf_tab(ng)%BUFLY(il)
2777 iadr = (il-1)*nel
2778 j = iadr + i
2779 vly(1:5) = zero
2780 DO it=1,nptt
2781 vg(1:5)= zero
2782 DO ir=1,nptr
2783 DO is=1,npts
2784 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
2785 dmax(i) = one/pm(64,matly(j))
2786 wpmax(i)= one/pm(41,matly(j))
2787 epst1(i)= pm(60,matly(j))
2788 epst2(i)= pm(61,matly(j))
2789 epsf1(i)= one/pm(98,matly(j))
2790 epsf2(i)= one/pm(99,matly(j))
2791C
2792 vg(1) = max(vg(1),lbuf%DMG(jj(2)+i)*dmax(i))
2793 vg(2) = max(vg(2),lbuf%DMG(jj(3)+i)*dmax(i))
2794 vg(3)=max(vg(3),abs(lbuf%PLA(i))*wpmax(i))
2795 IF(lbuf%CRAK(jj(1)+i) > zero) vg(4)= max(vg(4),
2796 . (lbuf%CRAK(jj(1)+i)+epst1(i))*epsf1(i))
2797 IF(lbuf%CRAK(jj(2)+i) > zero )vg(5) = max(vg(5),
2798 . (lbuf%CRAK(jj(2)+i)+epst2(i))*epsf2(i))
2799 ENDDO
2800 ENDDO
2801 vly(1) = vly(1) + vg(1)
2802 vly(2) = vly(2) + vg(2)
2803 vly(3) = vly(3) + vg(3)
2804 vly(4) = vly(4) + vg(4)
2805 vly(5) = vly(5) + vg(5)
2806 ENDDO ! NPTT
2807 ve(1) = ve(1) + vly(1)/nptt
2808 ve(2) = ve(2) + vly(2)/nptt
2809 ve(3) = ve(3) + vly(3)/nptt
2810 ve(4) = ve(4) + vly(4)/nptt
2811 ve(5) = ve(5) + vly(5)/nptt
2812 ENDDO ! DO IL=1,NLAY
2813 ve(1) = ve(1)/nlay
2814 ve(2) = ve(2)/nlay
2815 ve(3) = ve(3)/nlay
2816 ve(4) = ve(4)/nlay
2817 ve(5) = ve(5)/nlay
2818 evar(i) = max(evar(i),ve(1),ve(2),ve(3),
2819 . ve(4),ve(5))
2820 ENDDO ! I=LFT,JLT
2821 ELSEIF(ipt <= nlay) THEN
2822!! DO IL=1,NLAY ! IPT = IL (layer)
2823 DO i=lft,llt
2824 nptt = elbuf_tab(ng)%BUFLY(ipt)%NPTT
2825 bufly => elbuf_tab(ng)%BUFLY(ipt)
2826 iadr = (ipt - 1)*nel
2827 j = iadr + i
2828 vly(1:5) = zero
2829 DO it=1,nptt
2830 vg(1:5)= zero
2831 DO ir=1,nptr
2832 DO is=1,npts
2833 lbuf => elbuf_tab(ng)%BUFLY(ipt)%LBUF(ir,is,it)
2834 dmax(i) = one/pm(64,matly(j))
2835 wpmax(i)= one/pm(41,matly(j))
2836 epst1(i)= pm(60,matly(j))
2837 epst2(i)= pm(61,matly(j))
2838 epsf1(i)= one/pm(98,matly(j))
2839 epsf2(i)= one/pm(99,matly(j))
2840C
2841 vg(1) = max(vg(1),lbuf%DMG(jj(2)+i)*dmax(i))
2842 vg(2) = max(vg(2),lbuf%DMG(jj(3)+i)*dmax(i))
2843 vg(3)= max(vg(3),abs(lbuf%PLA(i))*wpmax(i))
2844 IF(lbuf%CRAK(jj(1)+i) > zero) vg(4)= max(vg(4),
2845 . (lbuf%CRAK(jj(1)+i)+epst1(i))*epsf1(i))
2846 IF(lbuf%CRAK(jj(2)+i) > zero )vg(5) = max(vg(5),
2847 . (lbuf%CRAK(jj(2)+i)+epst2(i))*epsf2(i))
2848 ENDDO
2849 ENDDO
2850 vly(1) =vly(1) + vg(1)
2851 vly(2) =vly(2) + vg(2)
2852 vly(3) =vly(3) + vg(3)
2853 vly(4) =vly(4) + vg(4)
2854 vly(5) =vly(5) + vg(5)
2855 ENDDO ! NPTT
2856 vly(1) =vly(1)/nptt
2857 vly(2) =vly(2)/nptt
2858 vly(3) =vly(3)/nptt
2859 vly(4) =vly(4)/nptt
2860 vly(5) =vly(5)/nptt
2861C
2862 evar(i) = max(evar(i),vly(1),vly(2),vly(3),
2863 . vly(4),vly(5))
2864 ENDDO ! I=LFT,JLT
2865 ENDIF
2866 ENDIF ! law 25 + SHell Composite PID
2867
2868C-------------------------------------
2869 ELSE IF (ifunc == 10670) THEN
2870C FAIL TIME : ELEMENT DELETED
2871C-------------------------------------
2872 DO i = lft, llt
2873 evar(i) = zero
2874 ENDDO
2875c
2876 DO il=1,nlay
2877 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
2878 DO is=1,npts
2879 DO it=1,nptt
2880 DO ir=1,nptr
2881 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,it)
2882 DO ifail=1,nfail
2883 DO i=lft,llt
2884 evar(i) = max(evar(i),fbuf%FLOC(ifail)%TDEL(i))
2885 ENDDO
2886 ENDDO
2887 ENDDO
2888 ENDDO
2889 ENDDO
2890 ENDDO
2891C-------------------------------------
2892 ELSE IF (ifunc == 10671) THEN
2893C SOUND SPEED
2894 ! /ANIM/ELEM/SSP
2895C-------------------------------------
2896 l = elbuf_tab(ng)%BUFLY(1)%L_SSP
2897 IF(l==0)THEN
2898 DO i=lft,llt
2899 evar(i) = zero
2900 ENDDO
2901 ELSE
2902 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
2903 DO i=lft,llt
2904 evar(i) = lbuf%SSP(i)
2905 ENDDO
2906 ENDIF
2907C-------------------------------------
2908 ELSE IF (ifunc == 10672) THEN
2909C ALE SCHLIEREN N/A
2910 !/ANIM/ELEM/SCHLIEREN
2911C-------------------------------------
2912 DO i=lft,llt
2913 evar(i) = zero
2914 ENDDO
2915C-----------------------------------------------
2916 ELSE IF (ifunc == 2156) THEN
2917C-----------------------------------------------
2918 IF (ity == 3) THEN
2919 DO i=lft,llt
2920 evar(i) = err_thk_sh4(nft+i)
2921 END DO
2922 ELSE
2923 DO i=lft,llt
2924 evar(i) = err_thk_sh3(nft+i)
2925 END DO
2926 ENDIF
2927C-------------------------------------
2928 ELSE IF (ifunc == 10676) THEN
2929C SPMD DOMAIN
2930C-------------------------------------
2931 DO i=lft,llt
2932 evar(i) = ispmd
2933 ENDDO
2934C-------------------------------------
2935 ELSEIF (ifunc == 10677) THEN ! /ANIM/SHELL/SIGEQ
2936 ! equivalent stress - other than VON MISES
2937C-------------------------------------
2938 IF (gbuf%G_SEQ > 0) THEN
2939C------------------
2940 ! Total number of integration points
2941 npgt = 0
2942 DO il=1,nlay
2943 bufly => elbuf_tab(ng)%BUFLY(il)
2944 npgt = npgt + bufly%NPTT*nptr*npts
2945 ENDDO
2946 ! Average equivalent stress on integration points
2947 DO i=lft,llt
2948 evar_tmp = zero
2949 DO il=1,nlay
2950 bufly => elbuf_tab(ng)%BUFLY(il)
2951 DO it=1,bufly%NPTT
2952 DO ir=1,nptr
2953 DO is=1,npts
2954 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
2955 evar_tmp = evar_tmp + lbuf%SEQ(i)/npgt
2956 ENDDO
2957 ENDDO
2958 ENDDO
2959 ENDDO
2960 evar(i) = evar_tmp
2961 ENDDO
2962C------------------
2963 ELSE ! VON MISES
2964 DO i=lft,llt
2965 s1 = gbuf%FOR(jj(1)+i)
2966 s2 = gbuf%FOR(jj(2)+i)
2967 s12= gbuf%FOR(jj(3)+i)
2968 vonm2 = s1*s1 + s2*s2 - s1*s2 + three*s12*s12
2969 evar(i) = sqrt(vonm2)
2970 ENDDO
2971 ENDIF ! IF (GBUF%G_SEQ > 0)
2972c------------------------------------
2973 ELSEIF (ifunc > 10677 .AND. ifunc < 10778 .AND.
2974 . (igtyp == 51 .OR. igtyp == 52).AND.
2975 . mlw /= 15 .AND. mlw /= 25 ) THEN
2976 ! EPSP/ILAY/UPPER
2977c------------------------------------
2978C available for IGTYP = 51 and 52 only (to be generalized)
2979 ilay = mod((ifunc - 10677), 100)
2980 IF (ilay == 0) ilay = 100
2981 IF (nlay > 1) THEN
2982 il = max(1,ilay)
2983 ELSE
2984 il = 1
2985 ENDIF
2986 bufly => elbuf_tab(ng)%BUFLY(il)
2987 nptt = bufly%NPTT
2988 ipt = max(1,nptt)
2989 IF (bufly%L_PLA > 0 .AND.
2990 . (il <= nlay .AND. ipt <= nptt))THEN
2991 IF (npg > 1) THEN
2992 DO i=lft,llt
2993 DO ir=1,nptr
2994 DO is=1,npts
2995 lbuf => bufly%LBUF(ir,is,ipt)
2996 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
2997 ENDDO
2998 ENDDO
2999 ENDDO
3000 ELSE ! (NPG == 1)
3001 lbuf => bufly%LBUF(1,1,ipt)
3002 DO i=lft,llt
3003 evar(i) = abs(lbuf%PLA(i))
3004 ENDDO
3005 ENDIF ! IF (NPG > 1) THEN
3006 ENDIF ! IF (BUFLY%L_PLA > 0) THEN
3007c------------------------------------
3008 ELSEIF (ifunc > 10777 .AND. ifunc < 10878 .AND.
3009 . (igtyp == 51 .OR. igtyp == 52) .AND.
3010 . mlw /= 15 .AND. mlw /= 25) THEN
3011 ! EPSP/ILAY/LOWER
3012c------------------------------------
3013C available for IGTYP = 51 and 52 only (to be generalized)
3014 ilay = mod((ifunc - 10777), 100)
3015 IF (ilay == 0) ilay = 100
3016 IF (nlay > 1) THEN
3017 il = max(1,ilay)
3018 ELSE
3019 il = 1
3020 ENDIF
3021 ipt = 1
3022 bufly => elbuf_tab(ng)%BUFLY(il)
3023 nptt = bufly%NPTT
3024 IF (bufly%L_PLA > 0 .AND.
3025 . (il <= nlay .AND. ipt <= nptt))THEN
3026 IF (npg > 1) THEN
3027 DO i=lft,llt
3028 DO ir=1,nptr
3029 DO is=1,npts
3030 lbuf => bufly%LBUF(ir,is,ipt)
3031 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
3032 ENDDO
3033 ENDDO
3034 ENDDO
3035 ELSE ! (NPG == 1)
3036 lbuf => bufly%LBUF(1,1,ipt)
3037 DO i=lft,llt
3038 evar(i) = abs(lbuf%PLA(i))
3039 ENDDO
3040 ENDIF ! IF (NPG > 1) THEN
3041 ENDIF ! IF (BUFLY%L_PLA > 0) THEN
3042c------------------------------------
3043 ELSEIF (ifunc > 10877 .AND. ifunc < 11888 .AND.
3044 . (igtyp == 51 .OR. igtyp == 52).AND.
3045 . mlw /= 15 .AND. mlw /= 25) THEN
3046 ! EPSP/ILAY/NIP
3047c------------------------------------
3048C
3049C available for IGTYP = 51 and 52 only (to be generalized)
3050C
3051 ius = ifunc - 10877
3052 il = int((ius - 1)/10)
3053 ipt = ius - 10*il
3054 IF (il <= nlay ) THEN
3055 bufly => elbuf_tab(ng)%BUFLY(il)
3056 nptt = bufly%NPTT
3057 IF (bufly%L_PLA > 0 .AND. ipt <= nptt) THEN
3058 IF (npg > 1) THEN
3059 DO i=lft,llt
3060 DO ir=1,nptr
3061 DO is=1,npts
3062 lbuf => bufly%LBUF(ir,is,ipt)
3063 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
3064 ENDDO
3065 ENDDO
3066 ENDDO
3067 ELSE ! (NPG == 1)
3068 lbuf => bufly%LBUF(1,1,ipt)
3069 DO i=lft,llt
3070 evar(i) = abs(lbuf%PLA(i))
3071 ENDDO
3072 ENDIF ! IF (NPG > 1) THEN
3073 ENDIF ! IF (BUFLY%L_PLA > 0) THEN
3074 ENDIF ! IF (IL <= NLAY )
3075c------------------------------------
3076 ELSEIF(ifunc == 11888)THEN
3077 !/ANIM/ELEM/BULK (QVIS)
3078c------------------------------------
3079 IF (gbuf%G_QVIS > 0) THEN
3080 DO i=lft,llt
3081 func(el2fa(nn3+nft+i)) = gbuf%QVIS(i)
3082 ENDDO
3083 ELSE
3084 DO i=lft,llt
3085 func(el2fa(nn3+nft+i)) = zero
3086 ENDDO
3087 ENDIF
3088c------------------------------------
3089 ELSEIF (ifunc == 11889) THEN
3090 IF (mlw /= 51 .AND. gbuf%G_TB > 0) THEN
3091 DO i=lft,llt
3092 func(el2fa(nn3+nft+i)) = -gbuf%TB(i)
3093 ENDDO
3094 ELSEIF (mlw == 51)THEN
3095 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
3096 ipos = 15
3097 itrimat = 4
3098 llt = iparg(2,ng)
3099 k = llt * ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1)
3100 DO i=lft,llt
3101 func(el2fa(nn3+nft+i)) = -mbuf%VAR(k+i)
3102 ENDDO
3103 ELSE
3104 DO i=lft,llt
3105 func(el2fa(nn3+nft+i)) = zero
3106 ENDDO
3107 ENDIF
3108C--------------------------------------------------
3109 ELSE IF (ifunc>11925 .AND. ifunc < 11925+mx_ply_anim+1)THEN
3110 !/ANIM/SHELL/IDPLY
3111c------------------------------------
3112 iply = ifunc - 11925
3113 IF (igtyp == 17 .OR. igtyp == 51) THEN
3114 IF (ply_anim( 3 * (iply - 1) + 2) == 1 )THEN
3115 DO j=1,nlay
3116 bufly => elbuf_tab(ng)%BUFLY(j)
3117 nptt = bufly%NPTT
3118 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
3119 IF (id_ply == ply_anim( 3 * (iply - 1) + 1) ) THEN
3120 DO i=lft,llt
3121 nb_plyoff = 0
3122 DO ir=1,nptr
3123 DO is=1,npts
3124 DO it=1,nptt
3125 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
3126 IF (lbuf%OFF(i) == zero) nb_plyoff = nb_plyoff + 1
3127 ENDDO
3128 ENDDO
3129 ENDDO
3130 IF ( nb_plyoff == nptr * npts * nptt ) THEN
3131 evar(i) = -one
3132 ELSE
3133 evar(i) = one
3134 ENDIF
3135 ENDDO
3136 ENDIF
3137 ENDDO
3138 ENDIF
3139 ELSEIF (igtyp == 52) THEN
3140 IF (ply_anim( 3 * (iply - 1) + 2) == 1 )THEN
3141 DO j=1,nlay
3142 bufly => elbuf_tab(ng)%BUFLY(j)
3143 nptt = bufly%NPTT
3144 id_ply = ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
3145 IF (id_ply == ply_anim( 3 * (iply - 1) + 1) ) THEN
3146 DO i=lft,llt
3147 nb_plyoff = 0
3148 DO ir=1,nptr
3149 DO is=1,npts
3150 DO it=1,nptt
3151 lbuf => elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,it)
3152 IF (lbuf%OFF(i) == zero) nb_plyoff = nb_plyoff + 1
3153 ENDDO
3154 ENDDO
3155 ENDDO
3156 IF ( nb_plyoff == nptr * npts * nptt ) THEN
3157 evar(i) = -one
3158 ELSE
3159 evar(i) = one
3160 ENDIF
3161 ENDDO
3162 ENDIF
3163 ENDDO
3164 ENDIF
3165 ENDIF
3166C--------------------------------------------------
3167 ELSE IF (ifunc> 11925+mx_ply_anim .AND. ifunc < 11925+(2*mx_ply_anim)+1)THEN
3168 !/ANIM/SHELL/IDPLY/PHI
3169c------------------------------------
3170 ivar = ifunc - (11925+mx_ply_anim)
3171 iply = ply_anim_phi( 3 * (ivar - 1) + 1)
3172 ipt = ply_anim_phi( 3 * (ivar - 1) + 3)
3173c
3174 DO j=1,nlay
3175 id_ply = 0
3176 IF (igtyp == 17 .OR. igtyp == 51) THEN
3177 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
3178 ELSEIF (igtyp == 52) THEN
3179 id_ply = ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
3180 ENDIF
3181c
3182 IF (id_ply == iply ) THEN
3183 bufly => elbuf_tab(ng)%BUFLY(j)
3184 IF (ity == 3) THEN
3185 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
3186 IF(ipt <= bufly%NPTT ) THEN
3187 lbuf_dir => elbuf_tab(ng)%BUFLY(j)%LBUF_DIR(ipt)
3188 ELSE
3189 lbuf_dir => elbuf_tab(ng)%BUFLY(j)%LBUF_DIR(1)
3190 ENDIF
3191 IF (mlw /= 0 .AND. mlw /= 13) THEN
3192 DO i=lft,llt
3193 n = i + nft
3194 x21 = x(1,ixc(3,n))-x(1,ixc(2,n))
3195 x32 = x(1,ixc(4,n))-x(1,ixc(3,n))
3196 x34 = x(1,ixc(4,n))-x(1,ixc(5,n))
3197 x41 = x(1,ixc(5,n))-x(1,ixc(2,n))
3198
3199 y21 = x(2,ixc(3,n))-x(2,ixc(2,n))
3200 y32 = x(2,ixc(4,n))-x(2,ixc(3,n))
3201 y34 = x(2,ixc(4,n))-x(2,ixc(5,n))
3202 y41 = x(2,ixc(5,n))-x(2,ixc(2,n))
3203
3204 z21 = x(3,ixc(3,n))-x(3,ixc(2,n))
3205 z32 = x(3,ixc(4,n))-x(3,ixc(3,n))
3206 z34 = x(3,ixc(4,n))-x(3,ixc(5,n))
3207 z41 = x(3,ixc(5,n))-x(3,ixc(2,n))
3208
3209 e1x = (x21+x34)
3210 e1y = (y21+y34)
3211 e1z = (z21+z34)
3212
3213 e2x = (x32+x41)
3214 e2y = (y32+y41)
3215 e2z = (z32+z41)
3216
3217 e3x = e1y*e2z-e1z*e2y
3218 e3y = e1z*e2x-e1x*e2z
3219 e3z = e1x*e2y-e1y*e2x
3220 IF (irep > 0) THEN
3221 rx = e1x
3222 ry = e1y
3223 rz = e1z
3224 sx = e2x
3225 sy = e2y
3226 sz = e2z
3227 ENDIF
3228 IF (ishfram == 0 .OR. igtyp == 16 ) THEN
3229C----- Symmetrical convected user-version 5 (DEFAULT)
3230 suma = e3x*e3x+e3y*e3y+e3z*e3z
3231 suma = one / max(sqrt(suma),em20)
3232 e3x = e3x * suma
3233 e3y = e3y * suma
3234 e3z = e3z * suma
3235C
3236 s1 = e1x*e1x+e1y*e1y+e1z*e1z
3237 s2 = e2x*e2x+e2y*e2y+e2z*e2z
3238 suma = sqrt(s1/s2)
3239 e1x = e1x + (e2y*e3z-e2z*e3y)*suma
3240 e1y = e1y + (e2z*e3x-e2x*e3z)*suma
3241 e1z = e1z + (e2x*e3y-e2y*e3x)*suma
3242C
3243 suma = e1x*e1x+e1y*e1y+e1z*e1z
3244 suma = one / max(sqrt(suma),em20)
3245 e1x = e1x * suma
3246 e1y = e1y * suma
3247 e1z = e1z * suma
3248C
3249 e2x = e3y * e1z - e3z * e1y
3250 e2y = e3z * e1x - e3x * e1z
3251 e2z = e3x * e1y - e3y * e1x
3252 ELSEIF (ishfram == 2) THEN
3253C----- Nononmetrical convected marker-version 4
3254 suma = e2x*e2x+e2y*e2y+e2z*e2z
3255 e1x = e1x*suma + e2y*e3z-e2z*e3y
3256 e1y = e1y*suma + e2z*e3x-e2x*e3z
3257 e1z = e1z*suma + e2x*e3y-e2y*e3x
3258 suma = e1x*e1x+e1y*e1y+e1z*e1z
3259 suma = one/max(sqrt(suma),em20)
3260 e1x = e1x*suma
3261 e1y = e1y*suma
3262 e1z = e1z*suma
3263C
3264 suma = e3x*e3x+e3y*e3y+e3z*e3z
3265 suma = one / max(sqrt(suma),em20)
3266 e3x = e3x * suma
3267 e3y = e3y * suma
3268 e3z = e3z * suma
3269C
3270 e2x = e3y*e1z-e3z*e1y
3271 e2y = e3z*e1x-e3x*e1z
3272 e2z = e3x*e1y-e3y*e1x
3273 suma = e2x*e2x+e2y*e2y+e2z*e2z
3274 suma = one/max(sqrt(suma),em20)
3275 e2x = e2x*suma
3276 e2y = e2y*suma
3277 e2z = e2z*suma
3278 ENDIF
3279 IF (irep >= 1) THEN
3280 aa = lbuf_dir%DIRA(i)
3281 bb = lbuf_dir%DIRA(i+nel)
3282 v1 = aa*rx + bb*sx
3283 v2 = aa*ry + bb*sy
3284 v3 = aa*rz + bb*sz
3285 vr = v1*e1x+ v2*e1y + v3*e1z
3286 vs = v1*e2x+ v2*e2y + v3*e2z
3287 suma=sqrt(vr*vr + vs*vs)
3288 dir1_1 = vr/suma
3289 dir1_2 = vs/suma
3290 ELSE
3291 dir1_1 = lbuf_dir%DIRA(i)
3292 dir1_2 = lbuf_dir%DIRA(i+nel)
3293 ENDIF
3294c
3295 phi =(hundred80/pi)*atan2(dir1_2,dir1_1)
3296 err = (abs(phi) - ninety)/ninety
3297 evar(i) = phi
3298 IF(abs(err) < em02) evar(i) = sign(ninety,phi)
3299 IF(abs(evar(i)) < one) evar(i) = zero
3300 ENDDO
3301 ENDIF ! MLW
3302 ELSEIF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52 ) THEN
3303 bufly => elbuf_tab(ng)%BUFLY(j)
3304 IF (mlw /= 0 .AND. mlw /= 13) THEN
3305 DO i=lft,llt
3306 n = i + nft
3307 x21 = x(1,ixc(3,n))-x(1,ixc(2,n))
3308 x32 = x(1,ixc(4,n))-x(1,ixc(3,n))
3309 x34 = x(1,ixc(4,n))-x(1,ixc(5,n))
3310 x41 = x(1,ixc(5,n))-x(1,ixc(2,n))
3311
3312 y21 = x(2,ixc(3,n))-x(2,ixc(2,n))
3313 y32 = x(2,ixc(4,n))-x(2,ixc(3,n))
3314 y34 = x(2,ixc(4,n))-x(2,ixc(5,n))
3315 y41 = x(2,ixc(5,n))-x(2,ixc(2,n))
3316
3317 z21 = x(3,ixc(3,n))-x(3,ixc(2,n))
3318 z32 = x(3,ixc(4,n))-x(3,ixc(3,n))
3319 z34 = x(3,ixc(4,n))-x(3,ixc(5,n))
3320 z41 = x(3,ixc(5,n))-x(3,ixc(2,n))
3321
3322 e1x = (x21+x34)
3323 e1y = (y21+y34)
3324 e1z = (z21+z34)
3325
3326 e2x = (x32+x41)
3327 e2y = (y32+y41)
3328 e2z = (z32+z41)
3329
3330 e3x = e1y*e2z-e1z*e2y
3331 e3y = e1z*e2x-e1x*e2z
3332 e3z = e1x*e2y-e1y*e2x
3333 IF (irep > 0) THEN
3334 rx = e1x
3335 ry = e1y
3336 rz = e1z
3337 sx = e2x
3338 sy = e2y
3339 sz = e2z
3340 ENDIF
3341 IF (ishfram == 0 .OR. igtyp == 16 ) THEN
3342C----- Symmetrical convected user-version 5 (DEFAULT)
3343 suma = e3x*e3x+e3y*e3y+e3z*e3z
3344 suma = one / max(sqrt(suma),em20)
3345 e3x = e3x * suma
3346 e3y = e3y * suma
3347 e3z = e3z * suma
3348C
3349 s1 = e1x*e1x+e1y*e1y+e1z*e1z
3350 s2 = e2x*e2x+e2y*e2y+e2z*e2z
3351 suma = sqrt(s1/s2)
3352 e1x = e1x + (e2y*e3z-e2z*e3y)*suma
3353 e1y = e1y + (e2z*e3x-e2x*e3z)*suma
3354 e1z = e1z + (e2x*e3y-e2y*e3x)*suma
3355C
3356 suma = e1x*e1x+e1y*e1y+e1z*e1z
3357 suma = one / max(sqrt(suma),em20)
3358 e1x = e1x * suma
3359 e1y = e1y * suma
3360 e1z = e1z * suma
3361C
3362 e2x = e3y * e1z - e3z * e1y
3363 e2y = e3z * e1x - e3x * e1z
3364 e2z = e3x * e1y - e3y * e1x
3365 ELSEIF (ishfram == 2) THEN
3366C----- Nononmetrical convected marker-version 4
3367 suma = e2x*e2x+e2y*e2y+e2z*e2z
3368 e1x = e1x*suma + e2y*e3z-e2z*e3y
3369 e1y = e1y*suma + e2z*e3x-e2x*e3z
3370 e1z = e1z*suma + e2x*e3y-e2y*e3x
3371 suma = e1x*e1x+e1y*e1y+e1z*e1z
3372 suma = one/max(sqrt(suma),em20)
3373 e1x = e1x*suma
3374 e1y = e1y*suma
3375 e1z = e1z*suma
3376C
3377 suma = e3x*e3x+e3y*e3y+e3z*e3z
3378 suma = one / max(sqrt(suma),em20)
3379 e3x = e3x * suma
3380 e3y = e3y * suma
3381 e3z = e3z * suma
3382C
3383 e2x = e3y*e1z-e3z*e1y
3384 e2y = e3z*e1x-e3x*e1z
3385 e2z = e3x*e1y-e3y*e1x
3386 suma = e2x*e2x+e2y*e2y+e2z*e2z
3387 suma = one/max(sqrt(suma),em20)
3388 e2x = e2x*suma
3389 e2y = e2y*suma
3390 e2z = e2z*suma
3391 ENDIF
3392 IF (irep >= 1) THEN
3393 aa = bufly%DIRA(i)
3394 bb = bufly%DIRA(i+nel)
3395 v1 = aa*rx + bb*sx
3396 v2 = aa*ry + bb*sy
3397 v3 = aa*rz + bb*sz
3398 vr = v1*e1x+ v2*e1y + v3*e1z
3399 vs = v1*e2x+ v2*e2y + v3*e2z
3400 suma=sqrt(vr*vr + vs*vs)
3401 dir1_1 = vr/suma
3402 dir1_2 = vs/suma
3403 ELSE
3404 dir1_1 = bufly%DIRA(i)
3405 dir1_2 = bufly%DIRA(i+nel)
3406 ENDIF
3407c
3408 phi =(hundred80/pi)*atan2(dir1_2,dir1_1)
3409 err = (abs(phi) - ninety)/ninety
3410 evar(i) = phi
3411 IF(abs(err) < em02) evar(i) = sign(ninety,phi)
3412 IF(abs(evar(i)) < one) evar(i) = zero
3413 ENDDO
3414 ENDIF ! MLW
3415 ENDIF ! IGTYP + IDRAPE
3416
3417 ELSEIF (ity == 7) THEN
3418 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
3419 IF(ipt <= bufly%NPTT ) THEN
3420 lbuf_dir => elbuf_tab(ng)%BUFLY(j)%LBUF_DIR(ipt)
3421 ELSE
3422 lbuf_dir => elbuf_tab(ng)%BUFLY(j)%LBUF_DIR(1)
3423 ENDIF
3424 IF (mlw /= 0 .AND. mlw /= 13) THEN
3425 DO i=lft,llt
3426 n = i + nft
3427 x21 = x(1,ixtg(3,n))-x(1,ixtg(2,n))
3428 x31 = x(1,ixtg(4,n))-x(1,ixtg(2,n))
3429 x32 = x(1,ixtg(4,n))-x(1,ixtg(3,n))
3430
3431 y21 = x(2,ixtg(3,n))-x(2,ixtg(2,n))
3432 y31 = x(2,ixtg(4,n))-x(2,ixtg(2,n))
3433 y32 = x(2,ixtg(4,n))-x(2,ixtg(3,n))
3434
3435 z21 = x(3,ixtg(3,n))-x(3,ixtg(2,n))
3436 z31 = x(3,ixtg(4,n))-x(3,ixtg(2,n))
3437 z32 = x(3,ixtg(4,n))-x(3,ixtg(3,n))
3438 IF (irep > 0) THEN
3439 e11 = x21
3440 e12 = y21
3441 e13 = z21
3442 e21 = x31
3443 e22 = y31
3444 e23 = z31
3445 ENDIF
3446 IF(ifram_old ==0 ) THEN
3447 CALL clsconv3(x21,y21,z21,x31,y31,z31,
3448 + e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z)
3449 ELSE
3450 e1x= x21
3451 e1y= y21
3452 e1z= z21
3453 x2l = sqrt(e1x*e1x+e1y*e1y+e1z*e1z)
3454 e1x=e1x/x2l
3455 e1y=e1y/x2l
3456 e1z=e1z/x2l
3457C
3458 e3x=y31*z32-z31*y32
3459 e3y=z31*x32-x31*z32
3460 e3z=x31*y32-y31*x32
3461 sum_ = sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
3462 e3x=e3x/sum_
3463 e3y=e3y/sum_
3464 e3z=e3z/sum_
3465 area = half * sum_
3466 e2x=e3y*e1z-e3z*e1y
3467 e2y=e3z*e1x-e3x*e1z
3468 e2z=e3x*e1y-e3y*e1x
3469 sum_ = sqrt(e2x*e2x+e2y*e2y+e2z*e2z)
3470 e2x=e2x/sum_
3471 e2y=e2y/sum_
3472 e2z=e2z/sum_
3473 END IF
3474 IF (irep >= 1) THEN
3475 aa = lbuf_dir%DIRA(i)
3476 bb = lbuf_dir%DIRA(i+nel)
3477 v1 = aa*e11 + bb*e21
3478 v2 = aa*e12 + bb*e22
3479 v3 = aa*e13 + bb*e23
3480 vr = v1*e1x + v2*e1y + v3*e1z
3481 vs = v1*e2x + v2*e2y + v3*e2z
3482 suma=sqrt(vr*vr + vs*vs)
3483 dir1_1 = vr/suma
3484 dir1_2 = vs/suma
3485 ELSE
3486 dir1_1 = lbuf_dir%DIRA(i)
3487 dir1_2 = lbuf_dir%DIRA(i+nel)
3488 ENDIF
3489 phi =(hundred80/pi)*atan2(dir1_2,dir1_1)
3490 err = (abs(phi) - ninety)/ninety
3491 evar(i) = phi
3492 IF(abs(err) < em02) evar(i) = sign(ninety,phi)
3493 IF(abs(evar(i)) < one) evar(i) = zero
3494 ENDDO
3495 ENDIF
3496 ELSEIF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52 ) THEN
3497 bufly => elbuf_tab(ng)%BUFLY(j)
3498 IF (mlw /= 0 .AND. mlw /= 13) THEN
3499 DO i=lft,llt
3500 n = i + nft
3501 x21 = x(1,ixtg(3,n))-x(1,ixtg(2,n))
3502 x31 = x(1,ixtg(4,n))-x(1,ixtg(2,n))
3503 x32 = x(1,ixtg(4,n))-x(1,ixtg(3,n))
3504
3505 y21 = x(2,ixtg(3,n))-x(2,ixtg(2,n))
3506 y31 = x(2,ixtg(4,n))-x(2,ixtg(2,n))
3507 y32 = x(2,ixtg(4,n))-x(2,ixtg(3,n))
3508
3509 z21 = x(3,ixtg(3,n))-x(3,ixtg(2,n))
3510 z31 = x(3,ixtg(4,n))-x(3,ixtg(2,n))
3511 z32 = x(3,ixtg(4,n))-x(3,ixtg(3,n))
3512 IF (irep > 0) THEN
3513 e11 = x21
3514 e12 = y21
3515 e13 = z21
3516 e21 = x31
3517 e22 = y31
3518 e23 = z31
3519 ENDIF
3520 IF(ifram_old ==0 ) THEN
3521 CALL clsconv3(x21,y21,z21,x31,y31,z31,
3522 + e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z)
3523 ELSE
3524 e1x= x21
3525 e1y= y21
3526 e1z= z21
3527 x2l = sqrt(e1x*e1x+e1y*e1y+e1z*e1z)
3528 e1x=e1x/x2l
3529 e1y=e1y/x2l
3530 e1z=e1z/x2l
3531C
3532 e3x=y31*z32-z31*y32
3533 e3y=z31*x32-x31*z32
3534 e3z=x31*y32-y31*x32
3535 sum_ = sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
3536 e3x=e3x/sum_
3537 e3y=e3y/sum_
3538 e3z=e3z/sum_
3539 area = half * sum_
3540 e2x=e3y*e1z-e3z*e1y
3541 e2y=e3z*e1x-e3x*e1z
3542 e2z=e3x*e1y-e3y*e1x
3543 sum_ = sqrt(e2x*e2x+e2y*e2y+e2z*e2z)
3544 e2x=e2x/sum_
3545 e2y=e2y/sum_
3546 e2z=e2z/sum_
3547 END IF
3548 IF (irep >= 1) THEN
3549 aa = bufly%DIRA(i)
3550 bb = bufly%DIRA(i+nel)
3551 v1 = aa*e11 + bb*e21
3552 v2 = aa*e12 + bb*e22
3553 v3 = aa*e13 + bb*e23
3554 vr = v1*e1x + v2*e1y + v3*e1z
3555 vs = v1*e2x + v2*e2y + v3*e2z
3556 suma=sqrt(vr*vr + vs*vs)
3557 dir1_1 = vr/suma
3558 dir1_2 = vs/suma
3559 ELSE
3560 dir1_1 = bufly%DIRA(i)
3561 dir1_2 = bufly%DIRA(i+nel)
3562 ENDIF
3563 phi =(hundred80/pi)*atan2(dir1_2,dir1_1)
3564 err = (abs(phi) - ninety)/ninety
3565 evar(i) = phi
3566 IF(abs(err) < em02) evar(i) = sign(ninety,phi)
3567 IF(abs(evar(i)) < one) evar(i) = zero
3568 ENDDO
3569 ENDIF
3570 ENDIF
3571 ENDIF
3572c
3573 ENDIF
3574 ENDDO
3575c------------------------------------
3576 ELSE IF (ifunc> 11925+(2*mx_ply_anim) .AND. ifunc < 11925+(3*mx_ply_anim)+1)THEN
3577 !/ANIM/SHELL/IDPLY/EPSP
3578c------------------------------------
3579 iply = ifunc - (11925+ 2*mx_ply_anim)
3580 ipt = ply_anim_epsp( 3 * (iply - 1) + 3)
3581c
3582 DO j=1,nlay
3583 id_ply = 0
3584 IF (igtyp == 17 .OR. igtyp == 51) THEN
3585 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
3586 ELSEIF (igtyp == 52) THEN
3587 id_ply = ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
3588 ENDIF
3589c
3590 IF (id_ply == ply_anim_epsp( 3 * (iply - 1) + 1) ) THEN
3591 bufly => elbuf_tab(ng)%BUFLY(j)
3592 IF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52 ) THEN
3593 nptt = bufly%NPTT
3594 IF( ipt <= nptt) THEN
3595 IF( npg > 1 ) THEN
3596 DO i=lft,llt
3597 DO ir=1,nptr
3598 DO is=1,npts
3599 evar(i) = evar(i) + abs(bufly%LBUF(ir,is,ipt)%PLA(i))/npg
3600 ENDDO
3601 ENDDO
3602 ENDDO
3603 ELSE
3604 DO i=lft,llt
3605 evar(i) = abs(bufly%LBUF(1,1,ipt)%PLA(i))
3606 ENDDO
3607 ENDIF
3608c
3609 ELSE
3610 DO i=lft,llt
3611 evar(i) = zero
3612 ENDDO
3613 ENDIF
3614 ELSE
3615 DO i=lft,llt
3616 evar(i) = zero
3617 ENDDO
3618 ENDIF
3619 ENDIF
3620 ENDDO
3621c------------------------------------
3622 ELSE IF (ifunc> 11925+(3*mx_ply_anim) .AND. ifunc < 11925+(4*mx_ply_anim)+1)THEN
3623 !/ANIM/SHELL/IDPLY/DAMA
3624c------------------------------------
3625 iply = ifunc - (11925+ 3*mx_ply_anim)
3626 ipt = ply_anim_dama( 3 * (iply - 1) + 3)
3627c
3628 IF(ifailure > 0) THEN
3629 DO j=1,nlay
3630 nptt = elbuf_tab(ng)%BUFLY(j)%NPTT
3631 id_ply = 0
3632 IF (igtyp == 17 .OR. igtyp == 51) THEN
3633 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
3634 ELSEIF (igtyp == 52) THEN
3635 id_ply=ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
3636 ENDIF
3637 IF (id_ply == ply_anim_dama( 3 *(iply - 1) + 1) )THEN
3638 IF (ipt <= nptt) THEN
3639 DO i=lft,llt
3640 DO ir = 1, nptr
3641 DO is = 1, npts
3642 fbuf => elbuf_tab(ng)%BUFLY(j)%FAIL(ir,is,ipt)
3643 DO ifail = 1, elbuf_tab(ng)%BUFLY(j)%NFAIL
3644 evar(i) = max(evar(i), fbuf%FLOC(ifail)%DAMMX(i))
3645 ENDDO
3646 ENDDO
3647 ENDDO
3648 ENDDO ! I=LFT,LLT
3649 ENDIF
3650 ENDIF
3651 ENDDO
3652 ENDIF
3653c
3654 IF(mlw == 25 .AND. (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52)) THEN
3655 IF(ity == 3)THEN
3656 DO i=lft,llt
3657 mat(i)=ixc(1,nft+i)
3658 pid(i)=ixc(6,nft+i)
3659 END DO
3660 ELSE
3661 DO i=lft,llt
3662 mat(i)=ixtg(1,nft+i)
3663 pid(i)=ixtg(5,nft+i)
3664 END DO
3665 END IF
3666c
3667 ipmat = 2 + nlay
3668 DO n=1,nlay
3669 iadr = (n-1)*nel
3670 DO i=lft,llt
3671 j = iadr+i
3672 matly(j) = stack%IGEO(ipmat+n,isubstack)
3673 END DO
3674 END DO
3675c
3676 DO j=1,nlay
3677 id_ply = 0
3678 IF (igtyp == 17 .OR. igtyp == 51) THEN
3679 id_ply = igeo(1,stack%IGEO(2+j,isubstack))
3680 ELSEIF (igtyp == 52) THEN
3681 id_ply=ply_info(1,stack%IGEO(2+j,isubstack)-numstack)
3682 ENDIF
3683c
3684 IF (id_ply == ply_anim_dama( 3 *(iply - 1) + 1) )THEN
3685 bufly => elbuf_tab(ng)%BUFLY(j)
3686 DO i=lft,llt
3687 nptt = elbuf_tab(ng)%BUFLY(j)%NPTT
3688 IF (ipt <= nptt) THEN
3689 iadr = (j - 1)*nel
3690 vly(1:5) = zero
3691 vg(1:5)= zero
3692 DO ir=1,nptr
3693 DO is=1,npts
3694 lbuf=> elbuf_tab(ng)%BUFLY(j)%LBUF(ir,is,ipt)
3695 dmax(i) = one/pm(64,matly(iadr + i))
3696 wpmax(i)= one/pm(41,matly(iadr + i))
3697 epst1(i)= pm(60,matly(iadr + i))
3698 epst2(i)= pm(61,matly(iadr + i))
3699 epsf1(i)= one/pm(98,matly(iadr + i))
3700 epsf2(i)= one/pm(99,matly(iadr + i))
3701C
3702 vg(1) = max(vg(1),lbuf%DMG(jj(2)+i)*dmax(i))
3703 vg(2) = max(vg(2),lbuf%DMG(jj(3)+i)*dmax(i))
3704 vg(3)= max(vg(3),abs(lbuf%PLA(i))*wpmax(i))
3705 IF(lbuf%CRAK(jj(1)+i) > zero) vg(4)= max(vg(4),
3706 . (lbuf%CRAK(jj(1)+i)+epst1(i))*epsf1(i))
3707 IF(lbuf%CRAK(jj(2)+i) > zero )vg(5) = max(vg(5),
3708 . (lbuf%CRAK(jj(2)+i)+epst2(i))*epsf2(i))
3709 ENDDO
3710 ENDDO
3711 ENDIF
3712 vly(1) =vg(1)
3713 vly(2) =vg(2)
3714 vly(3) =vg(3)
3715 vly(4) =vg(4)
3716 vly(5) =vg(5)
3717C
3718 evar(i) = max(evar(i),vly(1),vly(2),vly(3),vly(4),vly(5))
3719 ENDDO ! I=LFT,JLT
3720
3721 ENDIF
3722 ENDDO
3723 ENDIF
3724C-------------------
3725 ELSEIF (ifunc > 11925+4*mx_ply_anim .and.
3726 . ifunc < 11925+4*mx_ply_anim + 4) THEN ! FLD Damage factor
3727C-------------------
3728 idx = 11925+4*mx_ply_anim
3729 IF (ifunc == idx+1) THEN ! /ANIM/SHELL/FLDF/UPPER
3730 IF (nlay > 1) THEN
3731 il = nlay
3732 ipt = 1
3733 ELSE
3734 il = 1
3735 ipt = nptt
3736 ENDIF
3737 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
3738 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
3739 DO is=1,npts
3740 DO ir=1,nptr
3741 DO it=1,nptt
3742 ipt = it
3743 IF (nlay == 1) ipt = nptt
3744 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
3745 DO ifail=1,nfail
3746 IF (fbuf%FLOC(ifail)%ILAWF == 7) THEN ! check /FLD model
3747 DO i=lft,llt
3748 evar(i) = max(evar(i),fbuf%FLOC(ifail)%DAM(i))
3749 ENDDO
3750 ENDIF
3751 ENDDO
3752 ENDDO
3753 ENDDO
3754 ENDDO
3755c
3756 ELSEIF (ifunc == idx+2) THEN ! /ANIM/SHELL/FLDF/LOWER
3757 il = 1
3758 ipt = 1
3759 bufly => elbuf_tab(ng)%BUFLY(il)
3760 nptt = bufly%NPTT
3761 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
3762 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
3763 DO is=1,npts
3764 DO ir=1,nptr
3765 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
3766 DO ifail=1,nfail
3767 IF (fbuf%FLOC(ifail)%ILAWF == 7) THEN ! check /FLD model
3768 DO i=lft,llt
3769 evar(i) = max(evar(i),fbuf%FLOC(ifail)%DAM(i))
3770 ENDDO
3771 ENDIF
3772 ENDDO
3773 ENDDO
3774 ENDDO
3775c
3776 ELSEIF (ifunc == idx+3) THEN ! /ANIM/SHELL/FLDF/MEMB
3777 il = nlay / 2 + 1
3778 bufly => elbuf_tab(ng)%BUFLY(il)
3779 nptt = bufly%NPTT
3780 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
3781 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
3782 ipt = nptt/2 + 1
3783 DO is=1,npts
3784 DO ir=1,nptr
3785 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
3786 DO ifail=1,nfail
3787 IF (fbuf%FLOC(ifail)%ILAWF == 7) THEN ! check /FLD model
3788 DO i=lft,llt
3789 evar(i) = max(evar(i),fbuf%FLOC(ifail)%DAM(i))
3790 ENDDO
3791 ENDIF
3792 ENDDO
3793 ENDDO
3794 ENDDO
3795 ENDIF
3796C-------------------
3797 ELSEIF (ifunc > 11925+4*mx_ply_anim + 3.and.
3798 . ifunc < 11925+4*mx_ply_anim + 7) THEN ! FLD zone index
3799C-------------------
3800 idx = 11925+4*mx_ply_anim + 3
3801 IF (ifunc == idx+1) THEN ! /ANIM/SHELL/FLDZ/UPPER
3802 IF (nlay > 1) THEN
3803 il = nlay
3804 ipt = 1
3805 ELSE
3806 il = 1
3807 ipt = nptt
3808 ENDIF
3809 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
3810 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
3811 DO is=1,npts
3812 DO ir=1,nptr
3813 DO it=1,nptt
3814 ipt = it
3815 IF (nlay == 1) ipt = nptt
3816 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
3817 DO ifail=1,nfail
3818 IF (fbuf%FLOC(ifail)%ILAWF == 7) THEN ! check /FLD model
3819 DO i=lft,llt
3820 rindx = fbuf%FLOC(ifail)%INDX(i)
3821 evar(i) = max(evar(i),rindx)
3822 ENDDO
3823 ENDIF
3824 ENDDO
3825 ENDDO
3826 ENDDO
3827 ENDDO
3828c
3829 ELSEIF (ifunc == idx+2) THEN ! /ANIM/SHELL/FLDZ/LOWER
3830 il = 1
3831 ipt = 1
3832 bufly => elbuf_tab(ng)%BUFLY(il)
3833 nptt = bufly%NPTT
3834 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
3835 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
3836 DO is=1,npts
3837 DO ir=1,nptr
3838 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
3839 DO ifail=1,nfail
3840 IF (fbuf%FLOC(ifail)%ILAWF == 7) THEN ! check /FLD model
3841 DO i=lft,llt
3842 rindx = fbuf%FLOC(ifail)%INDX(i)
3843 evar(i) = max(evar(i),rindx)
3844 ENDDO
3845 ENDIF
3846 ENDDO
3847 ENDDO
3848 ENDDO
3849c
3850 ELSEIF (ifunc == idx+3) THEN ! /ANIM/SHELL/FLDZ/MEMB
3851 il = nlay / 2 + 1
3852 bufly => elbuf_tab(ng)%BUFLY(il)
3853 nptt = bufly%NPTT
3854 nfail = elbuf_tab(ng)%BUFLY(il)%NFAIL
3855 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
3856 ipt = nptt/2 + 1
3857 DO is=1,npts
3858 DO ir=1,nptr
3859 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,ipt)
3860 DO ifail=1,nfail
3861 IF (fbuf%FLOC(ifail)%ILAWF == 7) THEN ! check /FLD model
3862 DO i=lft,llt
3863 rindx = fbuf%FLOC(ifail)%INDX(i)
3864 evar(i) = max(evar(i),rindx)
3865 ENDDO
3866 ENDIF
3867 ENDDO
3868 ENDDO
3869 ENDDO
3870 ENDIF
3871C------------------------------------------------------------------------------
3872C------------------- Damage Output : ALL NPTT through ALL Layers --- PID_51, 52
3873C------------------------------------------------------------------------------
3874 ELSEIF (ifunc > 11925+4*mx_ply_anim+6 .AND. ifunc < 11925+4*mx_ply_anim+107
3875 . .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
3876C-------------------
3877 ! DAMA/ILAY/UPPER
3878C-------------------
3879 idx = 11925+4*mx_ply_anim+6
3880 ilay = mod((ifunc - idx), 100)
3881 IF (ilay == 0) ilay = 100
3882 IF (nlay > 1) THEN
3883 il = max(1,ilay)
3884 ELSE
3885 il = 1
3886 ENDIF
3887 bufly => elbuf_tab(ng)%BUFLY(il)
3888 nptt = bufly%NPTT
3889 it = max(1,nptt)
3890C
3891 DO i=lft,llt
3892 evar(i) = zero
3893 ENDDO
3894C
3895 IF (ifailure > 0) THEN
3896 IF (il <= nlay .AND. it <= nptt) THEN
3897 DO i = lft,llt
3898 DO ir = 1, nptr
3899 DO is = 1, npts
3900 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,it)
3901 DO ifail = 1, elbuf_tab(ng)%BUFLY(il)%NFAIL
3902 evar(i) = max(evar(i), fbuf%FLOC(ifail)%DAMMX(i))
3903 ENDDO
3904 ENDDO
3905 ENDDO
3906 ENDDO ! I = LFT,LLT
3907 ENDIF ! IF (IL <= NLAY .AND. IT <= NPTT)
3908 ENDIF ! IF (IFAILURE > 0)
3909C---
3910C for outp of dam inside law25
3911C---
3912 IF (mlw == 25 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
3913 IF (ity == 3) THEN
3914 DO i=lft,llt
3915 mat(i)=ixc(1,nft+i)
3916 pid(i)=ixc(6,nft+i)
3917 ENDDO
3918 ELSE
3919 DO i=lft,llt
3920 mat(i)=ixtg(1,nft+i)
3921 pid(i)=ixtg(5,nft+i)
3922 ENDDO
3923 ENDIF
3924C
3925 ipmat = 2 + nlay
3926 DO n=1,nlay
3927 iadr = (n-1)*nel
3928 DO i=lft,llt
3929 j = iadr+i
3930 matly(j) = stack%IGEO(ipmat+n,isubstack)
3931 ENDDO
3932 ENDDO
3933C
3934 IF (il <= nlay .AND. it <= nptt) THEN
3935 DO i=lft,llt
3936 iadr = (il - 1)*nel
3937 j = iadr + i
3938 vg(1:5)= zero
3939 DO ir=1,nptr
3940 DO is=1,npts
3941 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
3942 dmax(i) = one/pm(64,matly(j))
3943 wpmax(i)= one/pm(41,matly(j))
3944 epst1(i)= pm(60,matly(j))
3945 epst2(i)= pm(61,matly(j))
3946 epsf1(i)= one/pm(98,matly(j))
3947 epsf2(i)= one/pm(99,matly(j))
3948C
3949 vg(1) = max(vg(1),lbuf%DMG(jj(2)+i)*dmax(i))
3950 vg(2) = max(vg(2),lbuf%DMG(jj(3)+i)*dmax(i))
3951 vg(3) = max(vg(3),abs(lbuf%PLA(i))*wpmax(i))
3952 IF (lbuf%CRAK(jj(1)+i) > zero) vg(4) = max(vg(4),
3953 . (lbuf%CRAK(jj(1)+i)+epst1(i))*epsf1(i))
3954 IF (lbuf%CRAK(jj(2)+i) > zero) vg(5) = max(vg(5),
3955 . (lbuf%CRAK(jj(2)+i)+epst2(i))*epsf2(i))
3956 ENDDO
3957 ENDDO
3958 evar(i) = max(evar(i),vg(1),vg(2),vg(3),vg(4),vg(5))
3959 ENDDO ! I=LFT,JLT
3960 ENDIF ! IF (IL <= NLAY .AND. IT <= NPTT)
3961 ENDIF ! IF (MLW == 25 .AND. (IGTYP == 51 .OR. IGTYP == 52 ))
3962C-------------------
3963 ELSEIF (ifunc > 11925+4*mx_ply_anim+106 .AND. ifunc < 11925+4*mx_ply_anim+207
3964 . .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
3965C-------------------
3966 ! DAMA/ILAY/LOWER
3967C-------------------
3968 idx = 11925+4*mx_ply_anim+106
3969 ilay = mod((ifunc - idx), 100)
3970 IF (ilay == 0) ilay = 100
3971 IF (nlay > 1) THEN
3972 il = max(1,ilay)
3973 ELSE
3974 il = 1
3975 ENDIF
3976 it = 1
3977 bufly => elbuf_tab(ng)%BUFLY(il)
3978 nptt = bufly%NPTT
3979C
3980 DO i=lft,llt
3981 evar(i) = zero
3982 ENDDO
3983C
3984 IF (ifailure > 0) THEN
3985 IF (il <= nlay .AND. it <= nptt) THEN
3986 DO i = lft,llt
3987 DO ir = 1, nptr
3988 DO is = 1, npts
3989 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,it)
3990 DO ifail = 1, elbuf_tab(ng)%BUFLY(il)%NFAIL
3991 evar(i) = max(evar(i), fbuf%FLOC(ifail)%DAMMX(i))
3992 ENDDO
3993 ENDDO
3994 ENDDO
3995 ENDDO ! I = LFT,LLT
3996 ENDIF ! IF (IL <= NLAY .AND. IT <= NPTT)
3997 ENDIF ! IF (IFAILURE > 0) THEN
3998C---
3999C for outp of dam inside law25
4000C---
4001 IF (mlw == 25 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
4002 IF (ity == 3) THEN
4003 DO i=lft,llt
4004 mat(i)=ixc(1,nft+i)
4005 pid(i)=ixc(6,nft+i)
4006 ENDDO
4007 ELSE
4008 DO i=lft,llt
4009 mat(i)=ixtg(1,nft+i)
4010 pid(i)=ixtg(5,nft+i)
4011 ENDDO
4012 ENDIF
4013C
4014 ipmat = 2 + nlay
4015 DO n=1,nlay
4016 iadr = (n-1)*nel
4017 DO i=lft,llt
4018 j = iadr+i
4019 matly(j) = stack%IGEO(ipmat+n,isubstack)
4020 ENDDO
4021 ENDDO
4022C
4023 IF (il <= nlay .AND. it <= nptt) THEN
4024 DO i=lft,llt
4025 iadr = (il - 1)*nel
4026 j = iadr + i
4027 vg(1:5)= zero
4028 DO ir=1,nptr
4029 DO is=1,npts
4030 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
4031 dmax(i) = one/pm(64,matly(j))
4032 wpmax(i)= one/pm(41,matly(j))
4033 epst1(i)= pm(60,matly(j))
4034 epst2(i)= pm(61,matly(j))
4035 epsf1(i)= one/pm(98,matly(j))
4036 epsf2(i)= one/pm(99,matly(j))
4037C
4038 vg(1) = max(vg(1),lbuf%DMG(jj(2)+i)*dmax(i))
4039 vg(2) = max(vg(2),lbuf%DMG(jj(3)+i)*dmax(i))
4040 vg(3) = max(vg(3),abs(lbuf%PLA(i))*wpmax(i))
4041 IF (lbuf%CRAK(jj(1)+i) > zero) vg(4) = max(vg(4),
4042 . (lbuf%CRAK(jj(1)+i)+epst1(i))*epsf1(i))
4043 IF (lbuf%CRAK(jj(2)+i) > zero) vg(5) = max(vg(5),
4044 . (lbuf%CRAK(jj(2)+i)+epst2(i))*epsf2(i))
4045 ENDDO
4046 ENDDO
4047 evar(i) = max(evar(i),vg(1),vg(2),vg(3),vg(4),vg(5))
4048 ENDDO ! I=LFT,JLT
4049 ENDIF ! IF (IL <= NLAY .AND. IT <= NPTT)
4050 ENDIF ! IF (MLW == 25 .AND. (IGTYP == 51 .OR. IGTYP == 52 ))
4051C-------------------
4052 ELSEIF (ifunc > 11925+4*mx_ply_anim+206 .AND. ifunc < 11925+4*mx_ply_anim+307
4053 . .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
4054C-------------------
4055 ! DAMA/ILAY/MEMB
4056C-------------------
4057 idx = 11925+4*mx_ply_anim+206
4058 ilay = mod((ifunc - idx), 100)
4059 IF (nlay > 1) THEN
4060 il = max(1,ilay)
4061 ELSE
4062 il = 1
4063 ENDIF
4064 bufly => elbuf_tab(ng)%BUFLY(il)
4065 nptt = bufly%NPTT
4066 it = nptt/2 + 1 ! MEMB of layer ILAY
4067C
4068 DO i=lft,llt
4069 evar(i) = zero
4070 ENDDO
4071C
4072 IF (ifailure > 0) THEN
4073 IF (il <= nlay .AND. it <= nptt) THEN
4074 DO i = lft,llt
4075 DO ir = 1, nptr
4076 DO is = 1, npts
4077 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,it)
4078 DO ifail = 1, elbuf_tab(ng)%BUFLY(il)%NFAIL
4079 evar(i) = max(evar(i), fbuf%FLOC(ifail)%DAMMX(i))
4080 ENDDO
4081 ENDDO
4082 ENDDO
4083 ENDDO ! I = LFT,LLT
4084 ENDIF ! IF (IL <= NLAY .AND. IT <= NPTT)
4085 ENDIF ! IF (IFAILURE > 0) THEN
4086C---
4087C for outp of dam inside law25
4088C---
4089 IF (mlw == 25 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
4090 IF (ity == 3) THEN
4091 DO i=lft,llt
4092 mat(i)=ixc(1,nft+i)
4093 pid(i)=ixc(6,nft+i)
4094 ENDDO
4095 ELSE
4096 DO i=lft,llt
4097 mat(i)=ixtg(1,nft+i)
4098 pid(i)=ixtg(5,nft+i)
4099 ENDDO
4100 ENDIF
4101C
4102 ipmat = 2 + nlay
4103 DO n=1,nlay
4104 iadr = (n-1)*nel
4105 DO i=lft,llt
4106 j = iadr+i
4107 matly(j) = stack%IGEO(ipmat+n,isubstack)
4108 ENDDO
4109 ENDDO
4110C
4111 IF (il <= nlay .AND. it <= nptt) THEN
4112 DO i=lft,llt
4113 iadr = (il - 1)*nel
4114 j = iadr + i
4115 vg(1:5)= zero
4116 DO ir=1,nptr
4117 DO is=1,npts
4118 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
4119 dmax(i) = one/pm(64,matly(j))
4120 wpmax(i)= one/pm(41,matly(j))
4121 epst1(i)= pm(60,matly(j))
4122 epst2(i)= pm(61,matly(j))
4123 epsf1(i)= one/pm(98,matly(j))
4124 epsf2(i)= one/pm(99,matly(j))
4125C
4126 vg(1) = max(vg(1),lbuf%DMG(jj(2)+i)*dmax(i))
4127 vg(2) = max(vg(2),lbuf%DMG(jj(3)+i)*dmax(i))
4128 vg(3) = max(vg(3),abs(lbuf%PLA(i))*wpmax(i))
4129 IF (lbuf%CRAK(jj(1)+i) > zero) vg(4) = max(vg(4),
4130 . (lbuf%CRAK(jj(1)+i)+epst1(i))*epsf1(i))
4131 IF (lbuf%CRAK(jj(2)+i) > zero) vg(5) = max(vg(5),
4132 . (lbuf%CRAK(jj(2)+i)+epst2(i))*epsf2(i))
4133 ENDDO
4134 ENDDO
4135 evar(i) = max(evar(i),vg(1),vg(2),vg(3),vg(4),vg(5))
4136 ENDDO ! I=LFT,JLT
4137 ENDIF ! IF (IL <= NLAY .AND. IT <= NPTT)
4138 ENDIF ! IF (MLW == 25 .AND. (IGTYP == 51 .OR. IGTYP == 52 ))
4139C-------------------
4140 ELSEIF (ifunc > 11925+4*mx_ply_anim+306 .AND. ifunc < 11925+4*mx_ply_anim+1317
4141 . .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
4142C-------------------
4143 ! DAMA/ILAY/IPT
4144C-------------------
4145 idx = 11925+4*mx_ply_anim+306
4146 ius = ifunc - idx
4147 il = int((ius - 1)/10)
4148 it = ius - 10*il
4149C
4150 DO i=lft,llt
4151 evar(i) = zero
4152 ENDDO
4153C
4154 IF (ifailure > 0) THEN
4155 IF (il <= nlay) THEN
4156 bufly => elbuf_tab(ng)%BUFLY(il)
4157 nptt = bufly%NPTT
4158 IF (it <= nptt) THEN
4159 DO i = lft,llt
4160 DO ir = 1, nptr
4161 DO is = 1, npts
4162 fbuf => elbuf_tab(ng)%BUFLY(il)%FAIL(ir,is,it)
4163 DO ifail = 1, elbuf_tab(ng)%BUFLY(il)%NFAIL
4164 evar(i) = max(evar(i), fbuf%FLOC(ifail)%DAMMX(i))
4165 ENDDO
4166 ENDDO
4167 ENDDO
4168 ENDDO ! I = LFT,LLT
4169 ENDIF ! IF IT <= NPTT)
4170 ENDIF ! IF (IL <= NLAY )
4171 ENDIF ! IF(IFAILURE > 0) THEN
4172C---
4173C for outp of dam inside law25
4174C---
4175 IF (mlw == 25 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
4176 IF (ity == 3) THEN
4177 DO i=lft,llt
4178 mat(i)=ixc(1,nft+i)
4179 pid(i)=ixc(6,nft+i)
4180 ENDDO
4181 ELSE
4182 DO i=lft,llt
4183 mat(i)=ixtg(1,nft+i)
4184 pid(i)=ixtg(5,nft+i)
4185 ENDDO
4186 ENDIF
4187C
4188 ipmat = 2 + nlay
4189 DO n=1,nlay
4190 iadr = (n-1)*nel
4191 DO i=lft,llt
4192 j = iadr+i
4193 matly(j) = stack%IGEO(ipmat+n,isubstack)
4194 ENDDO
4195 ENDDO
4196C
4197 IF (il <= nlay) THEN
4198 bufly => elbuf_tab(ng)%BUFLY(il)
4199 nptt = bufly%NPTT
4200 IF (it <= nptt) THEN
4201 DO i=lft,llt
4202 iadr = (il - 1)*nel
4203 j = iadr + i
4204 vg(1:5)= zero
4205 DO ir=1,nptr
4206 DO is=1,npts
4207 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
4208 dmax(i) = one/pm(64,matly(j))
4209 wpmax(i)= one/pm(41,matly(j))
4210 epst1(i)= pm(60,matly(j))
4211 epst2(i)= pm(61,matly(j))
4212 epsf1(i)= one/pm(98,matly(j))
4213 epsf2(i)= one/pm(99,matly(j))
4214C
4215 vg(1) = max(vg(1),lbuf%DMG(jj(2)+i)*dmax(i))
4216 vg(2) = max(vg(2),lbuf%DMG(jj(3)+i)*dmax(i))
4217 vg(3) = max(vg(3),abs(lbuf%PLA(i))*wpmax(i))
4218 IF (lbuf%CRAK(jj(1)+i) > zero) vg(4) = max(vg(4),
4219 . (lbuf%CRAK(jj(1)+i)+epst1(i))*epsf1(i))
4220 IF (lbuf%CRAK(jj(2)+i) > zero) vg(5) = max(vg(5),
4221 . (lbuf%CRAK(jj(2)+i)+epst2(i))*epsf2(i))
4222 ENDDO
4223 ENDDO
4224 evar(i) = max(evar(i),vg(1),vg(2),vg(3),vg(4),vg(5))
4225 ENDDO ! I=LFT,JLT
4226 ENDIF ! IF (IT <= NPTT)
4227 ENDIF ! IF (IL <= NLAY)
4228 ENDIF ! IF (MLW == 25 .AND. (IGTYP == 51 .OR. IGTYP == 52 ))
4229
4230 ELSEIF(ifunc == 13242 + 4*mx_ply_anim )THEN
4231 IF(gbuf%G_DT>0)THEN
4232 DO i=lft,llt
4233 evar(i) = gbuf%DT(i)
4234 ENDDO
4235 ENDIF
4236C
4237 ELSEIF(ifunc == 13243 + 4*mx_ply_anim )THEN
4238 IF(gbuf%G_ISMS>0)THEN
4239 DO i=lft,llt
4240 evar(i) = gbuf%ISMS(i)
4241 ENDDO
4242 ENDIF
4243
4244 ELSEIF(ifunc == 13245 + 4*mx_ply_anim .AND. (mlw == 15 .OR. mlw == 25))THEN
4245! it's the plastic work /ANIM/ELEM/WPLA ( ifunc=4*MX_PLY_ANIM +13245)
4246 IF (gbuf%G_PLA > 0) THEN
4247 ! for law25, plastic work < 0 if the layer has reached failure-p !
4248 ilay = 1
4249 IF (nlay > 1) ilay = iabs(nlay)/2 + 1
4250 bufly => elbuf_tab(ng)%BUFLY(ilay)
4251 IF (bufly%L_PLA > 0) THEN
4252 IF (npg > 1) THEN
4253 IF(ity == 3) THEN
4254 IF(igtyp == 51 .OR. igtyp == 52) THEN
4255 nptt = bufly%NPTT
4256 DO is = 1,npts
4257 DO ir = 1,nptr
4258 DO it = 1, nptt
4259 DO i=1,nel
4260 evar(i) = evar(i) + fourth*bufly%LBUF(ir,is,it)%PLA(i)/nptt
4261 ENDDO
4262 ENDDO
4263 ENDDO
4264 ENDDO
4265 ELSE
4266 DO i=1,nel
4267 evar(i) = fourth*(bufly%LBUF(1,1,1)%PLA(i) + bufly%LBUF(2,1,1)%PLA(i) +
4268 . bufly%LBUF(1,2,1)%PLA(i) + bufly%LBUF(2,2,1)%PLA(i))
4269 ENDDO
4270 ENDIF ! igtyp
4271 ELSE ! ITY == 7
4272 IF(igtyp == 51 .OR. igtyp == 52) THEN
4273 nptt = bufly%NPTT
4274 DO it = 1,nptt
4275 DO ir =1,npg
4276 DO i=1,nel
4277 evar(i) = evar(i) + third*bufly%LBUF(ir,1,it)%PLA(i)/nptt
4278 ENDDO
4279 ENDDO
4280 ENDDO
4281 ELSE
4282 DO i=1,nel
4283 evar(i) = third*(bufly%LBUF(1,1,1)%PLA(i) + bufly%LBUF(1,1,1)%PLA(i) +
4284 . bufly%LBUF(1,1,1)%PLA(i))
4285 ENDDO
4286 ENDIF ! igtyp
4287 ENDIF ! ity
4288 ELSE ! NPG == 1
4289 IF(igtyp == 51 .OR. igtyp == 52) THEN
4290 nptt = bufly%NPTT ! needed for PID51 (new shell prop)
4291 DO it=1,nptt
4292 DO i=1,nel
4293 evar(i) = evar(i) + abs(bufly%LBUF(1,1,it)%PLA(i))/nptt
4294 ENDDO
4295 ENDDO
4296 ELSE
4297 nptt = bufly%NPTT !
4298 ipt = iabs(nptt)/2 + 1
4299 DO i=1,nel
4300 evar(i) = abs(bufly%LBUF(1,1,ipt)%PLA(i))/nptt
4301 ENDDO
4302 ENDIF
4303 ENDIF ! npg
4304 ENDIF ! BUFLY%L_PLA
4305 ENDIF ! gbuf
4306c------------------------------------
4307 ELSEIF (ifunc == 13246 + 4*mx_ply_anim .AND. (mlw == 15 .OR. mlw == 25)) THEN ! WPLA/UPPER
4308c------------------------------------
4309 IF (nlay > 1) THEN
4310 il = max(1,npt)
4311 ipt = 1
4312 ELSE
4313 il = 1
4314 ipt = max(1,npt)
4315 ENDIF
4316 bufly => elbuf_tab(ng)%BUFLY(il)
4317 IF (bufly%L_PLA > 0) THEN
4318 IF (npg > 1) THEN
4319 IF (igtyp == 51 .OR. igtyp == 52) ipt = bufly%NPTT
4320 DO i=lft,llt
4321 DO ir=1,nptr
4322 DO is=1,npts
4323 lbuf => bufly%LBUF(ir,is,ipt)
4324 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
4325 ENDDO
4326 ENDDO
4327 ENDDO
4328 ELSE
4329 IF (igtyp == 51 .OR. igtyp == 52) ipt = bufly%NPTT
4330 DO i=lft,llt
4331 evar(i) = abs(bufly%LBUF(1,1,ipt)%PLA(i))
4332 ENDDO
4333 ENDIF
4334 ENDIF
4335c------------------------------------
4336 ELSEIF (ifunc == 13247 + 4*mx_ply_anim .AND. (mlw == 15 .OR. mlw == 25 )) THEN ! WPLA/LOWER
4337c------------------------------------
4338 bufly => elbuf_tab(ng)%BUFLY(1)
4339 IF (bufly%L_PLA > 0) THEN
4340 IF (npg > 1) THEN
4341 DO i=lft,llt
4342 DO ir=1,nptr
4343 DO is=1,npts
4344 lbuf => bufly%LBUF(ir,is,1)
4345 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
4346 ENDDO
4347 ENDDO
4348 ENDDO
4349 ELSE
4350 DO i=lft,llt
4351 evar(i) = abs(bufly%LBUF(1,1,1)%PLA(i))
4352 ENDDO
4353 ENDIF
4354 ENDIF
4355c------------------------------------
4356 ELSEIF (ifunc > 13247 + 4*mx_ply_anim .AND. ifunc <= 13347 + 4*mx_ply_anim .AND.
4357 . (mlw == 15 .OR. mlw == 25)) THEN ! anim/shell/wpla/layer
4358c------------------------------------
4359 ilay = mod((ifunc - 13247 - 4*mx_ply_anim), 100)
4360 IF (ilay == 0) ilay = 100
4361 IF ((ilay <= nlay .or. ilay <= mpt) .and. gbuf%G_PLA > 0) THEN
4362 IF (npt == 0) THEN
4363 il = 1
4364 ipt = 1
4365 ELSEIF (nlay > 1) THEN
4366 il = ilay
4367 ipt = 1
4368 ELSE
4369 il = 1
4370 ipt = ilay
4371 ENDIF
4372 bufly => elbuf_tab(ng)%BUFLY(il)
4373 IF (bufly%L_PLA > 0) THEN
4374 IF (npg > 1) THEN
4375 IF (igtyp == 51 .OR. igtyp == 52) THEN
4376C PROP/51 : more than 1 integration point by layer: average value per layer
4377 nptt = bufly%NPTT
4378 npgt = npg*nptt
4379 DO i=lft,llt
4380 DO it=1,nptt
4381 DO ir=1,nptr
4382 DO is=1,npts
4383 lbuf => bufly%LBUF(ir,is,it)
4384 evar(i) = evar(i) + abs(lbuf%PLA(i))/npgt
4385 ENDDO
4386 ENDDO
4387 ENDDO
4388 ENDDO
4389 ELSE
4390 DO i=lft,llt
4391 DO ir=1,nptr
4392 DO is=1,npts
4393 lbuf => bufly%LBUF(ir,is,ipt)
4394 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
4395 ENDDO
4396 ENDDO
4397 ENDDO
4398 ENDIF
4399 ELSE
4400 IF (igtyp == 51 .OR. igtyp == 52) THEN
4401 nptt = bufly%NPTT
4402 DO i=lft,llt
4403 DO it=1,nptt
4404 evar(i) = evar(i) + abs(bufly%LBUF(1,1,it)%PLA(i))/nptt
4405 ENDDO
4406 ENDDO
4407 ELSE
4408 DO i=lft,llt
4409 evar(i) = abs(bufly%LBUF(1,1,ipt)%PLA(i))
4410 ENDDO
4411 ENDIF
4412 ENDIF
4413 ENDIF
4414 ENDIF
4415c------------------------------------
4416 ELSEIF (ifunc > 13347 + 4*mx_ply_anim .AND. ifunc <= 13447 + 4*mx_ply_anim .AND.
4417 . (igtyp == 51 .OR. igtyp == 52) .AND. (mlw == 15 .OR. mlw == 25) ) THEN
4418 ! WPLA/ILAY/UPPER
4419c------------------------------------
4420C available for IGTYP = 51 and 52 only (to be generalized)
4421 ilay = mod((ifunc - 13347 - 4*mx_ply_anim), 100)
4422 IF (ilay == 0) ilay = 100
4423 IF (nlay > 1) THEN
4424 il = max(1,ilay)
4425 ELSE
4426 il = 1
4427 ENDIF
4428 bufly => elbuf_tab(ng)%BUFLY(il)
4429 nptt = bufly%NPTT
4430 ipt = max(1,nptt)
4431 IF (bufly%L_PLA > 0 .AND.
4432 . (il <= nlay .AND. ipt <= nptt))THEN
4433 IF (npg > 1) THEN
4434 DO i=lft,llt
4435 DO ir=1,nptr
4436 DO is=1,npts
4437 lbuf => bufly%LBUF(ir,is,ipt)
4438 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
4439 ENDDO
4440 ENDDO
4441 ENDDO
4442 ELSE ! (NPG == 1)
4443 lbuf => bufly%LBUF(1,1,ipt)
4444 DO i=lft,llt
4445 evar(i) = abs(lbuf%PLA(i))
4446 ENDDO
4447 ENDIF ! IF (NPG > 1) THEN
4448 ENDIF ! IF (BUFLY%L_PLA > 0) THEN
4449c------------------------------------
4450 ELSEIF (ifunc > 13447 + 4*mx_ply_anim .AND. ifunc <= 13547 + 4*mx_ply_anim .AND.
4451 . (igtyp == 51 .OR. igtyp == 52) .AND. (mlw == 15 .OR. mlw == 25) ) THEN
4452 ! WPLA/ILAY/LOWER
4453c------------------------------------
4454C available for IGTYP = 51 and 52 only (to be generalized)
4455 ilay = mod((ifunc - 13447 - 4*mx_ply_anim), 100)
4456 IF (ilay == 0) ilay = 100
4457 IF (nlay > 1) THEN
4458 il = max(1,ilay)
4459 ELSE
4460 il = 1
4461 ENDIF
4462 ipt = 1
4463 bufly => elbuf_tab(ng)%BUFLY(il)
4464 nptt = bufly%NPTT
4465 IF (bufly%L_PLA > 0 .AND.
4466 . (il <= nlay .AND. ipt <= nptt))THEN
4467 IF (npg > 1) THEN
4468 DO i=lft,llt
4469 DO ir=1,nptr
4470 DO is=1,npts
4471 lbuf => bufly%LBUF(ir,is,ipt)
4472 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
4473 ENDDO
4474 ENDDO
4475 ENDDO
4476 ELSE ! (NPG == 1)
4477 lbuf => bufly%LBUF(1,1,ipt)
4478 DO i=lft,llt
4479 evar(i) = abs(lbuf%PLA(i))
4480 ENDDO
4481 ENDIF ! IF (NPG > 1) THEN
4482 ENDIF ! IF (BUFLY%L_PLA > 0) THEN
4483c------------------------------------
4484 ELSEIF (ifunc > 13547 + 4*mx_ply_anim .AND. ifunc <= 14547 + 4*mx_ply_anim .AND.
4485 . (igtyp == 51 .OR. igtyp == 52) .AND. (mlw == 15 .OR. mlw == 25) ) THEN
4486 ! WPLA/ILAY/NIP
4487c------------------------------------
4488C
4489C available for IGTYP = 51 and 52 only (to be generalized)
4490C
4491 ius = ifunc - 13547 - 4*mx_ply_anim
4492 il = int((ius - 1)/10)
4493 ipt = ius - 10*il
4494 il = il + 1
4495 IF (il <= nlay ) THEN
4496 bufly => elbuf_tab(ng)%BUFLY(il)
4497 nptt = bufly%NPTT
4498 IF (bufly%L_PLA > 0 .AND. ipt <= nptt) THEN
4499 IF (npg > 1) THEN
4500 DO i=lft,llt
4501 DO ir=1,nptr
4502 DO is=1,npts
4503 lbuf => bufly%LBUF(ir,is,ipt)
4504 evar(i) = evar(i) + abs(lbuf%PLA(i))/npg
4505 ENDDO
4506 ENDDO
4507 ENDDO
4508 ELSE ! (NPG == 1)
4509 lbuf => bufly%LBUF(1,1,ipt)
4510 DO i=lft,llt
4511 evar(i) = abs(lbuf%PLA(i))
4512 ENDDO
4513 ENDIF ! IF (NPG > 1) THEN
4514 ENDIF ! IF (BUFLY%L_PLA > 0) THEN
4515 ENDIF ! IF (IL <= NLAY )
4516 ! next IFUNC = 14548 + 4*MX_PLY_ANIM
4517c------------------------------------ OFF
4518 ELSEIF (ifunc == 13547 + 4*mx_ply_anim + 1000 + 1) THEN
4519 DO i=lft,llt
4520 IF (gbuf%G_OFF > 0) THEN
4521 IF(gbuf%OFF(i) > one) THEN
4522 evar(i) = gbuf%OFF(i) - one
4523 ELSEIF((gbuf%OFF(i) >= zero .AND. gbuf%OFF(i) <= one)) THEN
4524 evar(i) = gbuf%OFF(i)
4525 ELSE
4526 evar(i) = -one
4527 ENDIF
4528 ENDIF
4529 ENDDO
4530 ! next IFUNC = 13547 + 4*MX_PLY_ANIM + 1000 + 2
4531c------------------------------------ Mach Number
4532 ELSEIF(ifunc == 13547 + 4*mx_ply_anim + 1000 + 2)THEN
4533 IF (mlw == 151) THEN
4534 DO i = 1, nel
4535 vel(1) = multi_fvm%VEL(1, i + nft)
4536 vel(2) = multi_fvm%VEL(2, i + nft)
4537 vel(3) = multi_fvm%VEL(3, i + nft)
4538 vel(0) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))
4539 evar(i) = vel(0)/multi_fvm%SOUND_SPEED(i + nft)
4540 ENDDO
4541 ELSEIF(alefvm_param%ISOLVER>1)THEN
4542 l = elbuf_tab(ng)%BUFLY(1)%L_SSP
4543 IF(elbuf_tab(ng)%BUFLY(1)%L_SSP /= 0)THEN
4544 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
4545 DO i=1,nel
4546 vel(1) = gbuf%MOM(jj(1) + i) / gbuf%RHO(i)
4547 vel(2) = gbuf%MOM(jj(2) + i) / gbuf%RHO(i)
4548 vel(3) = gbuf%MOM(jj(3) + i) / gbuf%RHO(i)
4549 vel(0) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))
4550 evar(i) = vel(0)/lbuf%SSP(i)
4551 ENDDO
4552 ENDIF
4553 ELSE
4554 l = elbuf_tab(ng)%BUFLY(1)%L_SSP
4555 IF(n2d/=0.AND.elbuf_tab(ng)%BUFLY(1)%L_SSP /= 0)THEN
4556 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
4557 IF(jale/=0)THEN
4558 IF(ity==7) THEN
4559 DO i=1,nel
4560 tmp(1,1:3)=v(1,ixtg(2:4,i+nft))-w(1,ixtg(2:4,i+nft))
4561 tmp(2,1:3)=v(2,ixtg(2:4,i+nft))-w(2,ixtg(2:4,i+nft))
4562 tmp(3,1:3)=v(3,ixtg(2:4,i+nft))-w(3,ixtg(2:4,i+nft))
4563 vel(1) = sum(tmp(1,1:3))*third
4564 vel(2) = sum(tmp(2,1:3))*third
4565 vel(3) = sum(tmp(3,1:3))*third
4566 evar(i) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))/lbuf%SSP(i)
4567 ENDDO
4568 ELSEIF(ity==2) THEN
4569 DO i=1,nel
4570 tmp(1,1:4)=v(1,ixq(2:5,i+nft))-w(1,ixq(2:5,i+nft))
4571 tmp(2,1:4)=v(2,ixq(2:5,i+nft))-w(2,ixq(2:5,i+nft))
4572 tmp(3,1:4)=v(3,ixq(2:5,i+nft))-w(3,ixq(2:5,i+nft))
4573 vel(1) = sum(tmp(1,1:4))*fourth
4574 vel(2) = sum(tmp(2,1:4))*fourth
4575 vel(3) = sum(tmp(3,1:4))*fourth
4576 evar(i) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))/lbuf%SSP(i)
4577 ENDDO
4578 ENDIF
4579 ELSE
4580 IF(ity==7) THEN
4581 DO i=1,nel
4582 tmp(1,1:3)=v(1,ixtg(2:4,i+nft))
4583 tmp(2,1:3)=v(2,ixtg(2:4,i+nft))
4584 tmp(3,1:3)=v(3,ixtg(2:4,i+nft))
4585 vel(1) = sum(tmp(1,1:3))*third
4586 vel(2) = sum(tmp(2,1:3))*third
4587 vel(3) = sum(tmp(3,1:3))*third
4588 evar(i) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))/lbuf%SSP(i)
4589 ENDDO
4590 ELSEIF(ity==2) THEN
4591 DO i=1,nel
4592 tmp(1,1:4)=v(1,ixq(2:5,i+nft))
4593 tmp(2,1:4)=v(2,ixq(2:5,i+nft))
4594 tmp(3,1:4)=v(3,ixq(2:5,i+nft))
4595 vel(1) = sum(tmp(1,1:4))*fourth
4596 vel(2) = sum(tmp(2,1:4))*fourth
4597 vel(3) = sum(tmp(3,1:4))*fourth
4598 evar(i) = sqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3))/lbuf%SSP(i)
4599 ENDDO
4600 ENDIF
4601 ENDIF
4602 ENDIF
4603 ENDIF
4604c------------------------------------ DAMG
4605 ELSEIF((ifunc >= 13547 + 4*mx_ply_anim + 1000 + 4).AND.
4606 . (ifunc <= 13547 + 4*mx_ply_anim + 1000 + 18).AND.gbuf%G_DMG > 0)THEN
4607 idx = 13547 + 4*mx_ply_anim + 1000 + 4
4608c------------Mean
4609 IF (ifunc == idx) THEN
4610 DO i = lft, llt
4611 evar(i) = zero
4612 ENDDO
4613 npgt = npg*nptt
4614 DO il=1,nlay
4615 DO is=1,npts
4616 DO it=1,nptt
4617 DO ir=1,nptr
4618 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
4619 DO i=lft,llt
4620 evar(i) = evar(i) + lbuf%DMG(i)/npgt
4621 ENDDO
4622 ENDDO
4623 ENDDO
4624 ENDDO
4625 ENDDO
4626c------------Upper
4627 ELSEIF (ifunc == idx + 1) THEN
4628 DO i = lft, llt
4629 evar(i) = zero
4630 ENDDO
4631 DO il=1,nlay
4632 DO is=1,npts
4633 DO ir=1,nptr
4634 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,nptt)
4635 DO i=lft,llt
4636 evar(i) = evar(i) + lbuf%DMG(i)/(npg*nlay)
4637 ENDDO
4638 ENDDO
4639 ENDDO
4640 ENDDO
4641c------------Lower
4642 ELSEIF (ifunc == idx + 2) THEN
4643 DO i = lft, llt
4644 evar(i) = zero
4645 ENDDO
4646 DO il=1,nlay
4647 DO is=1,npts
4648 DO ir=1,nptr
4649 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,1)
4650 DO i=lft,llt
4651 evar(i) = evar(i) + lbuf%DMG(i)/(npg*nlay)
4652 ENDDO
4653 ENDDO
4654 ENDDO
4655 ENDDO
4656c----------Membrane
4657 ELSEIF (ifunc == idx + 3) THEN
4658 DO i = lft, llt
4659 evar(i) = zero
4660 ENDDO
4661 ! Odd number of thickness integration points
4662 IF ((mod(nptt,2)/=0).AND.(nptt>1)) THEN
4663 DO il=1,nlay
4664 DO is=1,npts
4665 DO ir=1,nptr
4666 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,ceiling(nptt/two))
4667 DO i=lft,llt
4668 evar(i) = evar(i) + lbuf%DMG(i)/(npg*nlay)
4669 ENDDO
4670 ENDDO
4671 ENDDO
4672 ENDDO
4673 ! Even number of thickness integration points
4674 ELSEIF ((mod(nptt,2)==0).AND.(nptt>1)) THEN
4675 DO il=1,nlay
4676 DO is=1,npts
4677 DO ir=1,nptr
4678 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,nptt/2)
4679 DO i=lft,llt
4680 evar(i) = evar(i) + lbuf%DMG(i)/(two*npg*nlay)
4681 ENDDO
4682 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,nptt/2+1)
4683 DO i=lft,llt
4684 evar(i) = evar(i) + lbuf%DMG(i)/(two*npg*nlay)
4685 ENDDO
4686 ENDDO
4687 ENDDO
4688 ENDDO
4689 ! NPTT = 1
4690 ELSE
4691 DO il=1,nlay
4692 DO is=1,npts
4693 DO ir=1,nptr
4694 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,1)
4695 DO i=lft,llt
4696 evar(i) = evar(i) + lbuf%DMG(i)/(npg*nlay)
4697 ENDDO
4698 ENDDO
4699 ENDDO
4700 ENDDO
4701 ENDIF
4702c------------At a given integration point in the thickness
4703 ELSEIF((ifunc >= idx + 3 + 1).AND.(ifunc <= idx + 3 + 11))THEN
4704 DO i = lft, llt
4705 evar(i) = zero
4706 ENDDO
4707 it = ifunc - (idx+3)
4708 IF (it<=nptt) THEN
4709 DO il=1,nlay
4710 DO is=1,npts
4711 DO ir=1,nptr
4712 lbuf => elbuf_tab(ng)%BUFLY(il)%LBUF(ir,is,it)
4713 DO i=lft,llt
4714 evar(i) = evar(i) + lbuf%DMG(i)/(npg*nlay)
4715 ENDDO
4716 ENDDO
4717 ENDDO
4718 ENDDO
4719 ENDIF
4720 ENDIF
4721c------------------------------------ NL_EPSP
4722 ELSEIF((ifunc >= 14567 + 4*mx_ply_anim).AND.
4723 . (ifunc <= 14580 + 4*mx_ply_anim).AND.
4724 . gbuf%G_PLANL > 0)THEN
4725 idx = 14567 + 4*mx_ply_anim
4726c------------Mean
4727 IF (ifunc == idx) THEN
4728 DO i = lft, llt
4729 evar(i) = zero
4730 ENDDO
4731 npgt = npg*nptt
4732 ! (Only 1 layer is supported with non-local)
4733 DO is=1,npts
4734 DO it=1,nptt
4735 DO ir=1,nptr
4736 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,it)
4737 DO i=lft,llt
4738 evar(i) = evar(i) + lbuf%PLANL(i)/npgt
4739 ENDDO
4740 ENDDO
4741 ENDDO
4742 ENDDO
4743c------------Upper
4744 ELSEIF (ifunc == idx + 1) THEN
4745 DO i = lft, llt
4746 evar(i) = zero
4747 ENDDO
4748 ! (Only 1 layer is supported with non-local)
4749 DO is=1,npts
4750 DO ir=1,nptr
4751 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,nptt)
4752 DO i=lft,llt
4753 evar(i) = evar(i) + lbuf%PLANL(i)/npg
4754 ENDDO
4755 ENDDO
4756 ENDDO
4757c------------Lower
4758 ELSEIF (ifunc == idx + 2) THEN
4759 DO i = lft, llt
4760 evar(i) = zero
4761 ENDDO
4762 ! (Only 1 layer is supported with non-local)
4763 DO is=1,npts
4764 DO ir=1,nptr
4765 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
4766 DO i=lft,llt
4767 evar(i) = evar(i) + lbuf%PLANL(i)/npg
4768 ENDDO
4769 ENDDO
4770 ENDDO
4771c------------At a given integration point in the thickness
4772 ELSEIF((ifunc >= idx + 2 + 1).AND.(ifunc <= idx + 2 + 11))THEN
4773 DO i = lft, llt
4774 evar(i) = zero
4775 ENDDO
4776 it = ifunc - (idx+2)
4777 ! (Only 1 layer is supported with non-local)
4778 IF (it<=nptt) THEN
4779 DO is=1,npts
4780 DO ir=1,nptr
4781 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,it)
4782 DO i=lft,llt
4783 evar(i) = evar(i) + lbuf%PLANL(i)/npg
4784 ENDDO
4785 ENDDO
4786 ENDDO
4787 ENDIF
4788 ENDIF
4789c------------------------------------ NL_EPSD
4790 ELSEIF((ifunc >= 14581 + 4*mx_ply_anim).AND.
4791 . (ifunc <= 14594 + 4*mx_ply_anim).AND.
4792 . gbuf%G_EPSDNL > 0)THEN
4793 idx = 14581 + 4*mx_ply_anim
4794c------------Mean
4795 IF (ifunc == idx) THEN
4796 DO i = lft, llt
4797 evar(i) = zero
4798 ENDDO
4799 npgt = npg*nptt
4800 ! (Only 1 layer is supported with non-local)
4801 DO is=1,npts
4802 DO it=1,nptt
4803 DO ir=1,nptr
4804 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,it)
4805 DO i=lft,llt
4806 evar(i) = evar(i) + lbuf%EPSDNL(i)/npgt
4807 ENDDO
4808 ENDDO
4809 ENDDO
4810 ENDDO
4811c------------Upper
4812 ELSEIF (ifunc == idx + 1) THEN
4813 DO i = lft, llt
4814 evar(i) = zero
4815 ENDDO
4816 ! (Only 1 layer is supported with non-local)
4817 DO is=1,npts
4818 DO ir=1,nptr
4819 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,nptt)
4820 DO i=lft,llt
4821 evar(i) = evar(i) + lbuf%EPSDNL(i)/npg
4822 ENDDO
4823 ENDDO
4824 ENDDO
4825c------------Lower
4826 ELSEIF (ifunc == idx + 2) THEN
4827 DO i = lft, llt
4828 evar(i) = zero
4829 ENDDO
4830 ! (Only 1 layer is supported with non-local)
4831 DO is=1,npts
4832 DO ir=1,nptr
4833 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
4834 DO i=lft,llt
4835 evar(i) = evar(i) + lbuf%EPSDNL(i)/npg
4836 ENDDO
4837 ENDDO
4838 ENDDO
4839c------------At a given integration point in the thickness
4840 ELSEIF((ifunc >= idx + 2 + 1).AND.(ifunc <= idx + 2 + 11))THEN
4841 DO i = lft, llt
4842 evar(i) = zero
4843 ENDDO
4844 it = ifunc - (idx+2)
4845 ! (Only 1 layer is supported with non-local)
4846 IF (it<=nptt) THEN
4847 DO is=1,npts
4848 DO ir=1,nptr
4849 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,it)
4850 DO i=lft,llt
4851 evar(i) = evar(i) + lbuf%EPSDNL(i)/npg
4852 ENDDO
4853 ENDDO
4854 ENDDO
4855 ENDIF
4856 ENDIF
4857c------------------------------------ TSAIWU
4858 ELSEIF (ifunc == 14595 + 4*mx_ply_anim .AND. (gbuf%G_TSAIWU > 0)) THEN
4859 ! /ANIM/ELEM/TSAIWU ( IFUNC = 4*MX_PLY_ANIM + 14595)
4860 IF (nlay > 1) THEN
4861 ipt = iabs(nlay)/2 + 1
4862 bufly => elbuf_tab(ng)%BUFLY(ipt)
4863 nptt = bufly%NPTT ! needed for PID51 (new shell prop)
4864 DO i=lft,llt
4865 DO ir=1,nptr
4866 DO is=1,npts
4867 DO it=1,nptt
4868 evar(i) = evar(i) + bufly%LBUF(ir,is,it)%TSAIWU(i)/(nptt*nptr*npts)
4869 ENDDO
4870 ENDDO
4871 ENDDO
4872 ENDDO
4873 ELSE
4874 bufly => elbuf_tab(ng)%BUFLY(1)
4875 IF (bufly%L_TSAIWU > 0) THEN
4876 nptt = bufly%NPTT ! needed for PID51 (new shell prop)
4877 ipt = iabs(nptt)/2 + 1
4878 DO i=lft,llt
4879 DO ir=1,nptr
4880 DO is=1,npts
4881 evar(i) = evar(i) + bufly%LBUF(ir,is,ipt)%TSAIWU(i)/(nptr*npts)
4882 ENDDO
4883 ENDDO
4884 ENDDO
4885 ENDIF
4886 ENDIF
4887c------------------------------------
4888 ELSEIF (ifunc == 14596 + 4*mx_ply_anim .AND. (gbuf%G_TSAIWU > 0)) THEN ! TSAIWU/UPPER
4889c------------------------------------
4890 IF (nlay > 1) THEN
4891 il = max(1,npt)
4892 ipt = 1
4893 ELSE
4894 il = 1
4895 ipt = max(1,npt)
4896 ENDIF
4897 bufly => elbuf_tab(ng)%BUFLY(il)
4898 IF (bufly%L_TSAIWU > 0) THEN
4899 IF (npg > 1) THEN
4900 IF (igtyp == 51 .OR. igtyp == 52) ipt = bufly%NPTT
4901 DO i=lft,llt
4902 DO ir=1,nptr
4903 DO is=1,npts
4904 lbuf => bufly%LBUF(ir,is,ipt)
4905 evar(i) = evar(i) + lbuf%TSAIWU(i)/npg
4906 ENDDO
4907 ENDDO
4908 ENDDO
4909 ELSE
4910 IF (igtyp == 51 .OR. igtyp == 52) ipt = bufly%NPTT
4911 DO i=lft,llt
4912 evar(i) = bufly%LBUF(1,1,ipt)%TSAIWU(i)
4913 ENDDO
4914 ENDIF
4915 ENDIF
4916c------------------------------------
4917 ELSEIF (ifunc == 14597 + 4*mx_ply_anim .AND. (gbuf%G_TSAIWU > 0)) THEN ! TSAIWU/LOWER
4918c------------------------------------
4919 bufly => elbuf_tab(ng)%BUFLY(1)
4920 IF (bufly%L_TSAIWU > 0) THEN
4921 IF (npg > 1) THEN
4922 DO i=lft,llt
4923 DO ir=1,nptr
4924 DO is=1,npts
4925 lbuf => bufly%LBUF(ir,is,1)
4926 evar(i) = evar(i) + lbuf%TSAIWU(i)/npg
4927 ENDDO
4928 ENDDO
4929 ENDDO
4930 ELSE
4931 DO i=lft,llt
4932 evar(i) = bufly%LBUF(1,1,1)%TSAIWU(i)
4933 ENDDO
4934 ENDIF
4935 ENDIF
4936c------------------------------------
4937 ELSEIF (ifunc > 14597 + 4*mx_ply_anim .AND. ifunc <= 14697 + 4*mx_ply_anim .AND.
4938 . (gbuf%G_TSAIWU > 0)) THEN ! anim/shell/tsaiwu/layer
4939c------------------------------------
4940 ilay = mod((ifunc - 14597 - 4*mx_ply_anim), 100)
4941 IF (ilay == 0) ilay = 100
4942 IF ((ilay <= nlay .OR. ilay <= mpt) .AND. gbuf%G_TSAIWU > 0) THEN
4943 IF (npt == 0) THEN
4944 il = 1
4945 ipt = 1
4946 ELSEIF (nlay > 1) THEN
4947 il = ilay
4948 ipt = 1
4949 ELSE
4950 il = 1
4951 ipt = ilay
4952 ENDIF
4953 bufly => elbuf_tab(ng)%BUFLY(il)
4954 IF (bufly%L_TSAIWU > 0) THEN
4955 IF (npg > 1) THEN
4956 IF (igtyp == 51 .OR. igtyp == 52) THEN
4957C PROP/51 : more than 1 integration point by layer: average value per layer
4958 nptt = bufly%NPTT
4959 npgt = npg*nptt
4960 DO i=lft,llt
4961 DO it=1,nptt
4962 DO ir=1,nptr
4963 DO is=1,npts
4964 lbuf => bufly%LBUF(ir,is,it)
4965 evar(i) = evar(i) + lbuf%TSAIWU(i)/npgt
4966 ENDDO
4967 ENDDO
4968 ENDDO
4969 ENDDO
4970 ELSE
4971 DO i=lft,llt
4972 DO ir=1,nptr
4973 DO is=1,npts
4974 lbuf => bufly%LBUF(ir,is,ipt)
4975 evar(i) = evar(i) + lbuf%TSAIWU(i)/npg
4976 ENDDO
4977 ENDDO
4978 ENDDO
4979 ENDIF
4980 ELSE
4981 IF (igtyp == 51 .OR. igtyp == 52) THEN
4982 nptt = bufly%NPTT
4983 DO i = lft,llt
4984 DO it = 1,nptt
4985 evar(i) = evar(i) + bufly%LBUF(1,1,it)%TSAIWU(i)/nptt
4986 ENDDO
4987 ENDDO
4988 ELSE
4989 DO i=lft,llt
4990 evar(i) = bufly%LBUF(1,1,ipt)%TSAIWU(i)
4991 ENDDO
4992 ENDIF
4993 ENDIF
4994 ENDIF
4995 ENDIF
4996c------------------------------------
4997 ELSEIF (ifunc > 14697 + 4*mx_ply_anim .AND. ifunc <= 14797 + 4*mx_ply_anim .AND.
4998 . (igtyp == 51 .OR. igtyp == 52) .AND. (gbuf%G_TSAIWU > 0) ) THEN
4999 ! TSAIWU/ILAY/UPPER
5000c------------------------------------
5001C available for IGTYP = 51 and 52 only (to be generalized)
5002 ilay = mod((ifunc - 14697 - 4*mx_ply_anim), 100)
5003 IF (ilay == 0) ilay = 100
5004 IF (nlay > 1) THEN
5005 il = max(1,ilay)
5006 ELSE
5007 il = 1
5008 ENDIF
5009 bufly => elbuf_tab(ng)%BUFLY(il)
5010 nptt = bufly%NPTT
5011 ipt = max(1,nptt)
5012 IF (bufly%L_TSAIWU > 0 .AND.
5013 . (il <= nlay .AND. ipt <= nptt))THEN
5014 IF (npg > 1) THEN
5015 DO i=lft,llt
5016 DO ir=1,nptr
5017 DO is=1,npts
5018 lbuf => bufly%LBUF(ir,is,ipt)
5019 evar(i) = evar(i) + lbuf%TSAIWU(i)/npg
5020 ENDDO
5021 ENDDO
5022 ENDDO
5023 ELSE ! (NPG == 1)
5024 lbuf => bufly%LBUF(1,1,ipt)
5025 DO i=lft,llt
5026 evar(i) = lbuf%TSAIWU(i)
5027 ENDDO
5028 ENDIF ! IF (NPG > 1) THEN
5029 ENDIF ! IF (BUFLY%L_TSAIWU > 0) THEN
5030c------------------------------------
5031 ELSEIF (ifunc > 14797 + 4*mx_ply_anim .AND. ifunc <= 14897 + 4*mx_ply_anim .AND.
5032 . (igtyp == 51 .OR. igtyp == 52) .AND. (gbuf%G_TSAIWU > 0) ) THEN
5033 ! WPLA/ILAY/LOWER
5034c------------------------------------
5035C available for IGTYP = 51 and 52 only (to be generalized)
5036 ilay = mod((ifunc - 14797 - 4*mx_ply_anim), 100)
5037 IF (ilay == 0) ilay = 100
5038 IF (nlay > 1) THEN
5039 il = max(1,ilay)
5040 ELSE
5041 il = 1
5042 ENDIF
5043 ipt = 1
5044 bufly => elbuf_tab(ng)%BUFLY(il)
5045 nptt = bufly%NPTT
5046 IF (bufly%L_TSAIWU > 0 .AND.
5047 . (il <= nlay .AND. ipt <= nptt))THEN
5048 IF (npg > 1) THEN
5049 DO i=lft,llt
5050 DO ir=1,nptr
5051 DO is=1,npts
5052 lbuf => bufly%LBUF(ir,is,ipt)
5053 evar(i) = evar(i) + lbuf%TSAIWU(i)/npg
5054 ENDDO
5055 ENDDO
5056 ENDDO
5057 ELSE ! (NPG == 1)
5058 lbuf => bufly%LBUF(1,1,ipt)
5059 DO i=lft,llt
5060 evar(i) = lbuf%TSAIWU(i)
5061 ENDDO
5062 ENDIF ! IF (NPG > 1) THEN
5063 ENDIF ! IF (BUFLY%L_TSAIWU > 0) THEN
5064c------------------------------------
5065 ELSEIF (ifunc > 14897 + 4*mx_ply_anim .AND. ifunc <= 15897 + 4*mx_ply_anim .AND.
5066 . (igtyp == 51 .OR. igtyp == 52) .AND. (gbuf%G_TSAIWU > 0) ) THEN
5067 ! TSAIWU/ILAY/NIP
5068c------------------------------------
5069C
5070C available for IGTYP = 51 and 52 only (to be generalized)
5071C
5072 ius = ifunc - 14897 - 4*mx_ply_anim
5073 il = int((ius - 1)/10)
5074 ipt = ius - 10*il
5075 il = il + 1
5076 IF (il <= nlay ) THEN
5077 bufly => elbuf_tab(ng)%BUFLY(il)
5078 nptt = bufly%NPTT
5079 IF (bufly%L_TSAIWU > 0 .AND. ipt <= nptt) THEN
5080 IF (npg > 1) THEN
5081 DO i=lft,llt
5082 DO ir=1,nptr
5083 DO is=1,npts
5084 lbuf => bufly%LBUF(ir,is,ipt)
5085 evar(i) = evar(i) + lbuf%TSAIWU(i)/npg
5086 ENDDO
5087 ENDDO
5088 ENDDO
5089 ELSE ! (NPG == 1)
5090 lbuf => bufly%LBUF(1,1,ipt)
5091 DO i=lft,llt
5092 evar(i) = lbuf%TSAIWU(i)
5093 ENDDO
5094 ENDIF ! IF (NPG > 1) THEN
5095 ENDIF ! IF (BUFLY%L_TSAIWU > 0) THEN
5096 ENDIF ! IF (IL <= NLAY )
5097C
5098 ! next IFUNC = 15898 + 4*MX_PLY_ANIM
5099C-------------------------------------
5100 ENDIF ! IFUNC SHELLS
5101C-------------------------------------
5102 IF (mlw == 0 .OR. mlw == 13) THEN
5103 IF (ity == 3) THEN
5104 DO i=lft,llt
5105 n = i + nft
5106 func(el2fa(nn4+n)) = zero
5107 ENDDO
5108 ELSE ! ITY = 7
5109 DO i=lft,llt
5110 n = i + nft
5111 func(el2fa(nn5+n)) = zero
5112 ENDDO
5113 ENDIF
5114C-------------------
5115 ELSEIF (ifunc == 3 .AND. mlw /= 151) THEN
5116C Specific energy
5117C-------------------
5118 IF (ity == 3) THEN
5119 DO i=lft,llt
5120 n = i + nft
5121 func(el2fa(nn4+n)) = evar(i)/
5122 . max(em30,mass(el2fa(nn4+n)))
5123 ENDDO
5124 ELSE ! ITY = 7
5125 DO i=lft,llt
5126 n = i + nft
5127 func(el2fa(nn5+n)) = evar(i)/
5128 . max(em30,mass(el2fa(nn5+n)))
5129 ENDDO
5130 ENDIF
5131C-------------------
5132 ELSEIF (ifunc == 25.AND.ity == 3) THEN
5133C energie hourglass
5134C-------------------
5135 DO i=lft,llt
5136 n = i + nft
5137 func(el2fa(nn4+n)) = ehour(n+numels)/
5138 . max(em30,mass(el2fa(nn4+n)))
5139 ENDDO
5140C-------------------
5141 ELSE ! IFUNC SHELLS
5142C cas general
5143C-------------------
5144 IF(ity == 3)THEN
5145 DO i=lft,llt
5146 n = i + nft
5147 func(el2fa(nn4+n)) = evar(i)
5148 ENDDO
5149 ELSE ! ITY = 7
5150 DO i=lft,llt
5151 n = i + nft
5152 func(el2fa(nn5+n)) = evar(i)
5153 ENDDO
5154 ENDIF
5155 ENDIF ! IFUNC
5156C-----------------------------------------------
5157 ENDIF ! ITY
5158C-----------------------------------------------
5159 END DO ! OFFSET
5160 ENDIF ! MLW /= 13
5161 enddo!next NG
5162C-----------------------------------------------
5163 IF (nspmd == 1) THEN
5164 DO n=1,nbf
5165 r4 = func(n)
5166 CALL write_r_c(r4,1)
5167 ENDDO
5168 ELSE
5169 DO n = 1, nbf_l
5170 wal(n) = func(n)
5171 ENDDO
5172C
5173 IF (ispmd == 0) THEN
5174 buf = (numelqg+numelcg+numeltgg)*4
5175 ELSE
5176 buf=1
5177 ENDIF
5178 CALL spmd_r4get_partn(1,nbf_l,nbpart,iadg,wal,buf)
5179 ENDIF
5180c-----------
5181 IF(ALLOCATED(wa_l))DEALLOCATE(wa_l)
5182 DEALLOCATE(wal)
5183 RETURN
if
if(complex_arithmetic) id
clsconv3
subroutine clsconv3(rx, ry, rz, sx, sy, sz, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)
Definition dfuncc.F:5195
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
ale_connectivity_mod
Definition ale_connectivity_mod.F:225
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_param
type(alefvm_param_), target alefvm_param
Definition alefvm_mod.F:121
initbuf_mod
Definition initbuf.F:160
initbuf_mod::initbuf
subroutine initbuf(iparg, ng, mtn, llt, nft, iad, ity, npt, jale, ismstr, jeul, jtur, jthe, jlag, jmult, jhbe, jivf, mid, jpor, jcvt, jclose, jpla, irep, iint, igtyp, israt, isrot, icsen, isorth, isorthg, ifailure, jsms)
Definition initbuf.F:261
schlieren_mod
Definition schlieren_mod.F:46
stack_mod
Definition stack_mod.F:75
stack_mod::ply_info
integer, dimension(:,:), allocatable ply_info
Definition stack_mod.F:133
output_schlieren
subroutine output_schlieren(evar, ix, x, iparg, wa_l, elbuf_tab, ale_connectivity, vol, ng, nix, ityp)
Definition output_schlieren.F:40
schlieren_buffer_gathering
subroutine schlieren_buffer_gathering(nercvois, nesdvois, lercvois, lesdvois, iparg, elbuf_tab, multi_fvm, itherm)
Definition schlieren_buffer_gathering.F:44
spmd_r4get_partn
subroutine spmd_r4get_partn(size, nbf_l, nbpart, iadg, wal, buf)
Definition spmd_r4get_partn.F:47
write_r_c
void write_r_c(float *w, int *len)
Definition write_routines.c:519