OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps33.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!|| sigeps33 ../engine/source/materials/mat/mat033/sigeps33.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
28!||--- calls -----------------------------------------------------
29!|| dreh ../engine/source/materials/mat/mat033/sigeps33.F
30!|| finter ../engine/source/tools/curve/finter.F
31!|| jacobiew ../engine/source/materials/mat/mat033/sigeps33.F
32!|| valpvec_v ../engine/source/materials/mat/mat033/sigeps33.F
33!|| valpvecdp_v ../engine/source/materials/mat/mat033/sigeps33.F
34!||====================================================================
35 SUBROUTINE sigeps33(
36 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC , NPF ,
37 2 TF , TIME , TIMESTEP, UPARAM, RHO0 , RHO ,
38 3 VOLUME , EINT ,
39 4 EPSPXX , EPSPYY , EPSPZZ , EPSPXY, EPSPYZ, EPSPZX,
40 5 DEPSXX , DEPSYY , DEPSZZ , DEPSXY, DEPSYZ, DEPSZX,
41 6 EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
42 7 SIGOXX , SIGOYY , SIGOZZ , SIGOXY, SIGOYZ, SIGOZX,
43 8 SIGNXX , SIGNYY , SIGNZZ , SIGNXY, SIGNYZ, SIGNZX,
44 9 SIGVXX , SIGVYY , SIGVZZ , SIGVXY, SIGVYZ, SIGVZX,
45 A SOUNDSP, VISCMAX, UVAR , OFF )
46
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51#include "scr05_c.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56#include "com01_c.inc"
57C----------------------------------------------------------------
58C I N P U T A R G U M E N T S
59C----------------------------------------------------------------
60 INTEGER NEL, NUPARAM, NUVAR
61 my_real
62 . TIME , TIMESTEP , UPARAM(NUPARAM),
63 . RHO (NEL), RHO0 (NEL), VOLUME(NEL), EINT(NEL),
64 . EPSPXX(NEL), EPSPYY(NEL), EPSPZZ(NEL),
65 . EPSPXY(NEL), EPSPYZ(NEL), EPSPZX(NEL),
66 . DEPSXX(NEL), DEPSYY(NEL), DEPSZZ(NEL),
67 . DEPSXY(NEL), DEPSYZ(NEL), DEPSZX(NEL),
68 . EPSXX (NEL), EPSYY (NEL), EPSZZ (NEL),
69 . EPSXY (NEL), EPSYZ (NEL), EPSZX (NEL),
70 . sigoxx(nel), sigoyy(nel), sigozz(nel),
71 . sigoxy(nel), sigoyz(nel), sigozx(nel)
72C----------------------------------------------------------------
73C O U T P U T A R G U M E N T S
74C----------------------------------------------------------------
75 my_real
76 . signxx(nel), signyy(nel), signzz(nel),
77 . signxy(nel), signyz(nel), signzx(nel),
78 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
79 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
80 . soundsp(nel), viscmax(nel)
81C----------------------------------------------------------------
82C I N P U T O U T P U T A R G U M E N T S
83C----------------------------------------------------------------
84 my_real uvar(nel,nuvar), off(nel)
85C----------------------------------------------------------------
86C VARIABLES FOR FUNCTION INTERPOLATION
87C----------------------------------------------------------------
88 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
89 my_real FINTER,TF(*)
90 EXTERNAL FINTER
91C----------------------------------------------------------------
92C L O C A L V A R I B L E S
93C----------------------------------------------------------------
94 INTEGER I,J,K,L,NROT,KEN
95 INTEGER VPVEC,ICASE
96 my_real
97 . EY,C1,C2,ET,VMU,VMU0
98 . , a,b,c,var
99 .
100 .
101 .
102 .
103 . , p0,phi,gama0,fac,fac1
104 . , s(3,3),sigpr(3),dirpr(3,3)
105 . ,sigt_cutoff,sigc_cutoff
106 my_real
107 . e(mvsiz),edot(mvsiz),dav,e1,e2,e3,e4,e5,e6,epsp
108 . , epet(mvsiz),emet(mvsiz),epets(mvsiz),emets(mvsiz)
109 .
110 .
111 .
112 .
113 .
114 .
115 .
116 .
117 . , sigair(mvsiz),gama(mvsiz),syield(mvsiz)
118 . , sigsxx(mvsiz),sigsyy(mvsiz),sigszz(mvsiz)
119 . , sigsxy(mvsiz),sigsyz(mvsiz),sigszx(mvsiz)
120 . , dexx(mvsiz),deyy(mvsiz),dezz(mvsiz)
121 . , dexy(mvsiz),deyz(mvsiz),dezx(mvsiz)
122 . , dedtxx(mvsiz),dedtyy(mvsiz),dedtzz(mvsiz)
123 . , dedtxy(mvsiz),dedtyz(mvsiz),dedtzx(mvsiz)
124 . , dsdtxx(mvsiz),dsdtyy(mvsiz),dsdtzz(mvsiz)
125 . , dsdtxy(mvsiz),dsdtyz(mvsiz),dsdtzx(mvsiz)
126C
127 my_real
128 . sigv(mvsiz,6), sigprv(mvsiz,3), dirprv(mvsiz,3,3)
129C
130
131C----------------------------------------------------------------
132
133C INITIALIZATION
134 IF(isigi == 0) THEN
135 IF(time==zero) THEN
136 DO 1000 j=1,nuvar
137 DO 1000 i=1,nel
138 uvar(i,j)=zero
139 1000 CONTINUE
140 ENDIF
141 ENDIF
142
143 ken = uparam(1)
144 ey = uparam(2)
145 a = uparam(3)
146 b = uparam(4)
147 c = uparam(5)
148 p0 = uparam(6)
149 phi = uparam(7)
150 gama0 = uparam(8)
151 !!
152 vpvec = 0
153C COMPUTE VOLUMETRIC STRAIN-AIR PRESSURE
154 DO 900 i=1,nel
155 gama(i) = rho0(i)/rho(i)-one+gama0
156 var = -(p0*gama(i))/(one+gama(i)-phi+em15)
157 sigair(i) = max(zero,var)
158 900 CONTINUE
159 icase = abs(ken) + 1
160 SELECT CASE (icase)
161C----------------------------------
162 CASE(1)
163C----------------------------------
164 fac = uparam(9)
165 fac1 = uparam(10)
166
167 DO 100 i=1,nel
168 IF(ifunc(1)/=0) THEN
169 syield(i)=fac*finter(ifunc(1),gama(i),npf,tf,b*c)
170 ELSE
171 syield(i) = abs(a+b*(one+c*gama(i)))
172 ENDIF
173 IF(ifunc(2)/=0)THEN
174 dav = (epspxx(i)+epspyy(i)+epspzz(i))*third
175 e1 = epspxx(i) - dav
176 e2 = epspyy(i) - dav
177 e3 = epspzz(i) - dav
178 e4 = half*epspxy(i)
179 e5 = half*epspyz(i)
180 e6 = half*epspzx(i)
181 epsp =half*(e1**2+e2**2+e3**2) +e4**2+e5**2+e6**2
182 epsp = sqrt(three*epsp)/three_half
183 syield(i) = fac1*syield(i)*(finter(ifunc(2),epsp,npf,tf,b*c))
184 ENDIF
185 100 CONTINUE
186
187 DO 110 i=1,nel
188 sigsxx(i)=sigoxx(i)+sigair(i)
189 sigsyy(i)=sigoyy(i)+sigair(i)
190 sigszz(i)=sigozz(i)+sigair(i)
191 sigsxy(i)=sigoxy(i)
192 sigsyz(i)=sigoyz(i)
193 sigszx(i)=sigozx(i)
194 110 CONTINUE
195
196C COMP UTE TRIAL STRESSES
197 DO 120 i=1,nel
198 sigsxx(i)=sigsxx(i)+ey*epspxx(i)*timestep
199 sigsyy(i)=sigsyy(i)+ey*epspyy(i)*timestep
200 sigszz(i)=sigszz(i)+ey*epspzz(i)*timestep
201 sigsxy(i)=sigsxy(i)+ey*epspxy(i)*timestep*half
202 sigsyz(i)=sigsyz(i)+ey*epspyz(i)*timestep*half
203 sigszx(i)=sigszx(i)+ey*epspzx(i)*timestep*half
204 120 CONTINUE
205C If old method then ....
206 IF (vpvec==1) THEN
207 DO 130 i=1,nel
208 s(1,1)=sigsxx(i)
209 s(2,1)=sigsxy(i)
210 s(2,2)=sigsyy(i)
211 s(3,1)=sigszx(i)
212 s(3,2)=sigsyz(i)
213 s(3,3)=sigszz(i)
214
215C TRANSFORM GLOBAL STRESSES (S_IJ) TO PRINCIPAL STRESSES (SIGPR_K)
216 CALL jacobiew(s,3,sigpr,dirpr,nrot)
217
218C YIELD CRITERIA - SCALING
219 sigpr(1)=min(syield(i),abs(sigpr(1)))*sign(one,sigpr(1))
220 sigpr(2)=min(syield(i),abs(sigpr(2)))*sign(one,sigpr(2))
221 sigpr(3)=min(syield(i),abs(sigpr(3)))*sign(one,sigpr(3))
222
223C TRANSFORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
224 DO 131 k=1,3
225 DO 131 l=1,3
226 s(k,l)=zero
227 IF(k==l) s(k,l)= sigpr(k)
228 131 CONTINUE
229 CALL dreh(s,dirpr,1,1,1)
230
231 sigsxx(i)=s(1,1)
232 sigsxy(i)=s(2,1)
233 sigsyy(i)=s(2,2)
234 sigszx(i)=s(3,1)
235 sigsyz(i)=s(3,2)
236 sigszz(i)=s(3,3)
237 130 CONTINUE
238C fin ancienne methode
239C sino n
240 ELSE
241 DO i = 1, nel
242 sigv(i,1) = sigsxx(i)
243 sigv(i,2) = sigsyy(i)
244 sigv(i,3) = sigszz(i)
245 sigv(i,4) = sigsxy(i)
246 sigv(i,5) = sigsyz(i)
247 sigv(i,6) = sigszx(i)
248 ENDDO
249C for a simple precision executing
250 IF (iresp==1) THEN
251 CALL valpvecdp_v(sigv,sigprv,dirprv,nel)
252 ELSE
253 CALL valpvec_v(sigv,sigprv,dirprv,nel)
254 ENDIF
255 DO i = 1, nel
256C YIEL D CRITERIA - SCALING
257 sigprv(i,1)=sign(min(syield(i),abs(sigprv(i,1))),sigprv(i,1))
258 sigprv(i,2)=sign(min(syield(i),abs(sigprv(i,2))),sigprv(i,2))
259 sigprv(i,3)=sign(min(syield(i),abs(sigprv(i,3))),sigprv(i,3))
260C TRANSFORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
261 sigsxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigprv(i,1)
262 . + dirprv(i,1,2)*dirprv(i,1,2)*sigprv(i,2)
263 . + dirprv(i,1,3)*dirprv(i,1,3)*sigprv(i,3)
264 sigsyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigprv(i,2)
265 . + dirprv(i,2,3)*dirprv(i,2,3)*sigprv(i,3)
266 . + dirprv(i,2,1)*dirprv(i,2,1)*sigprv(i,1)
267 sigszz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigprv(i,3)
268 . + dirprv(i,3,1)*dirprv(i,3,1)*sigprv(i,1)
269 . + dirprv(i,3,2)*dirprv(i,3,2)*sigprv(i,2)
270 sigsxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigprv(i,1)
271 . + dirprv(i,1,2)*dirprv(i,2,2)*sigprv(i,2)
272 . + dirprv(i,1,3)*dirprv(i,2,3)*sigprv(i,3)
273 sigsyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigprv(i,2)
274 . + dirprv(i,2,3)*dirprv(i,3,3)*sigprv(i,3)
275 . + dirprv(i,2,1)*dirprv(i,3,1)*sigprv(i,1)
276 sigszx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigprv(i,3)
277 . + dirprv(i,3,1)*dirprv(i,1,1)*sigprv(i,1)
278 . + dirprv(i,3,2)*dirprv(i,1,2)*sigprv(i,2)
279 ENDDO
280 ENDIF
281C
282C COMPUTE UPDATED STRESSES AND SOUND SPEED
283 DO 140 i=1,nel
284 signxx(i)=off(i)*sigsxx(i)-sigair(i)
285 signyy(i)=off(i)*sigsyy(i)-sigair(i)
286 signzz(i)=off(i)*sigszz(i)-sigair(i)
287 signxy(i)=off(i)*sigsxy(i)
288 signyz(i)=off(i)*sigsyz(i)
289 signzx(i)=off(i)*sigszx(i)
290
291 soundsp(i)=sqrt(ey/rho0(i))
292 140 CONTINUE
293
294C--------------------------
295 CASE(2)
296C--------------------------
297 c1 = uparam(9)
298 c2 = uparam(10)
299 et = uparam(11)
300 vmu = uparam(12)
301 vmu0 = uparam(13)
302 fac = uparam(14)
303 fac1 = uparam(15)
304
305C COMPUTE YOUNG MODULUS, ETC.
306 DO 200 i=1,nel
307 edot(i)=max(
308 & abs(epspxx(i)),abs(epspyy(i)),abs(epspzz(i)),
309 & abs(epspxy(i)),abs(epspyz(i)),abs(epspzx(i)))
310 e(i)=c1*edot(i)+c2
311 e(i)=max(e(i),ey)
312 epet(i)=(e(i)+et)/vmu
313 emet(i)=(e(i)*et)/vmu
314 epets(i)=(e(i)+et)/vmu0
315 emets(i)=two*(e(i)*et)/vmu0
316 200 CONTINUE
317
318C COMPUTE YIELD STRESS
319 DO 210 i=1,nel
320 IF(ifunc(1)/=0) THEN
321 syield(i)=fac*finter(ifunc(1),gama(i),npf,tf,b*c)
322 ELSE
323 syield(i) = abs(a+b*(one+c*gama(i)))
324 ENDIF
325 IF(ifunc(2)/=0)THEN
326 dav = (epspxx(i)+epspyy(i)+epspzz(i))*third
327 e1 = epspxx(i) - dav
328 e2 = epspyy(i) - dav
329 e3 = epspzz(i) - dav
330 e4 = half*epspxy(i)
331 e5 = half*epspyz(i)
332 e6 = half*epspzx(i)
333 epsp =half*(e1**2+e2**2+e3**2) +e4**2+e5**2+e6**2
334 epsp = sqrt(three*epsp)/three_half
335 syield(i) = fac1*syield(i)*(finter(ifunc(2),epsp,npf,tf,b*c))
336 ENDIF
337 210 CONTINUE
338
339 DO 220 i=1,nel
340 sigsxx(i)=sigoxx(i)+sigair(i)
341 sigsyy(i)=sigoyy(i)+sigair(i)
342 sigszz(i)=sigozz(i)+sigair(i)
343 sigsxy(i)=sigoxy(i)
344 sigsyz(i)=sigoyz(i)
345 sigszx(i)=sigozx(i)
346
347 dexx(i)=epsxx(i)
348 deyy(i)=epsyy(i)
349 dezz(i)=epszz(i)
350 dexy(i)=epsxy(i)* half
351 deyz(i)=epsyz(i)* half
352 dezx(i)=epszx(i)* half
353
354 dedtxx(i)=epspxx(i)
355 dedtyy(i)=epspyy(i)
356 dedtzz(i)=epspzz(i)
357 dedtxy(i)=epspxy(i) * half
358 dedtyz(i)=epspyz(i) * half
359 dedtzx(i)=epspzx(i) * half
360 220 CONTINUE
361C COMPUTE STRESS RATES
362 DO 230 i=1,nel
363 dsdtxx(i)=e(i)*dedtxx(i)-epet(i)*sigsxx(i)+emet(i)*dexx(i)
364 dsdtyy(i)=e(i)*dedtyy(i)-epet(i)*sigsyy(i)+emet(i)*deyy(i)
365 dsdtzz(i)=e(i)*dedtzz(i)-epet(i)*sigszz(i)+emet(i)*dezz(i)
366 dsdtxy(i)=e(i)*dedtxy(i)-epets(i)*sigsxy(i)+emets(i)*dexy(i)
367 dsdtyz(i)=e(i)*dedtyz(i)-epets(i)*sigsyz(i)+emets(i)*deyz(i)
368 dsdtzx(i)=e(i)*dedtzx(i)-epets(i)*sigszx(i)+emets(i)*dezx(i)
369 230 CONTINUE
370
371C COMPUTE TRIAL STRESSES
372 DO 240 i=1,nel
373 sigsxx(i)=sigsxx(i)+dsdtxx(i)*timestep
374 sigsyy(i)=sigsyy(i)+dsdtyy(i)*timestep
375 sigszz(i)=sigszz(i)+dsdtzz(i)*timestep
376 sigsxy(i)=sigsxy(i)+dsdtxy(i)*timestep
377 sigsyz(i)=sigsyz(i)+dsdtyz(i)*timestep
378 sigszx(i)=sigszx(i)+dsdtzx(i)*timestep
379 240 CONTINUE
380 IF(ken<0)GOTO 255
381C
382C If old method then ....
383 IF (vpvec==1) THEN
384 DO 250 i=1,nel
385 s(1,1)=sigsxx(i)
386 s(2,1)=sigsxy(i)
387 s(2,2)=sigsyy(i)
388 s(3,1)=sigszx(i)
389 s(3,2)=sigsyz(i)
390 s(3,3)=sigszz(i)
391
392C TRANSFORM GLOBAL STRESSES (S_IJ) TO PRINCIPAL STRESSES (SIGPR_K)
393 CALL jacobiew(s,3,sigpr,dirpr,nrot)
394
395C YIELD CRITERIA - SCALING
396 sigpr(1)=min(syield(i),abs(sigpr(1)))*sign(one,sigpr(1))
397 sigpr(2)=min(syield(i),abs(sigpr(2)))*sign(one,sigpr(2))
398 sigpr(3)=min(syield(i),abs(sigpr(3)))*sign(one,sigpr(3))
399
400C TRANSFORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
401 DO 251 k=1,3
402 DO 251 l=1,3
403 s(k,l)=zero
404 IF(k==l) s(k,l)= sigpr(k)
405 251 CONTINUE
406 CALL dreh(s,dirpr,1,1,1)
407
408 sigsxx(i)=s(1,1)
409 sigsxy(i)=s(2,1)
410 sigsyy(i)=s(2,2)
411 sigszx(i)=s(3,1)
412 sigsyz(i)=s(3,2)
413 sigszz(i)=s(3,3)
414 250 CONTINUE
415 ELSE
416C fin ancienne methode
417 DO i = 1, nel
418 sigv(i,1) = sigsxx(i)
419 sigv(i,2) = sigsyy(i)
420 sigv(i,3) = sigszz(i)
421 sigv(i,4) = sigsxy(i)
422 sigv(i,5) = sigsyz(i)
423 sigv(i,6) = sigszx(i)
424 ENDDO
425 CALL valpvec_v(sigv,sigprv,dirprv,nel)
426 DO i = 1, nel
427C YIELD CRITERIA - SCALING
428 sigprv(i,1)=sign(min(syield(i),abs(sigprv(i,1))),sigprv(i,1))
429 sigprv(i,2)=sign(min(syield(i),abs(sigprv(i,2))),sigprv(i,2))
430 sigprv(i,3)=sign(min(syield(i),abs(sigprv(i,3))),sigprv(i,3))
431C TRANSFORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
432 sigsxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigprv(i,1)
433 . + dirprv(i,1,2)*dirprv(i,1,2)*sigprv(i,2)
434 . + dirprv(i,1,3)*dirprv(i,1,3)*sigprv(i,3)
435 sigsyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigprv(i,2)
436 . + dirprv(i,2,3)*dirprv(i,2,3)*sigprv(i,3)
437 . + dirprv(i,2,1)*dirprv(i,2,1)*sigprv(i,1)
438 sigszz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigprv(i,3)
439 . + dirprv(i,3,1)*dirprv(i,3,1)*sigprv(i,1)
440 . + dirprv(i,3,2)*dirprv(i,3,2)*sigprv(i,2)
441 sigsxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigprv(i,1)
442 . + dirprv(i,1,2)*dirprv(i,2,2)*sigprv(i,2)
443 . + dirprv(i,1,3)*dirprv(i,2,3)*sigprv(i,3)
444 sigsyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigprv(i,2)
445 . + dirprv(i,2,3)*dirprv(i,3,3)*sigprv(i,3)
446 . + dirprv(i,2,1)*dirprv(i,3,1)*sigprv(i,1)
447 sigszx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigprv(i,3)
448 . + dirprv(i,3,1)*dirprv(i,1,1)*sigprv(i,1)
449 . + dirprv(i,3,2)*dirprv(i,1,2)*sigprv(i,2)
450 ENDDO
451 ENDIF
452C fin chgt
453 255 CONTINUE
454C COMPUTE UPDATED STRESSES AND SOUND SPEED
455 DO 260 i=1,nel
456 signxx(i)=off(i)*sigsxx(i)-sigair(i)
457 signyy(i)=off(i)*sigsyy(i)-sigair(i)
458 signzz(i)=off(i)*sigszz(i)-sigair(i)
459 signxy(i)=off(i)*sigsxy(i)
460 signyz(i)=off(i)*sigsyz(i)
461 signzx(i)=off(i)*sigszx(i)
462
463 soundsp(i)=sqrt(e(i)/rho0(i))
464 260 CONTINUE
465C--------------------------
466 CASE(3)
467C--------------------------
468 fac = uparam(9)
469 fac1 = uparam(10)
470 sigt_cutoff = uparam(11)
471 sigc_cutoff = sigt_cutoff
472 DO i=1,nel
473 IF(ifunc(1)/=0) THEN
474 syield(i)=fac*finter(ifunc(1),gama(i),npf,tf,b*c)
475 ELSE
476 syield(i) = abs(a+b*(one+c*gama(i)))
477 ENDIF
478 IF(gama(i) > zero) THEN
479 IF (gama(i)-uvar(i,1) < zero) sigc_cutoff = zero
480 ENDIF
481 IF(ifunc(2)/=0)THEN
482 dav = (epspxx(i)+epspyy(i)+epspzz(i))*third
483 e1 = epspxx(i) - dav
484 e2 = epspyy(i) - dav
485 e3 = epspzz(i) - dav
486 e4 = half*epspxy(i)
487 e5 = half*epspyz(i)
488 e6 = half*epspzx(i)
489 epsp =half*(e1**2+e2**2+e3**2) +e4**2+e5**2+e6**2
490 epsp = sqrt(three*epsp)/three_half
491 syield(i) = fac1*syield(i)*(finter(ifunc(2),epsp,npf,tf,b*c))
492 ENDIF
493 ENDDO
494
495 DO i=1,nel
496 sigsxx(i)=sigoxx(i)+sigair(i)
497 sigsyy(i)=sigoyy(i)+sigair(i)
498 sigszz(i)=sigozz(i)+sigair(i)
499 sigsxy(i)=sigoxy(i)
500 sigsyz(i)=sigoyz(i)
501 sigszx(i)=sigozx(i)
502 ENDDO
503
504C COMPUTE TRIAL STRESSES
505 DO i=1,nel
506 sigsxx(i)=sigsxx(i)+ey*epspxx(i)*timestep
507 sigsyy(i)=sigsyy(i)+ey*epspyy(i)*timestep
508 sigszz(i)=sigszz(i)+ey*epspzz(i)*timestep
509 sigsxy(i)=sigsxy(i)+ey*epspxy(i)*timestep*half
510 sigsyz(i)=sigsyz(i)+ey*epspyz(i)*timestep*half
511 sigszx(i)=sigszx(i)+ey*epspzx(i)*timestep*half
512 ENDDO
513 DO i = 1, nel
514 sigv(i,1) = sigsxx(i)
515 sigv(i,2) = sigsyy(i)
516 sigv(i,3) = sigszz(i)
517 sigv(i,4) = sigsxy(i)
518 sigv(i,5) = sigsyz(i)
519 sigv(i,6) = sigszx(i)
520 ENDDO
521C for a simple precision executing
522 IF (iresp==1) THEN
523 CALL valpvecdp_v(sigv,sigprv,dirprv,nel)
524 ELSE
525 CALL valpvec_v(sigv,sigprv,dirprv,nel)
526 ENDIF
527 DO i =1, nel
528C YIELD CRITERIA - SCALING
529 IF(abs(gama(i)) < em10) THEN
530 sigprv(i,1)=sign(min(abs(sigprv(i,1)), sigt_cutoff),sigprv(i,1))
531 sigprv(i,2)=sign(min(abs(sigprv(i,2)), sigt_cutoff),sigprv(i,2))
532 sigprv(i,3)=sign(min(abs(sigprv(i,3)), sigt_cutoff),sigprv(i,3))
533 ELSEIF(gama(i) < zero )then
534 !!IF(DELTA_GAMA(I) < ZERO)then
535 sigprv(i,1)=sign(min(syield(i),abs(sigprv(i,1))),sigprv(i,1))
536 sigprv(i,2)=sign(min(syield(i),abs(sigprv(i,2))),sigprv(i,2))
537 sigprv(i,3)=sign(min(syield(i),abs(sigprv(i,3))),sigprv(i,3))
538 !! ELSE
539 sigprv(i,1)=min(sigprv(i,1), sigt_cutoff)
540 sigprv(i,2)=min(sigprv(i,2), sigt_cutoff)
541 sigprv(i,3)=min(sigprv(i,3), sigt_cutoff)
542 !! ENDIF
543 ELSE
544 sigprv(i,1)=sign(min(abs(sigprv(i,1)), sigt_cutoff),sigprv(i,1))
545 sigprv(i,2)=sign(min(abs(sigprv(i,2)), sigt_cutoff),sigprv(i,2))
546 sigprv(i,3)=sign(min(abs(sigprv(i,3)), sigt_cutoff),sigprv(i,3))
547 !!SIGPRV(I,1)=MIN(SIGPRV(I,1), SIGT_CUTOFF)
548 !!SIGPRV(I,2)=MIN(SIGPRV(I,2), SIGT_CUTOFF)
549 !!SIGPRV(I,3)=MIN(SIGPRV(I,3), SIGT_CUTOFF)
550 !! IF(DELTA_GAMA(I) < ZERO) THEN
551 sigprv(i,1)=max(sigprv(i,1), zero)
552 sigprv(i,2)=max(sigprv(i,2), zero)
553 sigprv(i,3)=max(sigprv(i,3), zero)
554 !! ENDIF
555 ENDIF
556C TRANSF ORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
557 sigsxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigprv(i,1)
558 . + dirprv(i,1,2)*dirprv(i,1,2)*sigprv(i,2)
559 . + dirprv(i,1,3)*dirprv(i,1,3)*sigprv(i,3)
560!!
561 sigsyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigprv(i,2)
562 . + dirprv(i,2,3)*dirprv(i,2,3)*sigprv(i,3)
563 . + dirprv(i,2,1)*dirprv(i,2,1)*sigprv(i,1)
564!!
565 sigszz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigprv(i,3)
566 . + dirprv(i,3,1)*dirprv(i,3,1)*sigprv(i,1)
567 . + dirprv(i,3,2)*dirprv(i,3,2)*sigprv(i,2)
568!!
569 sigsxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigprv(i,1)
570 . + dirprv(i,1,2)*dirprv(i,2,2)*sigprv(i,2)
571 . + dirprv(i,1,3)*dirprv(i,2,3)*sigprv(i,3)
572!!
573 sigsyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigprv(i,2)
574 . + dirprv(i,2,3)*dirprv(i,3,3)*sigprv(i,3)
575 . + dirprv(i,2,1)*dirprv(i,3,1)*sigprv(i,1)
576!!
577 sigszx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigprv(i,3)
578 . + dirprv(i,3,1)*dirprv(i,1,1)*sigprv(i,1)
579 . + dirprv(i,3,2)*dirprv(i,1,2)*sigprv(i,2)
580!!
581 ENDDO
582C fin chgt
583C COMPUTE UPDATED STRESSES AND SOUND SPEED
584C
585 DO i=1,nel
586 signxx(i)=off(i)*sigsxx(i)-sigair(i)
587 signyy(i)=off(i)*sigsyy(i)-sigair(i)
588 signzz(i)=off(i)*sigszz(i)-sigair(i)
589 signxy(i)=off(i)*sigsxy(i)
590 signyz(i)=off(i)*sigsyz(i)
591 signzx(i)=off(i)*sigszx(i)
592 soundsp(i)=sqrt(ey/rho0(i))
593 ENDDO
594 END SELECT
595
596 RETURN
597
598 END
599
600!||====================================================================
601!|| jacobiew ../engine/source/materials/mat/mat033/sigeps33.F
602!||--- called by ------------------------------------------------------
603!|| fail_orthstrain ../engine/source/materials/fail/orthstrain/fail_orthstrain_s.F
604!|| hencky_strain ../engine/source/tools/seatbelts/shell_reactivation.F
605!|| sigeps33 ../engine/source/materials/mat/mat033/sigeps33.F
606!|| sigeps38 ../engine/source/materials/mat/mat038/sigeps38.F
607!||====================================================================
608 SUBROUTINE jacobiew(A,N,EW,EV,NROT)
609C-------------------------------------------------------KK141189-
610C COMPUTATION OF ALL EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC
611C MATRIX A BY THE JACOBI ALGORITHM
612C
613C A(N,N) EIGENWERTPROBLEM
614C N DIMENSION OF A
615C EW(N) EIGENVALUES
616C EV(N,N) EIGENVEKTORS
617C NROT NUMBER OF ROTATIONS
618C MAXA MAXIMUM ELEMENT OF A
619C-----------------------------------------------
620C I m p l i c i t T y p e s
621C-----------------------------------------------
622#include "implicit_f.inc"
623 INTEGER NN
624 PARAMETER (NN=9)
625
626 integer n,nrot
627 my_real
628 . a(n,n), ew(n), ev(n,n)
629 . , b(nn), z(nn)
630 INTEGER IZ,IS,ITER,J
631 my_real
632 . SUMRS,EPS,G,H,T,C,S,TAU,THETA
633C----------------------------------------------------------------
634 DO 130 iz=1,n
635 DO 120 is=1,n
636C PRACTICAL FOR RADIOSS (PASS ONLY LOWER DIAGONAL MATRIX)
637 IF(iz>is) a(is,iz) = a(iz,is)
638 ev(iz,is)=zero
639 120 CONTINUE
640 b(iz)=a(iz,iz)
641 ew(iz)=b(iz)
642 z(iz)=0.
643 ev(iz,iz)=one
644 130 CONTINUE
645
646 nrot=0
647
648C START ITERATION
649
650 DO 240 iter = 1,50
651 sumrs = zero
652
653C SUM OF THE OFF DIAGONALS
654 DO 150 iz=1,n-1
655 DO 140 is=iz+1,n
656 sumrs=sumrs+abs(a(iz,is))
657 140 CONTINUE
658 150 CONTINUE
659
660 IF (sumrs ==zero) GOTO 9000
661 IF (iter > 4) THEN
662 eps = zero
663 ELSE
664 eps = one_fifth*sumrs/n**2
665 ENDIF
666
667 DO 220 iz=1,n-1
668 DO 210 is=iz+1,n
669 g = 100. * abs(a(iz,is))
670 IF (iter>4 .AND. abs(ew(iz))+g==abs(ew(iz))
671 & .AND. abs(ew(is))+g==abs(ew(is))) THEN
672 a(iz,is)=zero
673 ELSE IF (abs(a(iz,is)) > eps) THEN
674 h = ew(is)-ew(iz)
675 IF (abs(h)+g==abs(h)) THEN
676 t = a(iz,is)/h
677 ELSE
678 theta = half*h/a(iz,is)
679 t=one/(abs(theta)+sqrt(one+theta**2))
680 IF (theta < zero) t=-t
681 ENDIF
682 c=one/sqrt(one+t**2)
683 s=t*c
684 tau=s/(one+c)
685 h=t*a(iz,is)
686 z(iz)=z(iz)-h
687 z(is)=z(is)+h
688 ew(iz)=ew(iz)-h
689 ew(is)=ew(is)+h
690 a(iz,is)=zero
691 DO 160 j=1,iz-1
692 g=a(j,iz)
693 h=a(j,is)
694 a(j,iz)=g-s*(h+g*tau)
695 a(j,is)=h+s*(g-h*tau)
696 160 CONTINUE
697 DO 170 j=iz+1,is-1
698 g=a(iz,j)
699 h=a(j,is)
700 a(iz,j)=g-s*(h+g*tau)
701 a(j,is)=h+s*(g-h*tau)
702 170 CONTINUE
703 DO 180 j=is+1,n
704 g=a(iz,j)
705 h=a(is,j)
706 a(iz,j)=g-s*(h+g*tau)
707 a(is,j)=h+s*(g-h*tau)
708 180 CONTINUE
709 DO 190 j=1,n
710 g=ev(j,iz)
711 h=ev(j,is)
712 ev(j,iz)=g-s*(h+g*tau)
713 ev(j,is)=h+s*(g-h*tau)
714 190 CONTINUE
715 nrot=nrot+1
716 ENDIF
717 210 CONTINUE
718 220 CONTINUE
719
720 DO 230 iz=1,n
721 b(iz)=b(iz)+z(iz)
722 ew(iz)=b(iz)
723 z(iz)=zero
724 230 CONTINUE
725 240 CONTINUE
726
727 9000 CONTINUE
728
729 RETURN
730
731 END
732
733
734!||====================================================================
735!|| dreh ../engine/source/materials/mat/mat033/sigeps33.F
736!||--- called by ------------------------------------------------------
737!|| sigeps33 ../engine/source/materials/mat/mat033/sigeps33.F
738!|| sigeps38 ../engine/source/materials/mat/mat038/sigeps38.F
739!||====================================================================
740 SUBROUTINE dreh(B,TR,II,JJ,KEN)
741C-------------------------------------------------------KK010587-
742C KEN=0 -> ROTATION OF MATRIX B BY TR = TR(T)*B*TR
743C KEN=1 -> ROTATION OF MATRIX B BY TR(T) = TR*B*TR(T)
744C II AND JJ POINT TO THE SUBMATRIX TO BE ROTATED; DEFAULT II=JJ=1
745C----------------------------------------------------------------
746C-----------------------------------------------
747C I m p l i c i t T y p e s
748C-----------------------------------------------
749#include "implicit_f.inc"
750 my_real
751 . b(3,3),tr(3,3),lk(3,3),x
752 INTEGER KEN,II,JJ
753 INTEGER I,J,K,I1,J1
754C----------------------------------------------------------------
755 IF(II<=0) ii=1
756 IF(jj<=0) jj=1
757
758 DO 20 i=1,3
759 DO 20 j=1,3
760 j1=jj+j-1
761 x=0.0
762 DO 10 k=1,3
763C IF(TR(K,I)==ZERO) GOTO 10
764 i1=k+ii-1
765C IF(B(I1,J1)==ZERO) GOTO 10
766 IF(ken/=1) x=x+tr(k,i)*b(i1,j1)
767 IF(ken==1) x=x+tr(i,k)*b(i1,j1)
768 10 CONTINUE
769 20 lk(i,j)=x
770
771 DO 40 i=1,3
772 DO 40 j=1,3
773 x=zero
774 DO 30 k=1,3
775C IF(TR(K,J)==ZERO) GOTO 30
776C IF(LK(I,K)==ZERO) GOTO 30
777 IF(ken/=1) x=x+lk(i,k)*tr(k,j)
778 IF(ken==1) x=x+lk(i,k)*tr(j,k)
779 30 CONTINUE
780 40 b(ii+i-1,jj+j-1)=x
781
782 RETURN
783
784 END
785C
786!||====================================================================
787!|| valpvec ../engine/source/materials/mat/mat033/sigeps33.F
788!||--- called by ------------------------------------------------------
789!|| alew7 ../engine/source/ale/grid/alew7.F
790!|| epsf2u ../engine/source/materials/mat/mat033/sigeps33.F
791!||--- calls -----------------------------------------------------
792!|| floatmin ../common_source/tools/math/precision.c
793!||====================================================================
794 SUBROUTINE valpvec(SIG,VAL,VEC,NEL)
795C-----------------------------------------------
796C I m p l i c i t T y p e s
797C-----------------------------------------------
798#include "implicit_f.inc"
799C-----------------------------------------------
800C G l o b a l P a r a m e t e r s
801C-----------------------------------------------
802#include "mvsiz_p.inc"
803C-----------------------------------------------
804C D u m m y A r g u m e n t s
805C-----------------------------------------------
806 my_real
807 . sig(6,*), val(3,*), vec(9,*)
808 INTEGER NEL
809C-----------------------------------------------
810C L o c a l V a r i a b l e s
811C-----------------------------------------------
812 INTEGER N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
813 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
814 my_real
815 . CS(6), STR(3,MVSIZ), A(3,3,MVSIZ), V(3,3,MVSIZ), B(3,3,MVSIZ),
816 . xmag(3,mvsiz), pr, aa, bb, aaa(mvsiz),
817 . cc, angp, dd, ftpi, ttpi, strmax,
818 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
819 . vmag,
820 .
821 .
822 . mdemi, xmaxinv, flm
823 REAL FLMIN
824C-----------------------------------------------
825C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
826C FTPI=(4/3)*PI, TTPI=(2/3)*PI
827C
828C principal deviator of stress
829C . . . . . . . . . . . . . . . . . . .
830 mdemi = -half
831 ttpi = acos(mdemi)
832 ftpi = two*ttpi
833C minimum precision depending on the type real or double
834 CALL floatmin(cs(1),cs(2),flmin)
835 flm = two*sqrt(flmin)
836 nindex3=0
837 DO nn = 1, nel
838 cs(1) = sig(1,nn)
839 cs(2) = sig(2,nn)
840 cs(3) = sig(3,nn)
841 cs(4) = sig(4,nn)
842 cs(5) = sig(5,nn)
843 cs(6) = sig(6,nn)
844 pr = -(cs(1)+cs(2)+cs(3))* third
845 cs(1) = cs(1) + pr
846 cs(2) = cs(2) + pr
847 cs(3) = cs(3) + pr
848 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
849 & - cs(2)*cs(3) - cs(1)*cs(3)
850 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
851 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
852 norminf(nn) = em10*norminf(nn)
853C cas racine triple
854c AA = MAX(AAA(NN),NORMINF(NN),EM20)
855 aa = max(aaa(nn),norminf(nn))
856C
857 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
858 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
859 & - two*cs(4)*cs(5)*cs(6)
860C
861 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
862 cc= min(cc,one)
863 cc= max(cc,-one)
864 angp=acos(cc) * third
865 dd=two*sqrt(aa * third)
866 str(1,nn)=dd*cos(angp)
867 str(2,nn)=dd*cos(angp+ftpi)
868 str(3,nn)=dd*cos(angp+ttpi)
869C
870 val(1,nn) = str(1,nn)-pr
871 val(2,nn) = str(2,nn)-pr
872 val(3,nn) = str(3,nn)-pr
873C accuracy reinforcement in compression ----
874 IF(abs(str(3,nn))>abs(str(1,nn))
875 & .AND.aaa(nn)>norminf(nn)) THEN
876 aa=str(1,nn)
877 str(1,nn)=str(3,nn)
878 str(3,nn)=aa
879 nindex3 = nindex3+1
880 index3(nindex3) = nn
881 ENDIF
882C . . . . . . . . . . .
883C VECTEURS PROPRES
884C . . . . . . . . . . .
885 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
886 tol1(nn)= max(em20,flm*strmax**2)
887 tol2(nn)=flm*strmax/3
888 a(1,1,nn)=cs(1)-str(1,nn)
889 a(2,2,nn)=cs(2)-str(1,nn)
890 a(3,3,nn)=cs(3)-str(1,nn)
891 a(1,2,nn)=cs(4)
892 a(2,1,nn)=cs(4)
893 a(2,3,nn)=cs(5)
894 a(3,2,nn)=cs(5)
895 a(1,3,nn)=cs(6)
896 a(3,1,nn)=cs(6)
897C
898 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
899 . *a(2,2,nn)
900 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
901 . *a(2,3,nn)
902 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
903 . *a(2,1,nn)
904 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
905 . *a(3,2,nn)
906 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
907 . *a(3,3,nn)
908 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
909 . *a(3,1,nn)
910 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
911 . *a(1,2,nn)
912 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
913 . *a(1,3,nn)
914 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
915 . *a(1,1,nn)
916 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
917 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
918 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
919
920 ENDDO
921C
922 nindex1 = 0
923 nindex2 = 0
924 DO nn = 1, nel
925 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
926 IF(xmag(1,nn)==xmaxv(nn)) THEN
927 lmaxv(nn) = 1
928 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
929 lmaxv(nn) = 2
930 ELSE
931 lmaxv(nn) = 3
932 ENDIF
933 IF(aaa(nn)<norminf(nn)) THEN
934 val(1,nn) = sig(1,nn)
935 val(2,nn) = sig(2,nn)
936 val(3,nn) = sig(3,nn)
937 v(1,1,nn) = one
938 v(2,1,nn) = zero
939 v(3,1,nn) = zero
940 v(1,2,nn) = zero
941 v(2,2,nn) = one
942 v(3,2,nn) = zero
943
944 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
945 nindex1 = nindex1 + 1
946 index1(nindex1) = nn
947 ELSE
948 nindex2 = nindex2 + 1
949 index2(nindex2) = nn
950 ENDIF
951 ENDDO
952C
953#include "vectorize.inc"
954 DO n = 1, nindex1
955 nn = index1(n)
956 lmax = lmaxv(nn)
957 xmaxinv = one/xmaxv(nn)
958 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
959 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
960 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
961 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
962 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
963 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
964C
965 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
966 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
967 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
968 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
969 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
970 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
971 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
972 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
973 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
974 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
975 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
976 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
977C
978 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
979 ENDDO
980C
981#include "vectorize.inc"
982 DO n = 1, nindex1
983 nn = index1(n)
984 IF(xmag(1,nn)==xmaxv(nn)) THEN
985 lmaxv(nn) = 1
986 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
987 lmaxv(nn) = 2
988 ELSE
989 lmaxv(nn) = 3
990 ENDIF
991C
992 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
993 lmax = lmaxv(nn)
994 IF(xmaxv(nn)>tol2(nn))THEN
995 xmaxinv = one/xmaxv(nn)
996 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
997 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
998 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
999 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
1000 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
1001 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
1002 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
1003 v(1,2,nn)=v(1,2,nn)*vmag
1004 v(2,2,nn)=v(2,2,nn)*vmag
1005 v(3,2,nn)=v(3,2,nn)*vmag
1006 ELSEIF(vmag>tol2(nn))THEN
1007 v(1,2,nn)=-v(2,1,nn)/vmag
1008 v(2,2,nn)=v(1,1,nn)/vmag
1009 v(3,2,nn)=zero
1010 ELSE
1011 v(1,2,nn)=one
1012 v(2,2,nn)=zero
1013 v(3,2,nn)=zero
1014 ENDIF
1015 ENDDO
1016C . . . . . . . . . . . . .
1017C SOLUTION DOUBLE
1018C . . . . . . . . . . . . .
1019 DO n = 1, nindex2
1020 nn = index2(n)
1021 xmag(1,nn)=sqrt(a(1,1,nn)**2+a(2,1,nn)**2)
1022 xmag(2,nn)=sqrt(a(1,2,nn)**2+a(2,2,nn)**2)
1023 xmag(3,nn)=sqrt(a(1,3,nn)**2+a(2,3,nn)**2)
1024C
1025 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1026 ENDDO
1027C
1028#include "vectorize.inc"
1029 DO n = 1, nindex2
1030 nn = index2(n)
1031 IF(xmag(1,nn)==xmaxv(nn)) THEN
1032 lmaxv(nn) = 1
1033 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1034 lmaxv(nn) = 2
1035 ELSE
1036 lmaxv(nn) = 3
1037 ENDIF
1038C
1039 lmax = lmaxv(nn)
1040 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
1041 & <tol2(nn))THEN
1042 xmaxinv = one/xmaxv(nn)
1043 v(1,1,nn)= zero
1044 v(2,1,nn)= zero
1045 v(3,1,nn)= one
1046 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
1047 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
1048 v(3,2,nn)= zero
1049C
1050 ELSEIF(xmaxv(nn)>tol2(nn))THEN
1051 xmaxinv = one/xmaxv(nn)
1052 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
1053 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
1054 v(3,1,nn)= zero
1055 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
1056 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
1057 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
1058 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
1059 v(1,2,nn)=v(1,2,nn)*vmag
1060 v(2,2,nn)=v(2,2,nn)*vmag
1061 v(3,2,nn)=v(3,2,nn)*vmag
1062 ELSE
1063 v(1,1,nn) = one
1064 v(2,1,nn) = zero
1065 v(3,1,nn) = zero
1066 v(1,2,nn) = zero
1067 v(2,2,nn) = one
1068 v(3,2,nn) = zero
1069 ENDIF
1070 ENDDO
1071C
1072 DO nn = 1, nel
1073 vec(1,nn)=v(1,1,nn)
1074 vec(2,nn)=v(2,1,nn)
1075 vec(3,nn)=v(3,1,nn)
1076 vec(4,nn)=v(1,2,nn)
1077 vec(5,nn)=v(2,2,nn)
1078 vec(6,nn)=v(3,2,nn)
1079 vec(7,nn)=vec(2,nn)*vec(6,nn)-vec(3,nn)*vec(5,nn)
1080 vec(8,nn)=vec(3,nn)*vec(4,nn)-vec(1,nn)*vec(6,nn)
1081 vec(9,nn)=vec(1,nn)*vec(5,nn)-vec(2,nn)*vec(4,nn)
1082 ENDDO
1083C Rewriting to bypass problem on itanium2 comp 9. + latency = 16
1084 DO n = 1, nindex3
1085 nn = index3(n)
1086C Str uses temporary table instead of temporary scalars
1087 str(1,nn)=vec(7,nn)
1088 str(2,nn)=vec(8,nn)
1089 str(3,nn)=vec(9,nn)
1090 ENDDO
1091 DO n = 1, nindex3
1092 nn = index3(n)
1093 vec(7,nn)=vec(1,nn)
1094 vec(8,nn)=vec(2,nn)
1095 vec(9,nn)=vec(3,nn)
1096 vec(1,nn)=-str(1,nn)
1097 vec(2,nn)=-str(2,nn)
1098 vec(3,nn)=-str(3,nn)
1099 ENDDO
1100C
1101 RETURN
1102 END
1103
1104!||====================================================================
1105!|| valpvecdp ../engine/source/materials/mat/mat033/sigeps33.F
1106!||--- called by ------------------------------------------------------
1107!|| epsf2u ../engine/source/materials/mat/mat033/sigeps33.F
1108!||--- calls -----------------------------------------------------
1109!|| floatmin ../common_source/tools/math/precision.c
1110!||====================================================================
1111 SUBROUTINE valpvecdp(SIG,VAL,VEC,NEL)
1112C-----------------------------------------------
1113C I m p l i c i t T y p e s
1114C-----------------------------------------------
1115#include "implicit_f.inc"
1116C-----------------------------------------------
1117C G l o b a l P a r a m e t e r s
1118C-----------------------------------------------
1119#include "mvsiz_p.inc"
1120C-----------------------------------------------
1121C D u m m y A r g u m e n t s
1122C-----------------------------------------------
1123 my_real
1124 . sig(6,*), val(3,*), vec(9,*)
1125 INTEGER NEL
1126C-----------------------------------------------
1127C L o c a l V a r i a b l e s
1128C-----------------------------------------------
1129 INTEGER N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
1130 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
1131 DOUBLE PRECISION
1132 . CS(6), STR(3,MVSIZ), A(3,3,MVSIZ), V(3,3,MVSIZ), B(3,3,MVSIZ),
1133 . xmag(3,mvsiz), pr, aa, bb, aaa(mvsiz),
1134 . cc, angp, dd, ftpi, ttpi, strmax,
1135 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
1136 . vmag,
1137 .
1138 .
1139 . mdemi, xmaxinv, flm,
1140 . valdp(3,mvsiz),vecdp(9,mvsiz)
1141 REAL FLMIN
1142C-----------------------------------------------
1143C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
1144C FTPI=(4/3)*PI, TTPI=(2/3)*PI
1145C
1146C principal deviator of stress
1147C . . . . . . . . . . . . . . . . . . .
1148 MDEMI = -half
1149 ttpi = acos(mdemi)
1150 ftpi = two*ttpi
1151C minimum precision depending on the type real or double
1152 CALL floatmin(cs(1),cs(2),flmin)
1153 flm = two*sqrt(flmin)
1154 nindex3=0
1155 DO nn = 1, nel
1156 cs(1) = sig(1,nn)
1157 cs(2) = sig(2,nn)
1158 cs(3) = sig(3,nn)
1159 cs(4) = sig(4,nn)
1160 cs(5) = sig(5,nn)
1161 cs(6) = sig(6,nn)
1162 pr = -(cs(1)+cs(2)+cs(3)) * third
1163 cs(1) = cs(1) + pr
1164 cs(2) = cs(2) + pr
1165 cs(3) = cs(3) + pr
1166 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
1167 & - cs(2)*cs(3) - cs(1)*cs(3)
1168 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
1169 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
1170 norminf(nn) = em10*norminf(nn)
1171C cas racine triple
1172c AA = MAX(AAA(NN),NORMINF(NN),EM20)
1173 aa = max(aaa(nn),norminf(nn))
1174C
1175 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
1176 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
1177 & - two*cs(4)*cs(5)*cs(6)
1178C
1179 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
1180 cc= min(cc,one)
1181 cc= max(cc,-one)
1182 angp=acos(cc) * third
1183 dd=two*sqrt(aa * third)
1184 str(1,nn)=dd*cos(angp)
1185 str(2,nn)=dd*cos(angp+ftpi)
1186 str(3,nn)=dd*cos(angp+ttpi)
1187C
1188 valdp(1,nn) = str(1,nn)-pr
1189 valdp(2,nn) = str(2,nn)-pr
1190 valdp(3,nn) = str(3,nn)-pr
1191C precision reinforcement in simple compression ----
1192 IF(abs(str(3,nn))>abs(str(1,nn))
1193 & .AND.aaa(nn)>norminf(nn)) THEN
1194 aa=str(1,nn)
1195 str(1,nn)=str(3,nn)
1196 str(3,nn)=aa
1197 nindex3 = nindex3+1
1198 index3(nindex3) = nn
1199 ENDIF
1200C . . . . . . . . . . .
1201C VECTEURS PROPRES
1202C . . . . . . . . . . .
1203 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
1204 tol1(nn)= max(em20,flm*strmax**2)
1205 tol2(nn)=flm*strmax * third
1206 a(1,1,nn)=cs(1)-str(1,nn)
1207 a(2,2,nn)=cs(2)-str(1,nn)
1208 a(3,3,nn)=cs(3)-str(1,nn)
1209 a(1,2,nn)=cs(4)
1210 a(2,1,nn)=cs(4)
1211 a(2,3,nn)=cs(5)
1212 a(3,2,nn)=cs(5)
1213 a(1,3,nn)=cs(6)
1214 a(3,1,nn)=cs(6)
1215C
1216 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
1217 . *a(2,2,nn)
1218 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
1219 . *a(2,3,nn)
1220 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
1221 . *a(2,1,nn)
1222 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
1223 . *a(3,2,nn)
1224 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
1225 . *a(3,3,nn)
1226 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
1227 . *a(3,1,nn)
1228 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
1229 . *a(1,2,nn)
1230 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
1231 . *a(1,3,nn)
1232 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
1233 . *a(1,1,nn)
1234 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
1235 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
1236 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
1237
1238 ENDDO
1239C
1240 nindex1 = 0
1241 nindex2 = 0
1242 DO nn = 1, nel
1243 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1244 IF(xmag(1,nn)==xmaxv(nn)) THEN
1245 lmaxv(nn) = 1
1246 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1247 lmaxv(nn) = 2
1248 ELSE
1249 lmaxv(nn) = 3
1250 ENDIF
1251 IF(aaa(nn)<norminf(nn)) THEN
1252 valdp(1,nn) = sig(1,nn)
1253 valdp(2,nn) = sig(2,nn)
1254 valdp(3,nn) = sig(3,nn)
1255 v(1,1,nn) = one
1256 v(2,1,nn) = zero
1257 v(3,1,nn) = zero
1258 v(1,2,nn) = zero
1259 v(2,2,nn) = one
1260 v(3,2,nn) = zero
1261 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
1262 nindex1 = nindex1 + 1
1263 index1(nindex1) = nn
1264 ELSE
1265 nindex2 = nindex2 + 1
1266 index2(nindex2) = nn
1267 ENDIF
1268 ENDDO
1269C
1270#include "vectorize.inc"
1271 DO n = 1, nindex1
1272 nn = index1(n)
1273 lmax = lmaxv(nn)
1274 xmaxinv = one/xmaxv(nn)
1275 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
1276 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
1277 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
1278 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
1279 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
1280 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
1281C
1282 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
1283 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
1284 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
1285 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
1286 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
1287 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
1288 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
1289 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
1290 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
1291 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
1292 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
1293 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
1294C
1295 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1296 ENDDO
1297C
1298#include "vectorize.inc"
1299 DO n = 1, nindex1
1300 nn = index1(n)
1301 IF(xmag(1,nn)==xmaxv(nn)) THEN
1302 lmaxv(nn) = 1
1303 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1304 lmaxv(nn) = 2
1305 ELSE
1306 lmaxv(nn) = 3
1307 ENDIF
1308C
1309 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
1310 lmax = lmaxv(nn)
1311 IF(xmaxv(nn)>tol2(nn))THEN
1312 xmaxinv = one/xmaxv(nn)
1313 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
1314 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
1315 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
1316 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
1317 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
1318 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
1319 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
1320 v(1,2,nn)=v(1,2,nn)*vmag
1321 v(2,2,nn)=v(2,2,nn)*vmag
1322 v(3,2,nn)=v(3,2,nn)*vmag
1323 ELSEIF(vmag>tol2(nn))THEN
1324 v(1,2,nn)=-v(2,1,nn)/vmag
1325 v(2,2,nn)=v(1,1,nn)/vmag
1326 v(3,2,nn)=zero
1327 ELSE
1328 v(1,2,nn)=one
1329 v(2,2,nn)=zero
1330 v(3,2,nn)=zero
1331 ENDIF
1332 ENDDO
1333C . . . . . . . . . . . . .
1334C SOLUTION DOUBLE
1335C . . . . . . . . . . . . .
1336 DO n = 1, nindex2
1337 nn = index2(n)
1338 xmag(1,nn)=sqrt(a(1,1,nn)**2+a(2,1,nn)**2)
1339 xmag(2,nn)=sqrt(a(1,2,nn)**2+a(2,2,nn)**2)
1340 xmag(3,nn)=sqrt(a(1,3,nn)**2+a(2,3,nn)**2)
1341C
1342 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1343 ENDDO
1344C
1345#include "vectorize.inc"
1346 DO n = 1, nindex2
1347 nn = index2(n)
1348 IF(xmag(1,nn)==xmaxv(nn)) THEN
1349 lmaxv(nn) = 1
1350 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1351 lmaxv(nn) = 2
1352 ELSE
1353 lmaxv(nn) = 3
1354 ENDIF
1355C
1356 lmax = lmaxv(nn)
1357 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
1358 & <tol2(nn))THEN
1359 xmaxinv = one/xmaxv(nn)
1360 v(1,1,nn)= zero
1361 v(2,1,nn)= zero
1362 v(3,1,nn)= one
1363 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
1364 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
1365 v(3,2,nn)= zero
1366C
1367 ELSEIF(xmaxv(nn)>tol2(nn))THEN
1368 xmaxinv = one/xmaxv(nn)
1369 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
1370 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
1371 v(3,1,nn)= zero
1372 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
1373 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
1374 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
1375 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
1376 v(1,2,nn)=v(1,2,nn)*vmag
1377 v(2,2,nn)=v(2,2,nn)*vmag
1378 v(3,2,nn)=v(3,2,nn)*vmag
1379 ELSE
1380 v(1,1,nn) = one
1381 v(2,1,nn) = zero
1382 v(3,1,nn) = zero
1383 v(1,2,nn) = zero
1384 v(2,2,nn) = one
1385 v(3,2,nn) = zero
1386 ENDIF
1387 ENDDO
1388C
1389 DO nn = 1, nel
1390 vecdp(1,nn)=v(1,1,nn)
1391 vecdp(2,nn)=v(2,1,nn)
1392 vecdp(3,nn)=v(3,1,nn)
1393 vecdp(4,nn)=v(1,2,nn)
1394 vecdp(5,nn)=v(2,2,nn)
1395 vecdp(6,nn)=v(3,2,nn)
1396 vecdp(7,nn)=vecdp(2,nn)*vecdp(6,nn)-vecdp(3,nn)*vecdp(5,nn)
1397 vecdp(8,nn)=vecdp(3,nn)*vecdp(4,nn)-vecdp(1,nn)*vecdp(6,nn)
1398 vecdp(9,nn)=vecdp(1,nn)*vecdp(5,nn)-vecdp(2,nn)*vecdp(4,nn)
1399 ENDDO
1400C
1401 DO nn = 1, nel
1402 val(1,nn)=valdp(1,nn)
1403 val(2,nn)=valdp(2,nn)
1404 val(3,nn)=valdp(3,nn)
1405 vec(1,nn)=vecdp(1,nn)
1406 vec(2,nn)=vecdp(2,nn)
1407 vec(3,nn)=vecdp(3,nn)
1408 vec(4,nn)=vecdp(4,nn)
1409 vec(5,nn)=vecdp(5,nn)
1410 vec(6,nn)=vecdp(6,nn)
1411 vec(7,nn)=vecdp(7,nn)
1412 vec(8,nn)=vecdp(8,nn)
1413 vec(9,nn)=vecdp(9,nn)
1414 ENDDO
1415C Rewriting to bypass problem on itanium2 comp 9. + latency = 16
1416 DO n = 1, nindex3
1417 nn = index3(n)
1418C Str uses temporary table instead of temporary scalars
1419 str(1,nn)=vec(7,nn)
1420 str(2,nn)=vec(8,nn)
1421 str(3,nn)=vec(9,nn)
1422 ENDDO
1423 DO n = 1, nindex3
1424 nn = index3(n)
1425 vec(7,nn)=vec(1,nn)
1426 vec(8,nn)=vec(2,nn)
1427 vec(9,nn)=vec(3,nn)
1428 vec(1,nn)=-str(1,nn)
1429 vec(2,nn)=-str(2,nn)
1430 vec(3,nn)=-str(3,nn)
1431 ENDDO
1432 RETURN
1433 END
1434!||====================================================================
1435!|| epsf2u ../engine/source/materials/mat/mat033/sigeps33.F
1436!||--- called by ------------------------------------------------------
1437!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
1438!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
1439!||--- calls -----------------------------------------------------
1440!|| valpvec ../engine/source/materials/mat/mat033/sigeps33.F
1441!|| valpvecdp ../engine/source/materials/mat/mat033/sigeps33.F
1442!||====================================================================
1443 SUBROUTINE epsf2u(
1444 1 NEL ,FXX , FXY , FXZ , FYX ,
1445 2 FYY ,FYZ , FZX , FZY , FZZ ,
1446 3 UXX ,UYY , UZZ , UXY , UYZ ,
1447 4 UXZ )
1448C-----------------------------------------------
1449C I M P L I C I T T Y P E S
1450C-----------------------------------------------
1451#include "implicit_f.inc"
1452C-----------------------------------------------
1453C G L O B A L P A R A M E T E R S
1454C-----------------------------------------------
1455#include "mvsiz_p.inc"
1456C-----------------------------------------------
1457C C O M M O N
1458C-----------------------------------------------
1459#include "scr05_c.inc"
1460C----------------------------------------------------------------
1461C I N P U T A R G U M E N T S
1462C----------------------------------------------------------------
1463 INTEGER NEL
1464 my_real
1465 . FXX(NEL) , FXY(NEL), FXZ(NEL),
1466 . FYX(NEL) , FYY(NEL), FYZ(NEL),
1467 . FZX(NEL) , FZY(NEL), FZZ(NEL)
1468C----------------------------------------------------------------
1469C O U T P U T A R G U M E N T S
1470C----------------------------------------------------------------
1471 my_real
1472 . UXX(NEL) , UXY(NEL), UXZ(NEL),
1473 . UYY(NEL) , UYZ(NEL), UZZ(NEL)
1474C----------------------------------------------------------------
1475C I N P U T O U T P U T A R G U M E N T S
1476C----------------------------------------------------------------
1477C----------------------------------------------------------------
1478C L O C A L V A R I B L E S
1479C----------------------------------------------------------------
1480 INTEGER I
1481 my_real
1482 . EV(3,MVSIZ),AV(6,MVSIZ),EVV(3,MVSIZ),DIRPRV(3,3,MVSIZ)
1483C----------------------------------------------------------------
1484C--------- [C]=[F]^t[F] strain-----
1485 DO i=1,nel
1486 av(1,i)=fxx(i)*fxx(i)+fyx(i)*fyx(i)+fzx(i)*fzx(i)
1487 av(2,i)=fyy(i)*fyy(i)+fxy(i)*fxy(i)+fzy(i)*fzy(i)
1488 av(3,i)=fzz(i)*fzz(i)+fzx(i)*fzx(i)+fzy(i)*fzy(i)
1489 av(4,i)=fxx(i)*fxy(i)+fyx(i)*fyy(i)+fzx(i)*fzy(i)
1490 av(6,i)=fxx(i)*fxz(i)+fyx(i)*fyz(i)+fzx(i)*fzz(i)
1491 av(5,i)=fxz(i)*fxy(i)+fyz(i)*fyy(i)+fzz(i)*fzy(i)
1492 ENDDO
1493C for a simple precision executing
1494 IF (iresp==1) THEN
1495 CALL valpvecdp(av,evv,dirprv,nel)
1496 ELSE
1497 CALL valpvec(av,evv,dirprv,nel)
1498 ENDIF
1499 DO i=1,nel
1500 ev(1,i)=sqrt(evv(1,i))
1501 ev(2,i)=sqrt(evv(2,i))
1502 ev(3,i)=sqrt(evv(3,i))
1503 ENDDO
1504* TRANSFORM PRINCIPAL U TO GLOBAL DIRECTIONS
1505 DO i=1,nel
1506C
1507 uxx(i) = dirprv(1,1,i)*dirprv(1,1,i)*ev(1,i)
1508 . + dirprv(1,2,i)*dirprv(1,2,i)*ev(2,i)
1509 . + dirprv(1,3,i)*dirprv(1,3,i)*ev(3,i)
1510c
1511 uyy(i) = dirprv(2,2,i)*dirprv(2,2,i)*ev(2,i)
1512 . + dirprv(2,3,i)*dirprv(2,3,i)*ev(3,i)
1513 . + dirprv(2,1,i)*dirprv(2,1,i)*ev(1,i)
1514c
1515 uzz(i) = dirprv(3,3,i)*dirprv(3,3,i)*ev(3,i)
1516 . + dirprv(3,1,i)*dirprv(3,1,i)*ev(1,i)
1517 . + dirprv(3,2,i)*dirprv(3,2,i)*ev(2,i)
1518c
1519 uxy(i) = dirprv(1,1,i)*dirprv(2,1,i)*ev(1,i)
1520 . + dirprv(1,2,i)*dirprv(2,2,i)*ev(2,i)
1521 . + dirprv(1,3,i)*dirprv(2,3,i)*ev(3,i)
1522c
1523 uyz(i) = dirprv(2,2,i)*dirprv(3,2,i)*ev(2,i)
1524 . + dirprv(2,3,i)*dirprv(3,3,i)*ev(3,i)
1525 . + dirprv(2,1,i)*dirprv(3,1,i)*ev(1,i)
1526c
1527 uxz(i) = dirprv(3,3,i)*dirprv(1,3,i)*ev(3,i)
1528 . + dirprv(3,1,i)*dirprv(1,1,i)*ev(1,i)
1529 . + dirprv(3,2,i)*dirprv(1,2,i)*ev(2,i)
1530 ENDDO
1531C
1532 RETURN
1533 END
1534!||====================================================================
1535!|| valpvecop ../engine/source/materials/mat/mat033/sigeps33.F
1536!||--- calls -----------------------------------------------------
1537!|| floatmin ../common_source/tools/math/precision.c
1538!||====================================================================
1539 SUBROUTINE valpvecop(SIG,VAL,VEC,NEL)
1540C-----------------------------------------------
1541C I m p l i c i t T y p e s
1542C-----------------------------------------------
1543#include "implicit_f.inc"
1544C-----------------------------------------------
1545C G l o b a l P a r a m e t e r s
1546C-----------------------------------------------
1547#include "mvsiz_p.inc"
1548C-----------------------------------------------
1549C D u m m y A r g u m e n t s
1550C-----------------------------------------------
1551 my_real
1552 . sig(6,*), val(3,*), vec(9,*)
1553 INTEGER NEL
1554C-----------------------------------------------
1555C L o c a l V a r i a b l e s
1556C-----------------------------------------------
1557 INTEGER N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
1558 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
1559 DOUBLE PRECISION
1560 . CS(6), STR(3,MVSIZ), A(3,3,MVSIZ), V(3,3,MVSIZ), B(3,3,MVSIZ),
1561 . XMAG(3,MVSIZ), PR, AA, BB, AAA(MVSIZ),
1562 . CC, ANGP, DD, FTPI, TTPI, STRMAX,
1563 . TOL1(MVSIZ), TOL2(MVSIZ), XMAXV(MVSIZ), NORMINF(MVSIZ),
1564 . VMAG,
1565 .
1566 .
1567 . mdemi, xmaxinv, flm,
1568 . valdp(3,mvsiz),vecdp(9,mvsiz),s4_2,s5_2,s6_2,c1,c2,sq_aa
1569 REAL FLMIN
1570C-----------------------------------------------
1571C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
1572C FTPI=(4/3)*PI, TTPI=(2/3)*PI
1573C
1574C principal deviator of stress
1575C . . . . . . . . . . . . . . . . . . .
1576 mdemi = -half
1577 ttpi = acos(mdemi)
1578 ftpi = two*ttpi
1579C minimum precision depending on the type real or double
1580 CALL floatmin(cs(1),cs(2),flmin)
1581 flm = two*sqrt(flmin)
1582 nindex3=0
1583 c1 = half*sqrt(twenty7)
1584 c2 = two*sqrt(third)
1585 DO nn = 1, nel
1586 cs(1) = sig(1,nn)
1587 cs(2) = sig(2,nn)
1588 cs(3) = sig(3,nn)
1589 cs(4) = sig(4,nn)
1590 cs(5) = sig(5,nn)
1591 cs(6) = sig(6,nn)
1592 pr = -(cs(1)+cs(2)+cs(3)) * third
1593 cs(1) = cs(1) + pr
1594 cs(2) = cs(2) + pr
1595 cs(3) = cs(3) + pr
1596 s4_2 = cs(4)*cs(4)
1597 s5_2 = cs(5)*cs(5)
1598 s6_2 = cs(6)*cs(6)
1599 aaa(nn)=s4_2 + s5_2 + s6_2 - cs(1)*cs(2)
1600 & - cs(2)*cs(3) - cs(1)*cs(3)
1601 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
1602 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
1603 norminf(nn) = em10*norminf(nn)
1604C cas racine triple
1605 aa = max(aaa(nn),norminf(nn))
1606C
1607 bb=cs(1)*s5_2 + cs(2)*s6_2
1608 & + cs(3)*s4_2 - cs(1)*cs(2)*cs(3)
1609 & - two*cs(4)*cs(5)*cs(6)
1610C
1611 sq_aa = sqrt(aa)
1612 cc=-c1/max(em10,sq_aa)*bb/max(em20,aa)
1613c CC=-SQRT(TWENTY7/MAX(EM20,AA))*BB*HALF/MAX(EM20,AA)
1614 cc= min(cc,one)
1615 cc= max(cc,-one)
1616 angp=acos(cc) * third
1617 dd=c2*sq_aa
1618 str(1,nn)=dd*cos(angp)
1619 str(2,nn)=dd*cos(angp+ftpi)
1620 str(3,nn)=dd*cos(angp+ttpi)
1621C
1622 valdp(1,nn) = str(1,nn)-pr
1623 valdp(2,nn) = str(2,nn)-pr
1624 valdp(3,nn) = str(3,nn)-pr
1625C precision reinforcement in simple compression ----
1626 IF(abs(str(3,nn))>abs(str(1,nn))
1627 & .AND.aaa(nn)>norminf(nn)) THEN
1628 aa=str(1,nn)
1629 str(1,nn)=str(3,nn)
1630 str(3,nn)=aa
1631 nindex3 = nindex3+1
1632 index3(nindex3) = nn
1633 ENDIF
1634C . . . . . . . . . . .
1635C VECTEURS PROPRES
1636C . . . . . . . . . . .
1637 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
1638 tol1(nn)= max(em20,flm*strmax*strmax)
1639 tol2(nn)=flm*strmax * third
1640 a(1,1,nn)=cs(1)-str(1,nn)
1641 a(2,2,nn)=cs(2)-str(1,nn)
1642 a(3,3,nn)=cs(3)-str(1,nn)
1643 a(1,2,nn)=cs(4)
1644 a(2,1,nn)=cs(4)
1645 a(2,3,nn)=cs(5)
1646 a(3,2,nn)=cs(5)
1647 a(1,3,nn)=cs(6)
1648 a(3,1,nn)=cs(6)
1649C
1650 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
1651 . *a(2,2,nn)
1652 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
1653 . *a(2,3,nn)
1654 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
1655 . *a(2,1,nn)
1656 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
1657 . *a(3,2,nn)
1658 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
1659 . *a(3,3,nn)
1660 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
1661 . *a(3,1,nn)
1662 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
1663 . *a(1,2,nn)
1664 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
1665 . *a(1,3,nn)
1666 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
1667 . *a(1,1,nn)
1668 xmag(1,nn)=b(1,1,nn)*b(1,1,nn)+b(2,1,nn)*b(2,1,nn)+
1669 . b(3,1,nn)*b(3,1,nn)
1670 xmag(2,nn)=b(1,2,nn)*b(1,2,nn)+b(2,2,nn)*b(2,2,nn)+
1671 . b(3,2,nn)*b(3,2,nn)
1672 xmag(3,nn)=b(1,3,nn)*b(1,3,nn)+b(2,3,nn)*b(2,3,nn)+
1673 . b(3,3,nn)*b(3,3,nn)
1674
1675 ENDDO
1676C
1677 nindex1 = 0
1678 nindex2 = 0
1679 DO nn = 1, nel
1680 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1681 IF(xmag(1,nn)==xmaxv(nn)) THEN
1682 lmaxv(nn) = 1
1683 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1684 lmaxv(nn) = 2
1685 ELSE
1686 lmaxv(nn) = 3
1687 ENDIF
1688 IF(aaa(nn)<norminf(nn)) THEN
1689 valdp(1,nn) = sig(1,nn)
1690 valdp(2,nn) = sig(2,nn)
1691 valdp(3,nn) = sig(3,nn)
1692 v(1,1,nn) = one
1693 v(2,1,nn) = zero
1694 v(3,1,nn) = zero
1695 v(1,2,nn) = zero
1696 v(2,2,nn) = one
1697 v(3,2,nn) = zero
1698 ELSEIF(xmaxv(nn)>tol1(nn)*tol1(nn)) THEN
1699 nindex1 = nindex1 + 1
1700 index1(nindex1) = nn
1701 ELSE
1702 nindex2 = nindex2 + 1
1703 index2(nindex2) = nn
1704 ENDIF
1705 ENDDO
1706C
1707#include "vectorize.inc"
1708 DO n = 1, nindex1
1709 nn = index1(n)
1710 lmax = lmaxv(nn)
1711 xmaxinv = one/sqrt(xmaxv(nn))
1712 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
1713 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
1714 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
1715 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
1716 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
1717 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
1718C
1719 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
1720 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
1721 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
1722 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
1723 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
1724 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
1725 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
1726 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
1727 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
1728 xmag(1,nn)=b(1,1,nn)*b(1,1,nn)+b(2,1,nn)*b(2,1,nn)+
1729 . b(3,1,nn)*b(3,1,nn)
1730 xmag(2,nn)=b(1,2,nn)*b(1,2,nn)+b(2,2,nn)*b(2,2,nn)+
1731 . b(3,2,nn)*b(3,2,nn)
1732 xmag(3,nn)=b(1,3,nn)*b(1,3,nn)+b(2,3,nn)*b(2,3,nn)+
1733 . b(3,3,nn)*b(3,3,nn)
1734C
1735 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1736 ENDDO
1737C
1738#include "vectorize.inc"
1739 DO n = 1, nindex1
1740 nn = index1(n)
1741 IF(xmag(1,nn)==xmaxv(nn)) THEN
1742 lmaxv(nn) = 1
1743 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1744 lmaxv(nn) = 2
1745 ELSE
1746 lmaxv(nn) = 3
1747 ENDIF
1748C
1749 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
1750 lmax = lmaxv(nn)
1751 xmaxv(nn)=sqrt(xmaxv(nn))
1752 IF(xmaxv(nn)>tol2(nn))THEN
1753 xmaxinv = one/xmaxv(nn)
1754 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
1755 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
1756 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
1757 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
1758 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
1759 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
1760 c1 = v(1,2,nn)*v(1,2,nn)+v(2,2,nn)*v(2,2,nn)+
1761 & v(3,2,nn)*v(3,2,nn)
1762 vmag=one/sqrt(c1)
1763 v(1,2,nn)=v(1,2,nn)*vmag
1764 v(2,2,nn)=v(2,2,nn)*vmag
1765 v(3,2,nn)=v(3,2,nn)*vmag
1766 ELSEIF(vmag>tol2(nn))THEN
1767 v(1,2,nn)=-v(2,1,nn)/vmag
1768 v(2,2,nn)=v(1,1,nn)/vmag
1769 v(3,2,nn)=zero
1770 ELSE
1771 v(1,2,nn)=one
1772 v(2,2,nn)=zero
1773 v(3,2,nn)=zero
1774 ENDIF
1775 ENDDO
1776C . . . . . . . . . . . . .
1777C SOLUTION DOUBLE
1778C . . . . . . . . . . . . .
1779 DO n = 1, nindex2
1780 nn = index2(n)
1781 xmag(1,nn)=a(1,1,nn)*a(1,1,nn)+a(2,1,nn)*a(2,1,nn)
1782 xmag(2,nn)=a(1,2,nn)*a(1,2,nn)+a(2,2,nn)*a(2,2,nn)
1783 xmag(3,nn)=a(1,3,nn)*a(1,3,nn)+a(2,3,nn)*a(2,3,nn)
1784C
1785 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1786 ENDDO
1787C
1788#include "vectorize.inc"
1789 DO n = 1, nindex2
1790 nn = index2(n)
1791 IF(xmag(1,nn)==xmaxv(nn)) THEN
1792 lmaxv(nn) = 1
1793 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1794 lmaxv(nn) = 2
1795 ELSE
1796 lmaxv(nn) = 3
1797 ENDIF
1798C
1799 lmax = lmaxv(nn)
1800 xmaxv(nn)=sqrt(xmaxv(nn))
1801 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
1802 & <tol2(nn))THEN
1803 xmaxinv = one/xmaxv(nn)
1804 v(1,1,nn)= zero
1805 v(2,1,nn)= zero
1806 v(3,1,nn)= one
1807 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
1808 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
1809 v(3,2,nn)= zero
1810C
1811 ELSEIF(xmaxv(nn)>tol2(nn))THEN
1812 xmaxinv = one/xmaxv(nn)
1813 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
1814 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
1815 v(3,1,nn)= zero
1816 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
1817 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
1818 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
1819 c1 = v(1,2,nn)*v(1,2,nn)+v(2,2,nn)*v(2,2,nn)+
1820 & v(3,2,nn)*v(3,2,nn)
1821 vmag=one/sqrt(c1)
1822 v(1,2,nn)=v(1,2,nn)*vmag
1823 v(2,2,nn)=v(2,2,nn)*vmag
1824 v(3,2,nn)=v(3,2,nn)*vmag
1825 ELSE
1826 v(1,1,nn) = one
1827 v(2,1,nn) = zero
1828 v(3,1,nn) = zero
1829 v(1,2,nn) = zero
1830 v(2,2,nn) = one
1831 v(3,2,nn) = zero
1832 ENDIF
1833 ENDDO
1834C
1835 DO nn = 1, nel
1836 vecdp(1,nn)=v(1,1,nn)
1837 vecdp(2,nn)=v(2,1,nn)
1838 vecdp(3,nn)=v(3,1,nn)
1839 vecdp(4,nn)=v(1,2,nn)
1840 vecdp(5,nn)=v(2,2,nn)
1841 vecdp(6,nn)=v(3,2,nn)
1842 vecdp(7,nn)=vecdp(2,nn)*vecdp(6,nn)-vecdp(3,nn)*vecdp(5,nn)
1843 vecdp(8,nn)=vecdp(3,nn)*vecdp(4,nn)-vecdp(1,nn)*vecdp(6,nn)
1844 vecdp(9,nn)=vecdp(1,nn)*vecdp(5,nn)-vecdp(2,nn)*vecdp(4,nn)
1845 ENDDO
1846C
1847 DO nn = 1, nel
1848 val(1,nn)=valdp(1,nn)
1849 val(2,nn)=valdp(2,nn)
1850 val(3,nn)=valdp(3,nn)
1851 vec(1,nn)=vecdp(1,nn)
1852 vec(2,nn)=vecdp(2,nn)
1853 vec(3,nn)=vecdp(3,nn)
1854 vec(4,nn)=vecdp(4,nn)
1855 vec(5,nn)=vecdp(5,nn)
1856 vec(6,nn)=vecdp(6,nn)
1857 vec(7,nn)=vecdp(7,nn)
1858 vec(8,nn)=vecdp(8,nn)
1859 vec(9,nn)=vecdp(9,nn)
1860 ENDDO
1861C Rewriting to bypass problem on itanium2 comp 9. + latency = 16
1862 DO n = 1, nindex3
1863 nn = index3(n)
1864C Str uses temporary table instead of temporary scalars
1865 str(1,nn)=vec(7,nn)
1866 str(2,nn)=vec(8,nn)
1867 str(3,nn)=vec(9,nn)
1868 ENDDO
1869 DO n = 1, nindex3
1870 nn = index3(n)
1871 vec(7,nn)=vec(1,nn)
1872 vec(8,nn)=vec(2,nn)
1873 vec(9,nn)=vec(3,nn)
1874 vec(1,nn)=-str(1,nn)
1875 vec(2,nn)=-str(2,nn)
1876 vec(3,nn)=-str(3,nn)
1877 ENDDO
1878 RETURN
1879 END
1880
1881
1882!||====================================================================
1883!|| valpvec_v ../engine/source/materials/mat/mat033/sigeps33.F
1884!||--- called by ------------------------------------------------------
1885!|| sigeps124 ../engine/source/materials/mat/mat124/sigeps124.F
1886!|| sigeps130 ../engine/source/materials/mat/mat130/sigeps130.F90
1887!|| sigeps163 ../engine/source/materials/mat/mat163/sigeps163.F90
1888!|| sigeps190 ../engine/source/materials/mat/mat190/sigeps190.F
1889!|| sigeps33 ../engine/source/materials/mat/mat033/sigeps33.F
1890!|| sigeps38 ../engine/source/materials/mat/mat038/sigeps38.F
1891!|| sigeps42 ../engine/source/materials/mat/mat042/sigeps42.F
1892!|| sigeps62 ../engine/source/materials/mat/mat062/sigeps62.F
1893!|| sigeps69 ../engine/source/materials/mat/mat069/sigeps69.f
1894!|| sigeps71 ../engine/source/materials/mat/mat071/sigeps71.F
1895!|| sigeps82 ../engine/source/materials/mat/mat082/sigeps82.F
1896!|| sigeps88 ../engine/source/materials/mat/mat088/sigeps88.F90
1897!|| sigeps90 ../engine/source/materials/mat/mat090/sigeps90.F
1898!|| sigeps92 ../engine/source/materials/mat/mat092/sigeps92.F
1899!|| sigeps94 ../engine/source/materials/mat/mat094/sigeps94.F
1900!||--- calls -----------------------------------------------------
1901!|| floatmin ../common_source/tools/math/precision.c
1902!||====================================================================
1903 SUBROUTINE valpvec_v(SIG,VAL,VEC,NEL)
1904C-----------------------------------------------
1905C I m p l i c i t T y p e s
1906C-----------------------------------------------
1907#include "implicit_f.inc"
1908C-----------------------------------------------
1909C G l o b a l P a r a m e t e r s
1910C-----------------------------------------------
1911#include "mvsiz_p.inc"
1912C-----------------------------------------------
1913C D u m m y A r g u m e n t s
1914C-----------------------------------------------
1915 my_real
1916 . sig(mvsiz,6), val(mvsiz,3), vec(mvsiz,9)
1917 INTEGER NEL
1918C-----------------------------------------------
1919C L o c a l V a r i a b l e s
1920C-----------------------------------------------
1921 INTEGER N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
1922 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
1923 my_real
1924 . CS(6), STR(MVSIZ,3), A(MVSIZ,3,3), V(MVSIZ,3,3), B(MVSIZ,3,3),
1925 . XMAG(MVSIZ,3), PR, AA, BB, AAA(MVSIZ),
1926 . CC, ANGP, DD, FTPI, TTPI, STRMAX,
1927 . TOL1(MVSIZ), TOL2(MVSIZ), XMAXV(MVSIZ), NORMINF(MVSIZ),
1928 . VMAG,
1929 .
1930 .
1931 . mdemi, xmaxinv, flm
1932 REAL FLMIN
1933C-----------------------------------------------
1934C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
1935C FTPI=(4/3)*PI, TTPI=(2/3)*PI
1936C
1937C principal deviator of stress
1938C . . . . . . . . . . . . . . . . . . .
1939 mdemi = -half
1940 ttpi = acos(mdemi)
1941 ftpi = two*ttpi
1942C minimum precision depending on the type real or double
1943 CALL floatmin(cs(1),cs(2),flmin)
1944 flm = two*sqrt(flmin)
1945 nindex3=0
1946 DO nn = 1, nel
1947 cs(1) = sig(nn,1)
1948 cs(2) = sig(nn,2)
1949 cs(3) = sig(nn,3)
1950 cs(4) = sig(nn,4)
1951 cs(5) = sig(nn,5)
1952 cs(6) = sig(nn,6)
1953 pr = -(cs(1)+cs(2)+cs(3))* third
1954 cs(1) = cs(1) + pr
1955 cs(2) = cs(2) + pr
1956 cs(3) = cs(3) + pr
1957 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
1958 & - cs(2)*cs(3) - cs(1)*cs(3)
1959 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
1960 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
1961 norminf(nn) = em10*norminf(nn)
1962C cas racine triple
1963c AA = MAX(AAA(NN),NORMINF(NN),EM20)
1964 aa = max(aaa(nn),norminf(nn))
1965C
1966 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
1967 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
1968 & - two*cs(4)*cs(5)*cs(6)
1969C
1970 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
1971 cc= min(cc,one)
1972 cc= max(cc,-one)
1973 angp=acos(cc) * third
1974 dd=two*sqrt(aa * third)
1975 str(nn,1)=dd*cos(angp)
1976 str(nn,2)=dd*cos(angp+ftpi)
1977 str(nn,3)=dd*cos(angp+ttpi)
1978C
1979 val(nn,1) = str(nn,1)-pr
1980 val(nn,2) = str(nn,2)-pr
1981 val(nn,3) = str(nn,3)-pr
1982C accuracy reinforcement in compression ----
1983 IF(abs(str(nn,3))>abs(str(nn,1))
1984 & .AND.aaa(nn)>norminf(nn)) THEN
1985 aa=str(nn,1)
1986 str(nn,1)=str(nn,3)
1987 str(nn,3)=aa
1988 nindex3 = nindex3+1
1989 index3(nindex3) = nn
1990 ENDIF
1991C . . . . . . . . . . .
1992C VECTEURS PROPRES
1993C . . . . . . . . . . .
1994 strmax= max(abs(str(nn,1)),abs(str(nn,3)))
1995 tol1(nn)= max(em20,flm*strmax**2)
1996 tol2(nn)=flm*strmax/3
1997 a(nn,1,1)=cs(1)-str(nn,1)
1998 a(nn,2,2)=cs(2)-str(nn,1)
1999 a(nn,3,3)=cs(3)-str(nn,1)
2000 a(nn,1,2)=cs(4)
2001 a(nn,2,1)=cs(4)
2002 a(nn,2,3)=cs(5)
2003 a(nn,3,2)=cs(5)
2004 a(nn,1,3)=cs(6)
2005 a(nn,3,1)=cs(6)
2006C
2007 b(nn,1,1)=a(nn,2,1)*a(nn,3,2)-a(nn,3,1)
2008 . *a(nn,2,2)
2009 b(nn,1,2)=a(nn,2,2)*a(nn,3,3)-a(nn,3,2)
2010 . *a(nn,2,3)
2011 b(nn,1,3)=a(nn,2,3)*a(nn,3,1)-a(nn,3,3)
2012 . *a(nn,2,1)
2013 b(nn,2,1)=a(nn,3,1)*a(nn,1,2)-a(nn,1,1)
2014 . *a(nn,3,2)
2015 b(nn,2,2)=a(nn,3,2)*a(nn,1,3)-a(nn,1,2)
2016 . *a(nn,3,3)
2017 b(nn,2,3)=a(nn,3,3)*a(nn,1,1)-a(nn,1,3)
2018 . *a(nn,3,1)
2019 b(nn,3,1)=a(nn,1,1)*a(nn,2,2)-a(nn,2,1)
2020 . *a(nn,1,2)
2021 b(nn,3,2)=a(nn,1,2)*a(nn,2,3)-a(nn,2,2)
2022 . *a(nn,1,3)
2023 b(nn,3,3)=a(nn,1,3)*a(nn,2,1)-a(nn,2,3)
2024 . *a(nn,1,1)
2025 xmag(nn,1)=sqrt(b(nn,1,1)**2+b(nn,2,1)**2+b(nn,3,1)**2)
2026 xmag(nn,2)=sqrt(b(nn,1,2)**2+b(nn,2,2)**2+b(nn,3,2)**2)
2027 xmag(nn,3)=sqrt(b(nn,1,3)**2+b(nn,2,3)**2+b(nn,3,3)**2)
2028
2029 ENDDO
2030C
2031 nindex1 = 0
2032 nindex2 = 0
2033 DO nn = 1, nel
2034 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2035 IF(xmag(nn,1)==xmaxv(nn)) THEN
2036 lmaxv(nn) = 1
2037 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2038 lmaxv(nn) = 2
2039 ELSE
2040 lmaxv(nn) = 3
2041 ENDIF
2042 IF(aaa(nn)<norminf(nn)) THEN
2043 val(nn,1) = sig(nn,1)
2044 val(nn,2) = sig(nn,2)
2045 val(nn,3) = sig(nn,3)
2046 v(nn,1,1) = one
2047 v(nn,2,1) = zero
2048 v(nn,3,1) = zero
2049 v(nn,1,2) = zero
2050 v(nn,2,2) = one
2051 v(nn,3,2) = zero
2052
2053 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
2054 nindex1 = nindex1 + 1
2055 index1(nindex1) = nn
2056 ELSE
2057 nindex2 = nindex2 + 1
2058 index2(nindex2) = nn
2059 ENDIF
2060 ENDDO
2061C
2062#include "vectorize.inc"
2063 DO n = 1, nindex1
2064 nn = index1(n)
2065 lmax = lmaxv(nn)
2066 xmaxinv = one/xmaxv(nn)
2067 v(nn,1,1)=b(nn,1,lmax)*xmaxinv
2068 v(nn,2,1)=b(nn,2,lmax)*xmaxinv
2069 v(nn,3,1)=b(nn,3,lmax)*xmaxinv
2070 a(nn,1,1)=a(nn,1,1)+str(nn,1)-str(nn,3)
2071 a(nn,2,2)=a(nn,2,2)+str(nn,1)-str(nn,3)
2072 a(nn,3,3)=a(nn,3,3)+str(nn,1)-str(nn,3)
2073C
2074 b(nn,1,1)=a(nn,2,1)*v(nn,3,1)-a(nn,3,1)*v(nn,2,1)
2075 b(nn,1,2)=a(nn,2,2)*v(nn,3,1)-a(nn,3,2)*v(nn,2,1)
2076 b(nn,1,3)=a(nn,2,3)*v(nn,3,1)-a(nn,3,3)*v(nn,2,1)
2077 b(nn,2,1)=a(nn,3,1)*v(nn,1,1)-a(nn,1,1)*v(nn,3,1)
2078 b(nn,2,2)=a(nn,3,2)*v(nn,1,1)-a(nn,1,2)*v(nn,3,1)
2079 b(nn,2,3)=a(nn,3,3)*v(nn,1,1)-a(nn,1,3)*v(nn,3,1)
2080 b(nn,3,1)=a(nn,1,1)*v(nn,2,1)-a(nn,2,1)*v(nn,1,1)
2081 b(nn,3,2)=a(nn,1,2)*v(nn,2,1)-a(nn,2,2)*v(nn,1,1)
2082 b(nn,3,3)=a(nn,1,3)*v(nn,2,1)-a(nn,2,3)*v(nn,1,1)
2083 xmag(nn,1)=sqrt(b(nn,1,1)**2+b(nn,2,1)**2+b(nn,3,1)**2)
2084 xmag(nn,2)=sqrt(b(nn,1,2)**2+b(nn,2,2)**2+b(nn,3,2)**2)
2085 xmag(nn,3)=sqrt(b(nn,1,3)**2+b(nn,2,3)**2+b(nn,3,3)**2)
2086C
2087 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2088 ENDDO
2089C
2090#include "vectorize.inc"
2091 DO n = 1, nindex1
2092 nn = index1(n)
2093 IF(xmag(nn,1)==xmaxv(nn)) THEN
2094 lmaxv(nn) = 1
2095 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2096 lmaxv(nn) = 2
2097 ELSE
2098 lmaxv(nn) = 3
2099 ENDIF
2100C
2101 vmag=sqrt(v(nn,1,1)**2+v(nn,2,1)**2)
2102 lmax = lmaxv(nn)
2103 IF(xmaxv(nn)>tol2(nn))THEN
2104 xmaxinv = one/xmaxv(nn)
2105 v(nn,1,3)=b(nn,1,lmax)*xmaxinv
2106 v(nn,2,3)=b(nn,2,lmax)*xmaxinv
2107 v(nn,3,3)=b(nn,3,lmax)*xmaxinv
2108 v(nn,1,2)=v(nn,2,3)*v(nn,3,1)-v(nn,2,1)*v(nn,3,3)
2109 v(nn,2,2)=v(nn,3,3)*v(nn,1,1)-v(nn,3,1)*v(nn,1,3)
2110 v(nn,3,2)=v(nn,1,3)*v(nn,2,1)-v(nn,1,1)*v(nn,2,3)
2111 vmag=one/sqrt(v(nn,1,2)**2+v(nn,2,2)**2+v(nn,3,2)**2)
2112 v(nn,1,2)=v(nn,1,2)*vmag
2113 v(nn,2,2)=v(nn,2,2)*vmag
2114 v(nn,3,2)=v(nn,3,2)*vmag
2115 ELSEIF(vmag>tol2(nn))THEN
2116 v(nn,1,2)=-v(nn,2,1)/vmag
2117 v(nn,2,2)=v(nn,1,1)/vmag
2118 v(nn,3,2)=zero
2119 ELSE
2120 v(nn,1,2)=one
2121 v(nn,2,2)=zero
2122 v(nn,3,2)=zero
2123 ENDIF
2124 ENDDO
2125C . . . . . . . . . . . . .
2126C SOLUTION DOUBLE
2127C . . . . . . . . . . . . .
2128 DO n = 1, nindex2
2129 nn = index2(n)
2130 xmag(nn,1)=sqrt(a(nn,1,1)**2+a(nn,2,1)**2)
2131 xmag(nn,2)=sqrt(a(nn,1,2)**2+a(nn,2,2)**2)
2132 xmag(nn,3)=sqrt(a(nn,1,3)**2+a(nn,2,3)**2)
2133C
2134 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2135 ENDDO
2136C
2137#include "vectorize.inc"
2138 DO n = 1, nindex2
2139 nn = index2(n)
2140 IF(xmag(nn,1)==xmaxv(nn)) THEN
2141 lmaxv(nn) = 1
2142 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2143 lmaxv(nn) = 2
2144 ELSE
2145 lmaxv(nn) = 3
2146 ENDIF
2147C
2148 lmax = lmaxv(nn)
2149 IF(max(abs(a(nn,3,1)),abs(a(nn,3,2)),abs(a(nn,3,3)))
2150 & <tol2(nn))THEN
2151 xmaxinv = one/xmaxv(nn)
2152 v(nn,1,1)= zero
2153 v(nn,2,1)= zero
2154 v(nn,3,1)= one
2155 v(nn,1,2)=-a(nn,2,lmax)*xmaxinv
2156 v(nn,2,2)= a(nn,1,lmax)*xmaxinv
2157 v(nn,3,2)= zero
2158C
2159 ELSEIF(xmaxv(nn)>tol2(nn))THEN
2160 xmaxinv = one/xmaxv(nn)
2161 v(nn,1,1)=-a(nn,2,lmax)*xmaxinv
2162 v(nn,2,1)= a(nn,1,lmax)*xmaxinv
2163 v(nn,3,1)= zero
2164 v(nn,1,2)=-a(nn,3,lmax)*v(nn,2,1)
2165 v(nn,2,2)= a(nn,3,lmax)*v(nn,1,1)
2166 v(nn,3,2)= a(nn,1,lmax)*v(nn,2,1)-a(nn,2,lmax)*v(nn,1,1)
2167 vmag=one/sqrt(v(nn,1,2)**2+v(nn,2,2)**2+v(nn,3,2)**2)
2168 v(nn,1,2)=v(nn,1,2)*vmag
2169 v(nn,2,2)=v(nn,2,2)*vmag
2170 v(nn,3,2)=v(nn,3,2)*vmag
2171 ELSE
2172 v(nn,1,1) = one
2173 v(nn,2,1) = zero
2174 v(nn,3,1) = zero
2175 v(nn,1,2) = zero
2176 v(nn,2,2) = one
2177 v(nn,3,2) = zero
2178 ENDIF
2179 ENDDO
2180C
2181 DO nn = 1,nel
2182 vec(nn,1)=v(nn,1,1)
2183 vec(nn,2)=v(nn,2,1)
2184 vec(nn,3)=v(nn,3,1)
2185 vec(nn,4)=v(nn,1,2)
2186 vec(nn,5)=v(nn,2,2)
2187 vec(nn,6)=v(nn,3,2)
2188 vec(nn,7)=vec(nn,2)*vec(nn,6)-vec(nn,3)*vec(nn,5)
2189 vec(nn,8)=vec(nn,3)*vec(nn,4)-vec(nn,1)*vec(nn,6)
2190 vec(nn,9)=vec(nn,1)*vec(nn,5)-vec(nn,2)*vec(nn,4)
2191 ENDDO
2192C Rewriting to bypass problem on itanium2 comp 9. + latency = 16
2193 DO n = 1, nindex3
2194 nn = index3(n)
2195C Str uses temporary table instead of temporary scalars
2196 str(nn,1)=vec(nn,7)
2197 str(nn,2)=vec(nn,8)
2198 str(nn,3)=vec(nn,9)
2199 ENDDO
2200 DO n = 1, nindex3
2201 nn = index3(n)
2202 vec(nn,7)=vec(nn,1)
2203 vec(nn,8)=vec(nn,2)
2204 vec(nn,9)=vec(nn,3)
2205 vec(nn,1)=-str(nn,1)
2206 vec(nn,2)=-str(nn,2)
2207 vec(nn,3)=-str(nn,3)
2208 ENDDO
2209C
2210 RETURN
2211 END
2212!||====================================================================
2213!|| valpvecdp_v ../engine/source/materials/mat/mat033/sigeps33.F
2214!||--- called by ------------------------------------------------------
2215!|| fail_spalling_s ../engine/source/materials/fail/spalling/fail_spalling_s.F90
2216!|| sigeps124 ../engine/source/materials/mat/mat124/sigeps124.F
2217!|| sigeps130 ../engine/source/materials/mat/mat130/sigeps130.F90
2218!|| sigeps163 ../engine/source/materials/mat/mat163/sigeps163.F90
2219!|| sigeps190 ../engine/source/materials/mat/mat190/sigeps190.F
2220!|| sigeps33 ../engine/source/materials/mat/mat033/sigeps33.F
2221!|| sigeps38 ../engine/source/materials/mat/mat038/sigeps38.F
2222!|| sigeps42 ../engine/source/materials/mat/mat042/sigeps42.F
2223!|| sigeps62 ../engine/source/materials/mat/mat062/sigeps62.F
2224!|| sigeps69 ../engine/source/materials/mat/mat069/sigeps69.F
2225!|| sigeps71 ../engine/source/materials/mat/mat071/sigeps71.F
2226!|| sigeps82 ../engine/source/materials/mat/mat082/sigeps82.F
2227!|| sigeps88 ../engine/source/materials/mat/mat088/sigeps88.F90
2228!|| sigeps90 ../engine/source/materials/mat/mat090/sigeps90.F
2229!|| sigeps92 ../engine/source/materials/mat/mat092/sigeps92.F
2230!|| sigeps94 ../engine/source/materials/mat/mat094/sigeps94.F
2231!||--- calls -----------------------------------------------------
2232!|| floatmin ../common_source/tools/math/precision.c
2233!||====================================================================
2234 SUBROUTINE valpvecdp_v(SIG,VAL,VEC,NEL)
2235C-----------------------------------------------
2236C I m p l i c i t T y p e s
2237C-----------------------------------------------
2238#include "implicit_f.inc"
2239C-----------------------------------------------
2240C G l o b a l P a r a m e t e r s
2241C-----------------------------------------------
2242#include "mvsiz_p.inc"
2243C-----------------------------------------------
2244C D u m m y A r g u m e n t s
2245C-----------------------------------------------
2246 my_real
2247 . sig(mvsiz,6), val(mvsiz,3), vec(mvsiz,9)
2248 INTEGER NEL
2249C-----------------------------------------------
2250C L o c a l V a r i a b l e s
2251C-----------------------------------------------
2252 INTEGER N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
2253 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
2254 DOUBLE PRECISION
2255 . CS(6), STR(MVSIZ,3), A(MVSIZ,3,3), V(MVSIZ,3,3), B(MVSIZ,3,3),
2256 . XMAG(MVSIZ,3), PR, AA, BB, AAA(MVSIZ),
2257 . CC, ANGP, DD, FTPI, TTPI, STRMAX,
2258 . TOL1(MVSIZ), TOL2(MVSIZ), XMAXV(MVSIZ), NORMINF(MVSIZ),
2259 . VMAG,
2260 .
2261 .
2262 . mdemi, xmaxinv, flm,
2263 . valdp(mvsiz,3),vecdp(mvsiz,9)
2264 REAL FLMIN
2265C-----------------------------------------------
2266C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
2267C FTPI=(4/3)*PI, TTPI=(2/3)*PI
2268C
2269C principal deviator of stress
2270C . . . . . . . . . . . . . . . . . . .
2271 mdemi = -half
2272 ttpi = acos(mdemi)
2273 ftpi = two*ttpi
2274C minimum precision depending on the type real or double
2275 CALL floatmin(cs(1),cs(2),flmin)
2276 flm = two*sqrt(flmin)
2277 nindex3=0
2278 DO nn = 1, nel
2279 cs(1) = sig(nn,1)
2280 cs(2) = sig(nn,2)
2281 cs(3) = sig(nn,3)
2282 cs(4) = sig(nn,4)
2283 cs(5) = sig(nn,5)
2284 cs(6) = sig(nn,6)
2285 pr = -(cs(1)+cs(2)+cs(3)) * third
2286 cs(1) = cs(1) + pr
2287 cs(2) = cs(2) + pr
2288 cs(3) = cs(3) + pr
2289 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
2290 & - cs(2)*cs(3) - cs(1)*cs(3)
2291 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
2292 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
2293 norminf(nn) = em10*norminf(nn)
2294C cas racine triple
2295c AA = MAX(AAA(NN),NORMINF(NN),EM20)
2296 aa = max(aaa(nn),norminf(nn))
2297C
2298 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
2299 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
2300 & - two*cs(4)*cs(5)*cs(6)
2301C
2302 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
2303 cc= min(cc,one)
2304 cc= max(cc,-one)
2305 angp=acos(cc) * third
2306 dd=two*sqrt(aa * third)
2307 str(nn,1)=dd*cos(angp)
2308 str(nn,2)=dd*cos(angp+ftpi)
2309 str(nn,3)=dd*cos(angp+ttpi)
2310C
2311 valdp(nn,1) = str(nn,1)-pr
2312 valdp(nn,2) = str(nn,2)-pr
2313 valdp(nn,3) = str(nn,3)-pr
2314C precision reinforcement in simple compression ----
2315 IF(abs(str(nn,3))>abs(str(nn,1))
2316 & .AND.aaa(nn)>norminf(nn)) THEN
2317 aa=str(nn,1)
2318 str(nn,1)=str(nn,3)
2319 str(nn,3)=aa
2320 nindex3 = nindex3+1
2321 index3(nindex3) = nn
2322 ENDIF
2323C . . . . . . . . . . .
2324C VECTEURS PROPRES
2325C . . . . . . . . . . .
2326 strmax= max(abs(str(nn,1)),abs(str(nn,3)))
2327 tol1(nn)= max(em20,flm*strmax**2)
2328 tol2(nn)=flm*strmax * third
2329 a(nn,1,1)=cs(1)-str(nn,1)
2330 a(nn,2,2)=cs(2)-str(nn,1)
2331 a(nn,3,3)=cs(3)-str(nn,1)
2332 a(nn,1,2)=cs(4)
2333 a(nn,2,1)=cs(4)
2334 a(nn,2,3)=cs(5)
2335 a(nn,3,2)=cs(5)
2336 a(nn,1,3)=cs(6)
2337 a(nn,3,1)=cs(6)
2338C
2339 b(nn,1,1)=a(nn,2,1)*a(nn,3,2)-a(nn,3,1)
2340 . *a(nn,2,2)
2341 b(nn,1,2)=a(nn,2,2)*a(nn,3,3)-a(nn,3,2)
2342 . *a(nn,2,3)
2343 b(nn,1,3)=a(nn,2,3)*a(nn,3,1)-a(nn,3,3)
2344 . *a(nn,2,1)
2345 b(nn,2,1)=a(nn,3,1)*a(nn,1,2)-a(nn,1,1)
2346 . *a(nn,3,2)
2347 b(nn,2,2)=a(nn,3,2)*a(nn,1,3)-a(nn,1,2)
2348 . *a(nn,3,3)
2349 b(nn,2,3)=a(nn,3,3)*a(nn,1,1)-a(nn,1,3)
2350 . *a(nn,3,1)
2351 b(nn,3,1)=a(nn,1,1)*a(nn,2,2)-a(nn,2,1)
2352 . *a(nn,1,2)
2353 b(nn,3,2)=a(nn,1,2)*a(nn,2,3)-a(nn,2,2)
2354 . *a(nn,1,3)
2355 b(nn,3,3)=a(nn,1,3)*a(nn,2,1)-a(nn,2,3)
2356 . *a(nn,1,1)
2357 xmag(nn,1)=sqrt(b(nn,1,1)**2+b(nn,2,1)**2+b(nn,3,1)**2)
2358 xmag(nn,2)=sqrt(b(nn,1,2)**2+b(nn,2,2)**2+b(nn,3,2)**2)
2359 xmag(nn,3)=sqrt(b(nn,1,3)**2+b(nn,2,3)**2+b(nn,3,3)**2)
2360
2361 ENDDO
2362C
2363 nindex1 = 0
2364 nindex2 = 0
2365 DO nn = 1, nel
2366 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2367 IF(xmag(nn,1)==xmaxv(nn)) THEN
2368 lmaxv(nn) = 1
2369 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2370 lmaxv(nn) = 2
2371 ELSE
2372 lmaxv(nn) = 3
2373 ENDIF
2374 IF(aaa(nn)<norminf(nn)) THEN
2375 valdp(nn,1) = sig(nn,1)
2376 valdp(nn,2) = sig(nn,2)
2377 valdp(nn,3) = sig(nn,3)
2378 v(nn,1,1) = one
2379 v(nn,2,1) = zero
2380 v(nn,3,1) = zero
2381 v(nn,1,2) = zero
2382 v(nn,2,2) = one
2383 v(nn,3,2) = zero
2384 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
2385 nindex1 = nindex1 + 1
2386 index1(nindex1) = nn
2387 ELSE
2388 nindex2 = nindex2 + 1
2389 index2(nindex2) = nn
2390 ENDIF
2391 ENDDO
2392C
2393#include "vectorize.inc"
2394 DO n = 1, nindex1
2395 nn = index1(n)
2396 lmax = lmaxv(nn)
2397 xmaxinv = one/xmaxv(nn)
2398 v(nn,1,1)=b(nn,1,lmax)*xmaxinv
2399 v(nn,2,1)=b(nn,2,lmax)*xmaxinv
2400 v(nn,3,1)=b(nn,3,lmax)*xmaxinv
2401 a(nn,1,1)=a(nn,1,1)+str(nn,1)-str(nn,3)
2402 a(nn,2,2)=a(nn,2,2)+str(nn,1)-str(nn,3)
2403 a(nn,3,3)=a(nn,3,3)+str(nn,1)-str(nn,3)
2404C
2405 b(nn,1,1)=a(nn,2,1)*v(nn,3,1)-a(nn,3,1)*v(nn,2,1)
2406 b(nn,1,2)=a(nn,2,2)*v(nn,3,1)-a(nn,3,2)*v(nn,2,1)
2407 b(nn,1,3)=a(nn,2,3)*v(nn,3,1)-a(nn,3,3)*v(nn,2,1)
2408 b(nn,2,1)=a(nn,3,1)*v(nn,1,1)-a(nn,1,1)*v(nn,3,1)
2409 b(nn,2,2)=a(nn,3,2)*v(nn,1,1)-a(nn,1,2)*v(nn,3,1)
2410 b(nn,2,3)=a(nn,3,3)*v(nn,1,1)-a(nn,1,3)*v(nn,3,1)
2411 b(nn,3,1)=a(nn,1,1)*v(nn,2,1)-a(nn,2,1)*v(nn,1,1)
2412 b(nn,3,2)=a(nn,1,2)*v(nn,2,1)-a(nn,2,2)*v(nn,1,1)
2413 b(nn,3,3)=a(nn,1,3)*v(nn,2,1)-a(nn,2,3)*v(nn,1,1)
2414 xmag(nn,1)=sqrt(b(nn,1,1)**2+b(nn,2,1)**2+b(nn,3,1)**2)
2415 xmag(nn,2)=sqrt(b(nn,1,2)**2+b(nn,2,2)**2+b(nn,3,2)**2)
2416 xmag(nn,3)=sqrt(b(nn,1,3)**2+b(nn,2,3)**2+b(nn,3,3)**2)
2417C
2418 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2419 ENDDO
2420C
2421#include "vectorize.inc"
2422 DO n = 1, nindex1
2423 nn = index1(n)
2424 IF(xmag(nn,1)==xmaxv(nn)) THEN
2425 lmaxv(nn) = 1
2426 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2427 lmaxv(nn) = 2
2428 ELSE
2429 lmaxv(nn) = 3
2430 ENDIF
2431C
2432 vmag=sqrt(v(nn,1,1)**2+v(nn,2,1)**2)
2433 lmax = lmaxv(nn)
2434 IF(xmaxv(nn)>tol2(nn))THEN
2435 xmaxinv = one/xmaxv(nn)
2436 v(nn,1,3)=b(nn,1,lmax)*xmaxinv
2437 v(nn,2,3)=b(nn,2,lmax)*xmaxinv
2438 v(nn,3,3)=b(nn,3,lmax)*xmaxinv
2439 v(nn,1,2)=v(nn,2,3)*v(nn,3,1)-v(nn,2,1)*v(nn,3,3)
2440 v(nn,2,2)=v(nn,3,3)*v(nn,1,1)-v(nn,3,1)*v(nn,1,3)
2441 v(nn,3,2)=v(nn,1,3)*v(nn,2,1)-v(nn,1,1)*v(nn,2,3)
2442 vmag=one/sqrt(v(nn,1,2)**2+v(nn,2,2)**2+v(nn,3,2)**2)
2443 v(nn,1,2)=v(nn,1,2)*vmag
2444 v(nn,2,2)=v(nn,2,2)*vmag
2445 v(nn,3,2)=v(nn,3,2)*vmag
2446 ELSEIF(vmag>tol2(nn))THEN
2447 v(nn,1,2)=-v(nn,2,1)/vmag
2448 v(nn,2,2)=v(nn,1,1)/vmag
2449 v(nn,3,2)=zero
2450 ELSE
2451 v(nn,1,2)=one
2452 v(nn,2,2)=zero
2453 v(nn,3,2)=zero
2454 ENDIF
2455 ENDDO
2456C . . . . . . . . . . . . .
2457C SOLUTION DOUBLE
2458C . . . . . . . . . . . . .
2459#include "vectorize.inc"
2460 DO n = 1, nindex2
2461 nn = index2(n)
2462 xmag(nn,1)=sqrt(a(nn,1,1)**2+a(nn,2,1)**2)
2463 xmag(nn,2)=sqrt(a(nn,1,2)**2+a(nn,2,2)**2)
2464 xmag(nn,3)=sqrt(a(nn,1,3)**2+a(nn,2,3)**2)
2465C
2466 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2467 ENDDO
2468C
2469#include "vectorize.inc"
2470 DO n = 1, nindex2
2471 nn = index2(n)
2472 IF(xmag(nn,1)==xmaxv(nn)) THEN
2473 lmaxv(nn) = 1
2474 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2475 lmaxv(nn) = 2
2476 ELSE
2477 lmaxv(nn) = 3
2478 ENDIF
2479C
2480 lmax = lmaxv(nn)
2481 IF(max(abs(a(nn,3,1)),abs(a(nn,3,2)),abs(a(nn,3,3)))
2482 & <tol2(nn))THEN
2483 xmaxinv = one/xmaxv(nn)
2484 v(nn,1,1)= zero
2485 v(nn,2,1)= zero
2486 v(nn,3,1)= one
2487 v(nn,1,2)=-a(nn,2,lmax)*xmaxinv
2488 v(nn,2,2)= a(nn,1,lmax)*xmaxinv
2489 v(nn,3,2)= zero
2490C
2491 ELSEIF(xmaxv(nn)>tol2(nn))THEN
2492 xmaxinv = one/xmaxv(nn)
2493 v(nn,1,1)=-a(nn,2,lmax)*xmaxinv
2494 v(nn,2,1)= a(nn,1,lmax)*xmaxinv
2495 v(nn,3,1)= zero
2496 v(nn,1,2)=-a(nn,3,lmax)*v(nn,2,1)
2497 v(nn,2,2)= a(nn,3,lmax)*v(nn,1,1)
2498 v(nn,3,2)= a(nn,1,lmax)*v(nn,2,1)-a(nn,2,lmax)*v(nn,1,1)
2499 vmag=one/sqrt(v(nn,1,2)**2+v(nn,2,2)**2+v(nn,3,2)**2)
2500 v(nn,1,2)=v(nn,1,2)*vmag
2501 v(nn,2,2)=v(nn,2,2)*vmag
2502 v(nn,3,2)=v(nn,3,2)*vmag
2503 ELSE
2504 v(nn,1,1) = one
2505 v(nn,2,1) = zero
2506 v(nn,3,1) = zero
2507 v(nn,1,2) = zero
2508 v(nn,2,2) = one
2509 v(nn,3,2) = zero
2510 ENDIF
2511 ENDDO
2512C
2513 DO nn = 1, nel
2514 vecdp(nn,1)=v(nn,1,1)
2515 vecdp(nn,2)=v(nn,2,1)
2516 vecdp(nn,3)=v(nn,3,1)
2517 vecdp(nn,4)=v(nn,1,2)
2518 vecdp(nn,5)=v(nn,2,2)
2519 vecdp(nn,6)=v(nn,3,2)
2520 vecdp(nn,7)=vecdp(nn,2)*vecdp(nn,6)-vecdp(nn,3)*vecdp(nn,5)
2521 vecdp(nn,8)=vecdp(nn,3)*vecdp(nn,4)-vecdp(nn,1)*vecdp(nn,6)
2522 vecdp(nn,9)=vecdp(nn,1)*vecdp(nn,5)-vecdp(nn,2)*vecdp(nn,4)
2523 ENDDO
2524C
2525 DO nn = 1, nel
2526 val(nn,1)=valdp(nn,1)
2527 val(nn,2)=valdp(nn,2)
2528 val(nn,3)=valdp(nn,3)
2529 vec(nn,1)=vecdp(nn,1)
2530 vec(nn,2)=vecdp(nn,2)
2531 vec(nn,3)=vecdp(nn,3)
2532 vec(nn,4)=vecdp(nn,4)
2533 vec(nn,5)=vecdp(nn,5)
2534 vec(nn,6)=vecdp(nn,6)
2535 vec(nn,7)=vecdp(nn,7)
2536 vec(nn,8)=vecdp(nn,8)
2537 vec(nn,9)=vecdp(nn,9)
2538 ENDDO
2539C Rewriting to bypass problem on itanium2 comp 9. + latency = 16
2540#include "vectorize.inc"
2541 DO n = 1, nindex3
2542 nn = index3(n)
2543C Str uses temporary table instead of temporary scalars
2544 str(nn,1)=vec(nn,7)
2545 str(nn,2)=vec(nn,8)
2546 str(nn,3)=vec(nn,9)
2547 ENDDO
2548#include "vectorize.inc"
2549 DO n = 1, nindex3
2550 nn = index3(n)
2551 vec(nn,7)=vec(nn,1)
2552 vec(nn,8)=vec(nn,2)
2553 vec(nn,9)=vec(nn,3)
2554 vec(nn,1)=-str(nn,1)
2555 vec(nn,2)=-str(nn,2)
2556 vec(nn,3)=-str(nn,3)
2557 ENDDO
2558 RETURN
2559 END
2560!||====================================================================
2561!|| valpvecop_v ../engine/source/materials/mat/mat033/sigeps33.F
2562!||--- called by ------------------------------------------------------
2563!|| princ_u ../engine/source/materials/mat/mat001/princ_u.F
2564!|| princ_u1 ../engine/source/materials/mat/mat001/princ_u1.F
2565!||--- calls -----------------------------------------------------
2566!|| floatmin ../common_source/tools/math/precision.c
2567!||====================================================================
2568 SUBROUTINE valpvecop_v(SIG,VAL,VEC,NEL)
2569C-----------------------------------------------
2570C I m p l i c i t T y p e s
2571C-----------------------------------------------
2572#include "implicit_f.inc"
2573C-----------------------------------------------
2574C G l o b a l P a r a m e t e r s
2575C-----------------------------------------------
2576#include "mvsiz_p.inc"
2577C-----------------------------------------------
2578C D u m m y A r g u m e n t s
2579C-----------------------------------------------
2580 my_real
2581 . sig(mvsiz,6), val(mvsiz,3), vec(mvsiz,9)
2582 INTEGER NEL
2583C-----------------------------------------------
2584C L o c a l V a r i a b l e s
2585C-----------------------------------------------
2586 INTEGER N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
2587 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
2588 DOUBLE PRECISION
2589 . CS1,CS2,CS3,CS4,CS5,CS6,
2590 . A(MVSIZ,3,3), V(MVSIZ,3,3), B(MVSIZ,3,3),
2591 . PR, AA, BB, AAA(MVSIZ),
2592 . CC, ANGP, DD, FTPI, TTPI, STRMAX,
2593 . TOL1(MVSIZ), TOL2(MVSIZ), XMAXV(MVSIZ), NORMINF(MVSIZ),
2594 . vmag,
2595 .
2596 . s1,s2,s3,
2597 .
2598 . mdemi, xmaxinv, flm,
2599 . s4_2,s5_2,s6_2,c1,c2,sq_aa,
2600 . valdp1(mvsiz),valdp2(mvsiz),valdp3(mvsiz),
2601 . str1(mvsiz),str2(mvsiz),str3(mvsiz),
2602 . xmag1(mvsiz),xmag2(mvsiz),xmag3(mvsiz),
2603 . vecdp1(mvsiz),vecdp2(mvsiz),vecdp3(mvsiz),
2604 . vecdp4(mvsiz),vecdp5(mvsiz),vecdp6(mvsiz),
2605 . vecdp7(mvsiz),vecdp8(mvsiz),vecdp9(mvsiz)
2606
2607 DOUBLE PRECISION :: SQRT3DMI, SIN_ANGP, COS_ANGP
2608 PARAMETER (SQRT3DMI=sqrt(3.0d0) / 2.0d0 )
2609
2610 REAL FLMIN
2611C-----------------------------------------------
2612C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
2613C FTPI=(4/3)*PI, TTPI=(2/3)*PI
2614C
2615C principal deviator of stress
2616C . . . . . . . . . . . . . . . . . . .
2617 MDEMI = -half
2618 ttpi = acos(mdemi)
2619 ftpi = two*ttpi
2620C minimum precision depending on the type real or double
2621 CALL floatmin(cs1,cs2,flmin)
2622 flm = two*sqrt(flmin)
2623 nindex3=0
2624 index3(1:nel) =0
2625 c1 = half*sqrt(twenty7)
2626 c2 = two*sqrt(third)
2627 DO nn = 1, nel
2628 cs1 = sig(nn,1)
2629 cs2 = sig(nn,2)
2630 cs3 = sig(nn,3)
2631 cs4 = sig(nn,4)
2632 cs5 = sig(nn,5)
2633 cs6 = sig(nn,6)
2634 pr = -(cs1+cs2+cs3) * third
2635 cs1 = cs1 + pr
2636 cs2 = cs2 + pr
2637 cs3 = cs3 + pr
2638 s4_2 = cs4*cs4
2639 s5_2 = cs5*cs5
2640 s6_2 = cs6*cs6
2641 aaa(nn)=s4_2 + s5_2 + s6_2 - cs1*cs2
2642 & - cs2*cs3 - cs1*cs3
2643 norminf(nn) = max(abs(cs1),abs(cs2),abs(cs3),
2644 & abs(cs4),abs(cs5),abs(cs6))
2645 norminf(nn) = em10*norminf(nn)
2646C cas racine triple
2647 aa = max(aaa(nn),norminf(nn))
2648 sq_aa = sqrt(aa)
2649C
2650 bb=cs1*s5_2 + cs2*s6_2
2651 & + cs3*s4_2 - cs1*cs2*cs3
2652 & - two*cs4*cs5*cs6
2653C
2654 cc=-c1/max(em10,sq_aa)*bb/max(em20,aa)
2655c CC=-SQRT(TWENTY7/MAX(EM20,AA))*BB*HALF/MAX(EM20,AA)
2656 cc= min(cc,one)
2657 cc= max(cc,-one)
2658 angp=acos(cc) * third
2659 dd=c2*sq_aa
2660
2661 cos_angp = cos(angp)
2662 sin_angp = sin(angp)
2663c STR1(NN)=DD*COS_ANGP
2664c STR2(NN)=DD*((-HALF*COS_ANGP) - (SIN_ANGP*SQRT3DMI))
2665c STR3(NN)=DD*((-HALF*COS_ANGP) + (SIN_ANGP*SQRT3DMI))
2666 str3(nn)=dd*cos_angp
2667 str1(nn)=dd*((-half*cos_angp) - (sin_angp*sqrt3dmi))
2668 str2(nn)=-(str1(nn)+str3(nn))
2669C
2670 valdp1(nn) = str1(nn)-pr
2671 valdp2(nn) = str2(nn)-pr
2672 valdp3(nn) = str3(nn)-pr
2673C precision reinforcement in simple compression ----
2674 IF(abs(str3(nn))>abs(str1(nn))
2675 & .AND.aaa(nn)>norminf(nn)) THEN
2676 aa=str1(nn)
2677 str1(nn)=str3(nn)
2678 str3(nn)=aa
2679 index3(nn) = 1
2680 nindex3 = 1
2681 ENDIF
2682
2683C . . . . . . . . . . .
2684C VECTEURS PROPRES
2685C . . . . . . . . . . .
2686 strmax= max(abs(str1(nn)),abs(str3(nn)))
2687 tol1(nn)= max(em20,flm*strmax*strmax)
2688 tol2(nn)=flm*strmax * third
2689 a(nn,1,1)=cs1-str1(nn)
2690 a(nn,2,2)=cs2-str1(nn)
2691 a(nn,3,3)=cs3-str1(nn)
2692 a(nn,1,2)=cs4
2693 a(nn,2,1)=cs4
2694 a(nn,2,3)=cs5
2695 a(nn,3,2)=cs5
2696 a(nn,1,3)=cs6
2697 a(nn,3,1)=cs6
2698C
2699 b(nn,1,1)=a(nn,2,1)*a(nn,3,2)-a(nn,3,1)
2700 . *a(nn,2,2)
2701 b(nn,1,2)=a(nn,2,2)*a(nn,3,3)-a(nn,3,2)
2702 . *a(nn,2,3)
2703 b(nn,1,3)=a(nn,2,3)*a(nn,3,1)-a(nn,3,3)
2704 . *a(nn,2,1)
2705 b(nn,2,1)=a(nn,3,1)*a(nn,1,2)-a(nn,1,1)
2706 . *a(nn,3,2)
2707 b(nn,2,2)=a(nn,3,2)*a(nn,1,3)-a(nn,1,2)
2708 . *a(nn,3,3)
2709 b(nn,2,3)=a(nn,3,3)*a(nn,1,1)-a(nn,1,3)
2710 . *a(nn,3,1)
2711 b(nn,3,1)=a(nn,1,1)*a(nn,2,2)-a(nn,2,1)
2712 . *a(nn,1,2)
2713 b(nn,3,2)=a(nn,1,2)*a(nn,2,3)-a(nn,2,2)
2714 . *a(nn,1,3)
2715 b(nn,3,3)=a(nn,1,3)*a(nn,2,1)-a(nn,2,3)
2716 . *a(nn,1,1)
2717 xmag1(nn)=b(nn,1,1)*b(nn,1,1)+b(nn,2,1)*b(nn,2,1)+
2718 . b(nn,3,1)*b(nn,3,1)
2719 xmag2(nn)=b(nn,1,2)*b(nn,1,2)+b(nn,2,2)*b(nn,2,2)+
2720 . b(nn,3,2)*b(nn,3,2)
2721 xmag3(nn)=b(nn,1,3)*b(nn,1,3)+b(nn,2,3)*b(nn,2,3)+
2722 . b(nn,3,3)*b(nn,3,3)
2723
2724 ENDDO
2725C
2726 nindex1 = 0
2727 nindex2 = 0
2728 DO nn = 1, nel
2729 xmaxv(nn)=max(xmag1(nn),xmag2(nn),xmag3(nn))
2730 IF(xmag1(nn)==xmaxv(nn)) THEN
2731 lmaxv(nn) = 1
2732 ELSEIF(xmag2(nn)==xmaxv(nn)) THEN
2733 lmaxv(nn) = 2
2734 ELSE
2735 lmaxv(nn) = 3
2736 ENDIF
2737 IF(aaa(nn)<norminf(nn)) THEN
2738 valdp1(nn) = sig(nn,1)
2739 valdp2(nn) = sig(nn,2)
2740 valdp3(nn) = sig(nn,3)
2741 v(nn,1,1) = one
2742 v(nn,2,1) = zero
2743 v(nn,3,1) = zero
2744 v(nn,1,2) = zero
2745 v(nn,2,2) = one
2746 v(nn,3,2) = zero
2747 ELSEIF(xmaxv(nn)>tol1(nn)*tol1(nn)) THEN
2748 nindex1 = nindex1 + 1
2749 index1(nindex1) = nn
2750 ELSE
2751 nindex2 = nindex2 + 1
2752 index2(nindex2) = nn
2753 ENDIF
2754 ENDDO
2755C
2756#include "vectorize.inc"
2757 DO n = 1, nindex1
2758 nn = index1(n)
2759 lmax = lmaxv(nn)
2760 xmaxinv = one/sqrt(xmaxv(nn))
2761 v(nn,1,1)=b(nn,1,lmax)*xmaxinv
2762 v(nn,2,1)=b(nn,2,lmax)*xmaxinv
2763 v(nn,3,1)=b(nn,3,lmax)*xmaxinv
2764 a(nn,1,1)=a(nn,1,1)+str1(nn)-str3(nn)
2765 a(nn,2,2)=a(nn,2,2)+str1(nn)-str3(nn)
2766 a(nn,3,3)=a(nn,3,3)+str1(nn)-str3(nn)
2767C
2768 b(nn,1,1)=a(nn,2,1)*v(nn,3,1)-a(nn,3,1)*v(nn,2,1)
2769 b(nn,1,2)=a(nn,2,2)*v(nn,3,1)-a(nn,3,2)*v(nn,2,1)
2770 b(nn,1,3)=a(nn,2,3)*v(nn,3,1)-a(nn,3,3)*v(nn,2,1)
2771 b(nn,2,1)=a(nn,3,1)*v(nn,1,1)-a(nn,1,1)*v(nn,3,1)
2772 b(nn,2,2)=a(nn,3,2)*v(nn,1,1)-a(nn,1,2)*v(nn,3,1)
2773 b(nn,2,3)=a(nn,3,3)*v(nn,1,1)-a(nn,1,3)*v(nn,3,1)
2774 b(nn,3,1)=a(nn,1,1)*v(nn,2,1)-a(nn,2,1)*v(nn,1,1)
2775 b(nn,3,2)=a(nn,1,2)*v(nn,2,1)-a(nn,2,2)*v(nn,1,1)
2776 b(nn,3,3)=a(nn,1,3)*v(nn,2,1)-a(nn,2,3)*v(nn,1,1)
2777 xmag1(nn)=b(nn,1,1)*b(nn,1,1)+b(nn,2,1)*b(nn,2,1)+
2778 . b(nn,3,1)*b(nn,3,1)
2779 xmag2(nn)=b(nn,1,2)*b(nn,1,2)+b(nn,2,2)*b(nn,2,2)+
2780 . b(nn,3,2)*b(nn,3,2)
2781 xmag3(nn)=b(nn,1,3)*b(nn,1,3)+b(nn,2,3)*b(nn,2,3)+
2782 . b(nn,3,3)*b(nn,3,3)
2783C
2784 xmaxv(nn)=max(xmag1(nn),xmag2(nn),xmag3(nn))
2785 ENDDO
2786C
2787#include "vectorize.inc"
2788 DO n = 1, nindex1
2789 nn = index1(n)
2790 IF(xmag1(nn)==xmaxv(nn)) THEN
2791 lmaxv(nn) = 1
2792 ELSEIF(xmag2(nn)==xmaxv(nn)) THEN
2793 lmaxv(nn) = 2
2794 ELSE
2795 lmaxv(nn) = 3
2796 ENDIF
2797C
2798 vmag=sqrt(v(nn,1,1)**2+v(nn,2,1)**2)
2799 lmax = lmaxv(nn)
2800 xmaxv(nn)=sqrt(xmaxv(nn))
2801 IF(xmaxv(nn)>tol2(nn))THEN
2802 xmaxinv = one/xmaxv(nn)
2803 v(nn,1,3)=b(nn,1,lmax)*xmaxinv
2804 v(nn,2,3)=b(nn,2,lmax)*xmaxinv
2805 v(nn,3,3)=b(nn,3,lmax)*xmaxinv
2806 v(nn,1,2)=v(nn,2,3)*v(nn,3,1)-v(nn,2,1)*v(nn,3,3)
2807 v(nn,2,2)=v(nn,3,3)*v(nn,1,1)-v(nn,3,1)*v(nn,1,3)
2808 v(nn,3,2)=v(nn,1,3)*v(nn,2,1)-v(nn,1,1)*v(nn,2,3)
2809 c1 = v(nn,1,2)*v(nn,1,2)+v(nn,2,2)*v(nn,2,2)+
2810 & v(nn,3,2)*v(nn,3,2)
2811 vmag=one/sqrt(c1)
2812 v(nn,1,2)=v(nn,1,2)*vmag
2813 v(nn,2,2)=v(nn,2,2)*vmag
2814 v(nn,3,2)=v(nn,3,2)*vmag
2815 ELSEIF(vmag>tol2(nn))THEN
2816 v(nn,1,2)=-v(nn,2,1)/vmag
2817 v(nn,2,2)=v(nn,1,1)/vmag
2818 v(nn,3,2)=zero
2819 ELSE
2820 v(nn,1,2)=one
2821 v(nn,2,2)=zero
2822 v(nn,3,2)=zero
2823 ENDIF
2824 ENDDO
2825C . . . . . . . . . . . . .
2826C SOLUTION DOUBLE
2827C . . . . . . . . . . . . .
2828#include "vectorize.inc"
2829 DO n = 1, nindex2
2830 nn = index2(n)
2831 xmag1(nn)=a(nn,1,1)*a(nn,1,1)+a(nn,2,1)*a(nn,2,1)
2832 xmag2(nn)=a(nn,1,2)*a(nn,1,2)+a(nn,2,2)*a(nn,2,2)
2833 xmag3(nn)=a(nn,1,3)*a(nn,1,3)+a(nn,2,3)*a(nn,2,3)
2834C
2835 xmaxv(nn)=max(xmag1(nn),xmag2(nn),xmag3(nn))
2836 ENDDO
2837C
2838#include "vectorize.inc"
2839 DO n = 1, nindex2
2840 nn = index2(n)
2841 IF(xmag1(nn)==xmaxv(nn)) THEN
2842 lmaxv(nn) = 1
2843 ELSEIF(xmag2(nn)==xmaxv(nn)) THEN
2844 lmaxv(nn) = 2
2845 ELSE
2846 lmaxv(nn) = 3
2847 ENDIF
2848C
2849 lmax = lmaxv(nn)
2850 xmaxv(nn)=sqrt(xmaxv(nn))
2851 IF(max(abs(a(nn,3,1)),abs(a(nn,3,2)),abs(a(nn,3,3)))
2852 & <tol2(nn))THEN
2853 xmaxinv = one/xmaxv(nn)
2854 v(nn,1,1)= zero
2855 v(nn,2,1)= zero
2856 v(nn,3,1)= one
2857 v(nn,1,2)=-a(nn,2,lmax)*xmaxinv
2858 v(nn,2,2)= a(nn,1,lmax)*xmaxinv
2859 v(nn,3,2)= zero
2860C
2861 ELSEIF(xmaxv(nn)>tol2(nn))THEN
2862 xmaxinv = one/xmaxv(nn)
2863 v(nn,1,1)=-a(nn,2,lmax)*xmaxinv
2864 v(nn,2,1)= a(nn,1,lmax)*xmaxinv
2865 v(nn,3,1)= zero
2866 v(nn,1,2)=-a(nn,3,lmax)*v(nn,2,1)
2867 v(nn,2,2)= a(nn,3,lmax)*v(nn,1,1)
2868 v(nn,3,2)= a(nn,1,lmax)*v(nn,2,1)-a(nn,2,lmax)*v(nn,1,1)
2869 c1 = v(nn,1,2)*v(nn,1,2)+v(nn,2,2)*v(nn,2,2)+
2870 & v(nn,3,2)*v(nn,3,2)
2871 vmag=one/sqrt(c1)
2872 v(nn,1,2)=v(nn,1,2)*vmag
2873 v(nn,2,2)=v(nn,2,2)*vmag
2874 v(nn,3,2)=v(nn,3,2)*vmag
2875 ELSE
2876 valdp1(nn) = sig(nn,1)
2877 valdp2(nn) = sig(nn,2)
2878 valdp3(nn) = sig(nn,3)
2879 v(nn,1,1) = one
2880 v(nn,2,1) = zero
2881 v(nn,3,1) = zero
2882 v(nn,1,2) = zero
2883 v(nn,2,2) = one
2884 v(nn,3,2) = zero
2885 index3(nn) = 0
2886 ENDIF
2887 ENDDO
2888C
2889 DO nn = 1, nel
2890 vecdp1(nn)=v(nn,1,1)
2891 vecdp2(nn)=v(nn,2,1)
2892 vecdp3(nn)=v(nn,3,1)
2893 vecdp4(nn)=v(nn,1,2)
2894 vecdp5(nn)=v(nn,2,2)
2895 vecdp6(nn)=v(nn,3,2)
2896 vecdp7(nn)=vecdp2(nn)*vecdp6(nn)-vecdp3(nn)*vecdp5(nn)
2897 vecdp8(nn)=vecdp3(nn)*vecdp4(nn)-vecdp1(nn)*vecdp6(nn)
2898 vecdp9(nn)=vecdp1(nn)*vecdp5(nn)-vecdp2(nn)*vecdp4(nn)
2899 ENDDO
2900C
2901 DO nn = 1, nel
2902 val(nn,1)=valdp1(nn)
2903 val(nn,2)=valdp2(nn)
2904 val(nn,3)=valdp3(nn)
2905 vec(nn,1)=vecdp1(nn)
2906 vec(nn,2)=vecdp2(nn)
2907 vec(nn,3)=vecdp3(nn)
2908 vec(nn,4)=vecdp4(nn)
2909 vec(nn,5)=vecdp5(nn)
2910 vec(nn,6)=vecdp6(nn)
2911 vec(nn,7)=vecdp7(nn)
2912 vec(nn,8)=vecdp8(nn)
2913 vec(nn,9)=vecdp9(nn)
2914 ENDDO
2915C Rewriting to bypass problem on itanium2 comp 9. + latency = 16
2916 IF (nindex3 > 0) THEN
2917 DO nn = 1, nel
2918 IF(index3(nn) == 1) THEN
2919C Structilis com temporary table instead of temporary scalars
2920 s1=vec(nn,7)
2921 s2=vec(nn,8)
2922 s3=vec(nn,9)
2923 vec(nn,7)=vec(nn,1)
2924 vec(nn,8)=vec(nn,2)
2925 vec(nn,9)=vec(nn,3)
2926 vec(nn,1)=-s1
2927 vec(nn,2)=-s2
2928 vec(nn,3)=-s3
2929 ENDIF
2930 ENDDO
2931 END IF !(NINDEX3 > 0) THEN
2932 RETURN
2933 END
2934
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
floatmin
void floatmin(int *a, int *b, float *flm)
Definition precision.c:71
dreh
subroutine dreh(b, tr, ii, jj, ken)
Definition sigeps33.F:741
valpvecdp
subroutine valpvecdp(sig, val, vec, nel)
Definition sigeps33.F:1112
epsf2u
subroutine epsf2u(nel, fxx, fxy, fxz, fyx, fyy, fyz, fzx, fzy, fzz, uxx, uyy, uzz, uxy, uyz, uxz)
Definition sigeps33.F:1448
jacobiew
subroutine jacobiew(a, n, ew, ev, nrot)
Definition sigeps33.F:609
valpvecdp_v
subroutine valpvecdp_v(sig, val, vec, nel)
Definition sigeps33.F:2235
sigeps33
subroutine sigeps33(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off)
Definition sigeps33.F:46
valpvec
subroutine valpvec(sig, val, vec, nel)
Definition sigeps33.F:795
valpvecop
subroutine valpvecop(sig, val, vec, nel)
Definition sigeps33.F:1540
valpvecop_v
subroutine valpvecop_v(sig, val, vec, nel)
Definition sigeps33.F:2569
valpvec_v
subroutine valpvec_v(sig, val, vec, nel)
Definition sigeps33.F:1904
sigeps69
subroutine sigeps69(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, wxxdt, wyydt, wzzdt, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, et, ihet, nuvarv, uvarv, offg, epsth3, iexpan)
Definition sigeps69.F:46