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