OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps101.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/.
23C polypropylene - Bouvard law
24!||====================================================================
25!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
26!||--- called by ------------------------------------------------------
27!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
28!||--- calls -----------------------------------------------------
29!|| at_tpu ../engine/source/materials/mat/mat101/sigeps101.F
30!|| atxb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
31!|| axb_tpu ../engine/source/materials/mat/mat101/sigeps101.f
32!|| axbt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
33!|| det_tpu ../engine/source/materials/mat/mat101/sigeps101.F
34!|| dev_tpu ../engine/source/materials/mat/mat101/sigeps101.F
35!|| eigsrt ../engine/source/materials/mat/mat101/sigeps101.F
36!|| emod_tpu ../engine/source/materials/mat/mat101/sigeps101.F
37!|| err_tpu ../engine/source/materials/mat/mat101/sigeps101.F
38!|| finter ../engine/source/tools/curve/finter.F
39!|| inv_tpu ../engine/source/materials/mat/mat101/sigeps101.F
40!|| mat3x3tovec6x1 ../engine/source/materials/mat/mat101/sigeps101.F
41!|| qaqt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
42!|| qtaq_tpu ../engine/source/materials/mat/mat101/sigeps101.F
43!|| skw_tpu ../engine/source/materials/mat/mat101/sigeps101.F
44!|| sp_tpu ../engine/source/materials/mat/mat101/sigeps101.F
45!|| sst_tpu ../engine/source/materials/mat/mat101/sigeps101.F
46!|| st_tpu ../engine/source/materials/mat/mat101/sigeps101.F
47!|| sym_tpu ../engine/source/materials/mat/mat101/sigeps101.F
48!|| tr_tpu ../engine/source/materials/mat/mat101/sigeps101.F
49!|| valprop_tpu ../engine/source/materials/mat/mat101/sigeps101.F
50!|| vec6x1tomat3x3 ../engine/source/materials/mat/mat101/sigeps101.F
51!||====================================================================
52 SUBROUTINE sigeps101(
53 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC ,
54 2 NPF ,TF , TIME , TIMESTEP, UPARAM,
55 3 RHO0 , RHO ,VOLUME , EINT , NGL ,
56 5 DEPSXX , DEPSYY , DEPSZZ , DEPSXY, DEPSYZ, DEPSZX,
57 6 EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
58 7 SIGOXX , SIGOYY , SIGOZZ , SIGOXY, SIGOYZ, SIGOZX,
59 8 SIGNXX , SIGNYY , SIGNZZ , SIGNXY, SIGNYZ, SIGNZX,
60 9 SIGVXX , SIGVYY , SIGVZZ , SIGVXY, SIGVYZ, SIGVZX,
61 C SOUNDSP, VISCMAX, UVAR , OFF , ISMSTR, ET ,
62 D IHET ,OFFG , EPSTH, IEXPAN ,TEMPEL,
63 1 FPSXX , FPSXY , FPSXZ , FPSYX ,
64 2 FPSYY ,FPSYZ , FPSZX , FPSZY , FPSZZ ,
65 3 UPSXX ,UPSYY , UPSZZ , UPSXY , UPSYZ ,
66 4 UPSXZ )
67C
68C-----------------------------------------------
69C I M P L I C I T T Y P E S
70C-----------------------------------------------
71#include "implicit_f.inc"
72C----------------------------------------------------------------
73C I N P U T A R G U M E N T S
74C----------------------------------------------------------------
75 INTEGER NEL, NVARF, NUPARAM, NUVAR,ISMSTR,NGL(*),IHET,IEXPAN
76 my_real
77 . TIME , TIMESTEP , UPARAM(NUPARAM),UPARAMF(NUPARAM),
78 . RHO (NEL), RHO0 (NEL), VOLUME(NEL), EINT(NEL),
79 . DEPSXX(NEL), DEPSYY(NEL), DEPSZZ(NEL),
80 . DEPSXY(NEL), DEPSYZ(NEL), DEPSZX(NEL),
81 . EPSXX (NEL), EPSYY (NEL), EPSZZ (NEL),
82 . EPSXY (NEL), EPSYZ (NEL), EPSZX (NEL),
83 . SIGOXX(NEL), SIGOYY(NEL), SIGOZZ(NEL),
84 . SIGOXY(NEL), SIGOYZ(NEL), SIGOZX(NEL),OFFG(NEL),
85 . EPSTH(NEL) ,TEMPEL(NEL),
86 . FPSXX(NEL),FPSYY(NEL),FPSZZ(NEL),
87 . FPSXY(NEL),FPSYX(NEL),FPSXZ(NEL),
88 . FPSZX(NEL),FPSYZ(NEL),FPSZY(NEL),
89 . upsxx(nel),upsyy(nel),upszz(nel),
90 . upsxy(nel),upsyz(nel),upsxz(nel)
91C-----------------------------------------------
92
93C----------------------------------------------------------------
94C O U T P U T A R G U M E N T S
95C----------------------------------------------------------------
96 my_real
97 . signxx(nel), signyy(nel), signzz(nel),
98 . signxy(nel), signyz(nel), signzx(nel),
99 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
100 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
101 . soundsp(nel), viscmax(nel), et(nel)
102C----------------------------------------------------------------
103C I N P U T O U T P U T A R G U M E N T S
104C----------------------------------------------------------------
105 my_real
106 . uvar(nel,nuvar), off(nel)
107C----------------------------------------------------------------
108C VARIABLES FOR FUNCTION INTERPOLATION
109C----------------------------------------------------------------
110 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
111 my_real FINTER,FINTTE,TF(*),FINT2V
112 EXTERNAL FINTER,FINTTE
113C----------------------------------------------------------------
114C L O C A L V A R I B L E S
115C----------------------------------------------------------------
116 INTEGER I,J,KK,LL,N
117 INTEGER IE,NROT,IER,NDI,NSHR,NTENS,SIGN
118
119 my_real
120 . et1,et2,et3,g,rbulk,aa,cc,sb,
121C
122 . bi1(nel),bi2(nel),jdet(nel),i1(nel),ip1(nel),gammaold(nel),stiff(nel),
123 . temp(nel),
124C
125 . nb(nel,3),sninv(nel,3,3)
126C----------------------------------------------------------------
127
128 LOGICAL DEBUG, DEBUG1, DEBUG2, DEBUG3, DEBUG4, DEBUG5
129
130 my_real
131 . timea(2),dtime,f(3,3), f0(3,3), u0(3,3), f1(3,3), u1(3,3),
132 . tensor1(3,3), tensor2(3,3), tensor3(3,3), tensor4(3,3),
133 . tensor5(3,3),
134 . tens0(3,3), tens1(3,3),
135 . rdfgrd(3,3), vgrad(3,3),
136 . d_umat(3,3), w_umat(3,3),
137 . rote0(3,3), rote1(3,3), rot0(3,3), rot1(3,3),
138 . betam(3,3), betam0(3,3), alpham(3,3),
139 . dev_alpham(3,3), dev_xm(3,3),
140 . taum(3,3), sigma_q(3,3),
141 . fvm(3,3), fvm_inv(3,3), bp(3,3),
142 . ftheta(3,3), ftheta_inv(3,3),
143 . fe0(3,3), fi_inv(3,3), ce0(3,3), ue0(3,3),
144 . ue0_inv(3,3), ee0(3,3), xm(3,3), sigmam(3,3),
145 . dev_meff(3,3), dv(3,3), temp1(3,3), xmeff(3,3),
146 . temp1_inv(3,3), temp2(3,3), temp2_inv(3,3),
147 . exp_dvdt(3,3), fvm_iter(3,3),
148 . fe1(3,3), ce1(3,3), ue1(3,3),
149 . ue1_inv(3,3), ee1(3,3),
150 . vce0(3,3), dce0(3),
151 . vce1(3,3), dce1(3),
152 . a(3,3),
153c
154 . fv(6), sigma(6),
155 . beta(6),
156 . dev_meff_vec(6), nv_vec(6),
157 . dv_vec(6), tau(6),
158 . fv_iter(6), xm_vec(6), d_umat_vec(6),
159c
160 . j1, t1, t2, t3,
161 . deref, de0, de, dpoiss,
162 . dedot_ref, dve1, dve2,
163 . gamv_ref, dq, dv1, dc3, dc4, calphak1,
164 . calphak2, dc5, dc6, dc7, dc8, dc9, dc10,
165 . dc11, dc12, dc13, dc14, dc1, dc2,
166 . theta0, beta0, btheta, factor, cv, rho_p,
167 . theta, thetadot, dtheta_iter, thetadotiter,
168 . thetaiter, thetag,
169 . d_theta_opt, thetai,
170 . dmu, dk, alphap, dm, dy0,
171 . ckappa1, ckappa2 ,
172 . h0, des1_0, destar_0, destar_s, g0,
173 . h1, des2_0, des2_s, tr_bp,
174 . des10, des20, destar0, dlambdap_bar,
175 . dmu_r, dlambda_l, drs1,
176 . twomu0, bulka0,
177 . des1, destar, des2, dkappa1, dkappa2,
178 . gamv, gamveq, dkappa, tr_alpha,
179 . det_f0, det_f1, det, det_fe,
180 . dmu_b, ksi_b, det_beta,
181 . eqstress, pressure, norm_devmeff, tr_m,
182 . dgamv_iter, det_vce0,
183 . gamv_iter, det_vce1, d1, d2, d3,
184 . de1, de2, de3,
185 . tr_ee1, tr_ee0, coeff1_a, coeff1_b,
186 . coeff1, coeff2, coeff3, exchange, gamv0, rot, ps
187
188 my_real
189 . tseven, small, big, ds32, ds3, ds13, dkb,
190 . dr
191
192 my_real
193 . pident(6), factotwo(6), factohalf(6)
194 my_real
195 . pmident(3,3)
196C
197 DATA pmident /1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0,
198 & 0.d0, 0.d0, 1.d0/
199 DATA pident /1.d0, 1.d0, 1.d0, 0.d0, 0.d0, 0.d0/
200 DATA factotwo /1.d0, 1.d0, 1.d0, 2.d0, 2.d0, 2.d0/
201 DATA factohalf /1.d0, 1.d0, 1.d0, 0.5d0, 0.5d0, 0.5d0/
202C-----------------------------------------------
203C
204 ndi= 3
205C-----------------------------------------------
206C TIMESTEP
207
208
209 dtime = timestep
210 tseven = two*ten+seven
211 small = em6
212 big = ep16
213
214 ds32 = sqrt(three/two)
215 ds3 = sqrt(three)
216 ds13 = sqrt(one/three)
217C---- MATERIAL CONSTANTS
218C
219 dkb = 1.3806*em30 ! boltzmann constant
220 !DR = 8.31*EM3
221C
222C---- DEBUG
223C
224 debug = .false. !STEP0
225 debug1= .false. !STEP1 DATA
226 debug2= .false. !STEP1 BRANCHA_0
227 debug3= .false. !STEP1 VISCOPLASTICITY
228 debug4= .false. !STEP1 BRANCHB_1
229 debug5= .false. !STEP1 BRANCHA_1 & STRESS
230 deref = uparam(1)
231 de0 = uparam(2)
232 dpoiss = uparam(3)
233 dve1 = uparam(4)
234 dve2 = uparam(5)
235 dedot_ref = uparam(6)
236 gamv_ref = uparam(7)
237 alphap = uparam(8)
238 dq = uparam(9)
239 dv1 = uparam(10)
240 dm = uparam(11)
241 dc3 = uparam(12)
242 dc4 = uparam(13)
243 calphak1 = uparam(14)
244 calphak2 = uparam(15)
245 h0 = uparam(16)
246 des1_0 = uparam(17)
247 dc5 = uparam(18)
248 dc6 = uparam(19)
249 dc7 = uparam(20)
250 dc8 = uparam(21)
251 dc9 = uparam(22)
252 dc10 = uparam(23)
253 h1 = uparam(24)
254 des2_0 = uparam(25)
255 dc11 = uparam(26)
256 dc12 = uparam(27)
257 dc13 = uparam(28)
258 dc14 = uparam(29)
259 dc1 = uparam(30)
260 dc2 = uparam(31)
261 dlambda_l = uparam(32)
262 rho_p = uparam(33)
263 cv = uparam(34)
264 theta0 = uparam(35)
265 beta0 = uparam(36)
266 !BTHETA = UPARAM(37)
267 factor = uparam(38)
268 d_theta_opt = uparam(39)
269 thetai = uparam(40)
270 thetag = uparam(41)
271 dr = uparam(42)
272
273
274c------------------------------------------------------
275
276C-----------------------------------------------
277C RECUPERER CHAMP TEMPERATURE
278C-----------------------------------------------
279 temp(1:nel) = tempel(1:nel)
280C-----------------------------------------------
281C USER VARIABLES INITIALIZATION
282C-----------------------------------------------
283
284
285 IF ( (time == zero) ) THEN
286
287C
288C LOOP THROUGH ALL BLOCKS
289C
290 DO ie = 1, nel
291
292C
293C------ TEMPERATURE DEPENDENCE 1
294C
295 destar_0 = dc5*(thetai-theta0)+ dc6
296 dmu_r = dc1*(thetai-theta0)+ dc2
297C
298C-------- COMPUTE DMU_B
299C
300 dmu_b = dmu_r
301C
302C-------- INITIALIZE STATE VARIABLES
303C
304 uvar(ie, 1) = one !FV11
305 uvar(ie, 2) = one !FV22
306 uvar(ie, 3) = one !FV33
307 uvar(ie, 4) = zero !FV12
308 uvar(ie, 5) = zero !FV23
309 uvar(ie, 6) = zero !FV13
310 uvar(ie, 7) = des1_0 !DES1
311 uvar(ie, 8) = destar_0 !DESTAR
312 uvar(ie, 9) = des2_0 !DES2
313 uvar(ie,10) = one !beta11
314 uvar(ie,11) = one !BETA22
315 uvar(ie,12) = one !BETA33
316 uvar(ie,13) = zero !BETA12
317 uvar(ie,14) = zero !BETA23
318 uvar(ie,15) = zero !BETA13
319 uvar(ie,16) = dmu_b !DMU_B
320 uvar(ie,17) = zero !GAMV
321 uvar(ie,18) = zero !DGAMV
322 uvar(ie,19) = zero !GAMVEQ
323 uvar(ie,20) = thetai
324 IF (abs(d_theta_opt-one) < em06) THEN
325 uvar(ie,20) = temp(ie)
326 ENDIF
327 uvar(ie,21) = zero !THETADOT (ADIABATIC)
328
329 uvar(ie,22) = one !F011
330 uvar(ie,23) = one !F022
331 uvar(ie,24) = one !F033
332 uvar(ie,25) = zero !F012
333 uvar(ie,26) = zero !F021
334 uvar(ie,27) = zero !F013
335 uvar(ie,28) = zero !F031
336 uvar(ie,29) = zero !F023
337 uvar(ie,30) = zero !F032
338
339 uvar(ie,31) = one !U011
340 uvar(ie,32) = one !U022
341 uvar(ie,33) = one !U033
342 uvar(ie,34) = zero !U012
343 uvar(ie,35) = zero !U021
344 uvar(ie,36) = zero !U013
345 uvar(ie,37) = zero !U031
346 uvar(ie,38) = zero !U023
347 uvar(ie,39) = zero !U032
348 uvar(ie,40) = zero !E133 true strain xx
349 uvar(ie,41) = zero !E133 true strain yy
350 uvar(ie,42) = zero !E133 true strain zz
351
352C
353C------ INTERNAL STATE VARIABLES AT BEGINNING OF STEP
354
355 fv(1) = uvar(ie, 1)
356 fv(2) = uvar(ie, 2)
357 fv(3) = uvar(ie, 3)
358 fv(4) = uvar(ie, 4)
359 fv(5) = uvar(ie, 5)
360 fv(6) = uvar(ie, 6)
361 theta = uvar(ie, 20)
362C
363C----- setup F1; U1
364C
365 f1(1,1) = fpsxx(ie)
366 f1(2,2) = fpsyy(ie)
367 f1(3,3) = fpszz(ie)
368 f1(1,2) = fpsxy(ie)
369 f1(2,1) = fpsyx(ie)
370 f1(1,3) = fpsxz(ie)
371 f1(3,1) = fpszx(ie)
372 f1(2,3) = fpsyz(ie)
373 f1(3,2) = fpszy(ie)
374
375 u1(1,1) = upsxx(ie)
376 u1(2,2) = upsyy(ie)
377 u1(3,3) = upszz(ie)
378 u1(1,2) = upsxy(ie)
379 u1(2,1) = upsxy(ie)
380 u1(1,3) = upsxz(ie)
381 u1(3,1) = upsxz(ie)
382 u1(2,3) = upsyz(ie)
383 u1(3,2) = upsyz(ie)
384C---- COMPUTE DROT_ABA
385C
386C----- --- ROTATION TENSORS FROM F=RU, I.E., R=F*U^(-1)
387C
388 CALL st_tpu(tensor1, zero, 3*3)
389 CALL inv_tpu(u1, tensor1, det)
390 CALL axb_tpu(f1, tensor1, rot1, 3)
391
392 !--------------------------------------------------------------
393 !-------- COMPUTE ELASTIC STRESSES
394 !--------------------------------------------------------------
395 !-------- TEMPERATURE DEPENDENCE OF PARAMETERS
396 !---(STRAIN RATE DEPENDENCE IS NEGLECTED FOR TIME=0)
397 !
398
399 CALL emod_tpu(deref,de0,theta,theta0,dve1,dve2,
400 & em04,dedot_ref,de)
401
402 dmu = de / (two+two*dpoiss)
403 dk = two*dmu*(one+dpoiss)/(three*(one-two*dpoiss))
404
405 !CALCULATE SOME CONSTANTS
406
407 twomu0 = two*dmu
408 bulka0 = dk - (two/3.d0)*dmu
409
410 !COMPUTE: FE->CE->UE->RE->EE
411 CALL st_tpu(fvm, zero, 3*3)
412 CALL st_tpu(fvm_inv,zero, 3*3)
413 CALL st_tpu(tens0, zero, 3*3)
414 CALL st_tpu(fe1, zero, 3*3)
415 CALL st_tpu(ce1, zero, 3*3)
416 CALL st_tpu(ue1, zero, 3*3)
417 CALL st_tpu(ue1_inv,zero, 3*3)
418 CALL st_tpu(rote1, zero, 3*3)
419 CALL st_tpu(ee1, zero, 3*3)
420
421
422 !DETERMINATION OF FE AND CE
423 CALL vec6x1tomat3x3(fv, fvm,3)
424 CALL inv_tpu(fvm, fvm_inv, det)
425 CALL axb_tpu(f1, fvm_inv, fe1, 3)
426 CALL atxb_tpu(fe1, fe1, ce1, 3)
427C
428 DO i=1,ndi
429 DO j=1,ndi
430 a(i,j)= ce1(i,j)
431 ENDDO
432 ENDDO
433C
434 CALL valprop_tpu(a, 3, 3, dce1, vce1, nrot, ier)
435C
436 CALL eigsrt(dce1, vce1, 3, 3)
437C
438 IF (ier==1) THEN
439 CALL err_tpu('ERROR IN EIG VAL')
440 ENDIF
441 !EIGENVALUES (AND ASSOC EIGENVECTORS) ARE ORDERED FROM
442 !LARGER TO SMALLER.
443 !REDEFINE AXIS(2) TO BE THE LARGEST IN ORDER TO IMPROVE
444 !ACCURACY IN THE
445 !CALCULATION OF THE ESHELBY TENSOR.
446 !IF DET(B)<0 MEANS THAT THE SYSTEM IS LEFT HANDED. IT IS MADE
447 !RIGHT
448 !HANDED BY EXCHANGING 1 AND 2.
449 sign=-one
450 CALL det_tpu(vce1, det_vce1)
451 IF(det_vce1 <= zero) sign=one
452 DO i=1,3
453 exchange=vce1(i,1)
454 vce1(i,1)=vce1(i,2)
455 vce1(i,2)=exchange*sign
456 ENDDO
457 exchange=dce1(1)
458 dce1(1)=dce1(2)
459 dce1(2)=exchange
460 !COMPUTE UE1
461 d1=sqrt(dce1(1))
462 d2=sqrt(dce1(2))
463 d3=sqrt(dce1(3))
464 de1=log(sqrt(dce1(1)))
465 de2=log(sqrt(dce1(2)))
466 de3=log(sqrt(dce1(3)))
467 DO i=1,3
468 DO j=1,3
469 ue1(i,j)=d1*vce1(i,1)*vce1(j,1)
470 & +d2*vce1(i,2)*vce1(j,2)+d3*vce1(i,3)*vce1(j,3)
471 END DO
472 END DO
473 DO i=1,3
474 DO j=1,3
475 ee1(i,j)=de1*vce1(i,1)*vce1(j,1)
476 & +de2*vce1(i,2)*vce1(j,2)+de3*vce1(i,3)*vce1(j,3)
477 END DO
478 END DO
479C
480 CALL inv_tpu(ue1, ue1_inv, det)
481 CALL axb_tpu(fe1, ue1_inv, rote1, 3)
482 CALL det_tpu(fe1, det_fe)
483C
484 !DETERMINATION OF TAU
485C
486 CALL st_tpu(xm, zero, 3*3)
487 CALL st_tpu(taum, zero, 3*3)
488C
489 CALL tr_tpu(ee1, tr_ee1, 3)
490 coeff1_a = twomu0
491 coeff1_b = bulka0*tr_ee1
492 CALL at_tpu(coeff1_a, ee1, coeff1_b, pmident, xm, 9)
493 CALL qaqt_tpu(rote1, xm, taum, 3)
494
495C
496 !DETERMINATION OF SIGMA
497C
498 CALL st_tpu(sigmam, zero, 3*3)
499 CALL st_tpu(sigma_q, zero, 3*3)
500C
501 CALL sst_tpu(one/det_fe, taum, sigmam, 3*3)
502C
503C------- DETERMINATION OF SIGMA (IN THE ROTATED CONFIGURATION)
504 CALL qtaq_tpu(rot1, sigmam, sigma_q, 3)
505C
506C------- SIGMA_Q(I,J)
507 uvar(ie,22) = f1(1,1) !F011
508 uvar(ie,23) = f1(2,2) !F022
509 uvar(ie,24) = f1(3,3) !F033
510 uvar(ie,25) = f1(1,2) !F012
511 uvar(ie,26) = f1(2,1) !F021
512 uvar(ie,27) = f1(1,3) !F013
513 uvar(ie,28) = f1(3,1) !F031
514 uvar(ie,29) = f1(2,3) !F023
515 uvar(ie,30) = f1(3,2) !F032
516 uvar(ie,31) = u1(1,1) !U011
517 uvar(ie,32) = u1(2,2) !U022
518 uvar(ie,33) = u1(3,3) !U033
519 uvar(ie,34) = u1(1,2) !U012
520 uvar(ie,35) = u1(2,1) !U021
521 uvar(ie,36) = u1(1,3) !U013
522 uvar(ie,37) = u1(3,1) !U031
523 uvar(ie,38) = u1(2,3) !U023
524 uvar(ie,39) = u1(3,2) !U032
525 uvar(ie,40) = log(u1(1,1)) !E111
526 uvar(ie,41) = log(u1(2,2)) !E122
527 uvar(ie,42) = log(u1(3,3)) !E133
528 signxx(ie) = sigoxx(ie)
529 signyy(ie) = sigoyy(ie)
530 signzz(ie) = sigozz(ie)
531 signxy(ie) = sigoxy(ie)
532 signyz(ie) = sigoyz(ie)
533 signzx(ie) = sigozx(ie)
534
535 soundsp(ie) = sqrt((dk+4*dmu/3)/rho0(ie))
536 viscmax(ie) = zero
537
538 ENDDO
539C
540 RETURN
541C
542 ENDIF ! ( (TIME == ZERO) )
543C
544
545C
546C---------------------------------------------
547C---- EVOLVE MATERIAL STATE - COMPUTE STRESSES
548C---------------------------------------------
549C
550C---- LOOP OVER ELEMENTS/IPS TO UPDATE STRESS
551C
552
553 DO 20 ie = 1, nel
554C
555 !INTERNAL STATE VARIABLES AT BEGINNING OF STEP
556 fv(1) = uvar(ie, 1)
557 fv(2) = uvar(ie, 2)
558 fv(3) = uvar(ie, 3)
559 fv(4) = uvar(ie, 4)
560 fv(5) = uvar(ie, 5)
561 fv(6) = uvar(ie, 6)
562 des10 = uvar(ie, 7)
563 destar0 = uvar(ie, 8)
564 des20 = uvar(ie, 9)
565 beta(1) = uvar(ie,10)
566 beta(2) = uvar(ie,11)
567 beta(3) = uvar(ie,12)
568 beta(4) = uvar(ie,13)
569 beta(5) = uvar(ie,14)
570 beta(6) = uvar(ie,15)
571 dmu_b = uvar(ie,16)
572 gamv = uvar(ie, 17)
573 gamveq = uvar(ie, 19)
574 theta = uvar(ie, 20)
575 thetadot = uvar(ie, 21)
576 IF (abs(d_theta_opt-one) < small) theta=temp(ie)
577C
578C----- SETUP F0; U0; F1; U1;
579 CALL st_tpu(f0, zero, 3*3)
580 CALL st_tpu(u0, zero, 3*3)
581 CALL st_tpu(f1, zero, 3*3)
582 CALL st_tpu(u1, zero, 3*3)
583C
584
585C
586C----- setup F0; U0, F1, U1
587C
588 f0(1,1) = uvar(ie,22)
589 f0(2,2) = uvar(ie,23)
590 f0(3,3) = uvar(ie,24)
591 f0(1,2) = uvar(ie,25)
592 f0(2,1) = uvar(ie,26)
593 f0(1,3) = uvar(ie,27)
594 f0(3,1) = uvar(ie,28)
595 f0(2,3) = uvar(ie,29)
596 f0(3,2) = uvar(ie,30)
597
598 u0(1,1) = uvar(ie,31)
599 u0(2,2) = uvar(ie,32)
600 u0(3,3) = uvar(ie,33)
601 u0(1,2) = uvar(ie,34)
602 u0(2,1) = uvar(ie,35)
603 u0(1,3) = uvar(ie,36)
604 u0(3,1) = uvar(ie,37)
605 u0(2,3) = uvar(ie,38)
606 u0(3,2) = uvar(ie,39)
607
608 f1(1,1) = fpsxx(ie)
609 f1(2,2) = fpsyy(ie)
610 f1(3,3) = fpszz(ie)
611 f1(1,2) = fpsxy(ie)
612 f1(2,1) = fpsyx(ie)
613 f1(1,3) = fpsxz(ie)
614 f1(3,1) = fpszx(ie)
615 f1(2,3) = fpsyz(ie)
616 f1(3,2) = fpszy(ie)
617
618 u1(1,1) = upsxx(ie)
619 u1(2,2) = upsyy(ie)
620 u1(3,3) = upszz(ie)
621 u1(1,2) = upsxy(ie)
622 u1(2,1) = upsxy(ie)
623 u1(1,3) = upsxz(ie)
624 u1(3,1) = upsxz(ie)
625 u1(2,3) = upsyz(ie)
626 u1(3,2) = upsyz(ie)
627! C
628C------ !COMPUTE L,D,W,DROT FROM F_N & F_(N+1)
629C
630C-------!RELATIVE DEFORMATION GRADIENT: F=F_(N+1)*F_N^{-1}
631C
632 CALL st_tpu(tensor1, zero, 3*3)
633 CALL st_tpu(rdfgrd, zero, 3*3)
634 CALL st_tpu(tensor2, zero, 3*3)
635 CALL st_tpu(tensor3, zero, 3*3)
636 CALL st_tpu(vgrad, zero, 3*3)
637C
638 CALL inv_tpu(f0, tensor1, det)
639 CALL axb_tpu(f1, tensor1, rdfgrd, 3)
640C
641C------- APPROXIMATION TO VELOCITY GRADIENT: L=2/DT*(F-I)*(F+I)^{-1}-
642C
643 CALL at_tpu(one, rdfgrd, -one, pmident, tensor2, 3*3)
644 CALL at_tpu(one, rdfgrd, one, pmident, tensor1, 3*3)
645 CALL inv_tpu(tensor1, tensor3, det)
646 CALL axb_tpu(tensor2, tensor3, vgrad, 3)
647 CALL sst_tpu(two/dtime, vgrad, vgrad, 3*3)
648C
649C------- RATE OF DEFORMATION AND SPIN: D_UMAT=SYM(L), W_UMAT=SKW(L) -
650C
651 CALL st_tpu(d_umat, zero, 3*3)
652 CALL st_tpu(w_umat, zero, 3*3)
653C
654 CALL sym_tpu(vgrad, d_umat, 3)
655 CALL skw_tpu(vgrad, w_umat, 3)
656C
657C------- COMPUTE EDOT=SQRT(2/3*D:D) -
658C
659 CALL st_tpu(d_umat_vec, zero, 6)
660C
661 CALL mat3x3tovec6x1(d_umat,d_umat_vec,3)
662 CALL sp_tpu(d_umat_vec,d_umat_vec,coeff3,6)
663 coeff2 = sqrt((two/three)*coeff3)
664C
665C-------- ROTATION TENSORS FROM F=RU, I.E., R=F*U^(-1)
666C
667 CALL st_tpu(rot1, zero, 3*3)
668 CALL st_tpu(tensor1, zero, 3*3)
669C
670 CALL inv_tpu(u1, tensor1, det)
671 CALL axb_tpu(f1, tensor1, rot1, 3)
672C
673C------- TEMPERATURE DEPENDENCE OF PARAMETERS -
674 IF (coeff2 < small) THEN
675 coeff2 = em04
676 ENDIF
677C
678 CALL emod_tpu(deref,de0,theta,theta0,dve1,dve2,
679 & coeff2,dedot_ref,de)
680 dmu = de / (two+two*dpoiss)
681 dk = two*dmu*(one+dpoiss)/(three*(one-two*dpoiss))
682 dy0 = dc3*(theta-theta0)+dc4
683 ckappa1 = calphak1*dmu
684 destar_0 = dc5*(theta-theta0)+dc6
685 destar_s = dc7*(theta-theta0)+dc8
686 g0 = dc9*(theta-theta0)+dc10
687 ckappa2 = calphak2*dmu
688 des2_s = dc11*(theta-theta0)+dc12
689 dmu_r = dc1*(theta-theta0)+dc2
690 drs1 = dc13*(theta-theta0)+ dc14
691C
692 !-----------------------------------------------------
693 !CHECK PARAM
694 !-----------------------------------------------------
695 twomu0 = two*dmu
696 bulka0 = dk - two_third*dmu
697 !-------------------------------------------------------------72
698 !- AT THE BEGINNING OF THE STEP
699 !-------------------------------------------------------------72
700C
701 CALL st_tpu(fvm, zero, 3*3)
702 CALL st_tpu(fvm_inv, zero, 3*3)
703 CALL st_tpu(ftheta, zero, 3*3)
704 CALL st_tpu(ftheta_inv, zero, 3*3)
705 CALL st_tpu(tens0, zero, 3*3)
706 CALL st_tpu(tens1, zero, 3*3)
707 CALL st_tpu(fe0, zero, 3*3)
708 CALL st_tpu(ce0, zero, 3*3)
709 CALL st_tpu(ue0, zero, 3*3)
710 CALL st_tpu(ue0_inv, zero, 3*3)
711 CALL st_tpu(rote0, zero, 3*3)
712 CALL st_tpu(ee0, zero, 3*3)
713 CALL st_tpu(xm, zero, 3*3)
714C
715C DETERMINATION OF BP
716C
717 CALL vec6x1tomat3x3(fv,fvm,3)
718 CALL axbt_tpu(fvm, fvm, bp, 3)
719 CALL tr_tpu(bp, tr_bp, 3)
720C
721C DETERMINATION OF FE AND CE
722C
723 CALL sst_tpu(one+beta0*(theta-thetai), pmident, ftheta, 3*3)
724 CALL inv_tpu(ftheta, ftheta_inv, det)
725 CALL inv_tpu(fvm, fvm_inv, det)
726 CALL axb_tpu(f0, ftheta_inv, tens0, 3)
727 CALL axb_tpu(tens0, fvm_inv, fe0, 3)
728 CALL atxb_tpu(fe0, fe0, ce0, 3)
729C
730 DO i=1,ndi
731 DO j=1,ndi
732 a(i,j)= ce0(i,j)
733 ENDDO
734 ENDDO
735 CALL valprop_tpu(a, 3, 3, dce0, vce0, nrot, ier)
736 CALL eigsrt(dce0, vce0, 3, 3)
737 IF (ier==1) THEN
738 CALL err_tpu('ERROR IN EIG VAL')
739 ENDIF
740 !EIGENVALUES (AND ASSOC EIGENVECTORS) ARE ORDERED FROM LARGER
741 !TO SMALLER.
742 !REDEFINE AXIS(2) TO BE THE LARGEST IN ORDER TO IMPROVE ACCURACY
743 !IN THE CALCULATION OF THE ESHELBY TENSOR.
744 !IF DET(B)<0 MEANS THAT THE SYSTEM IS LEFT HANDED. IT IS MADE
745 !RIGHT HANDED BY EXCHANGING 1 AND 2.
746 sign=-one
747 CALL det_tpu(vce0, det_vce0)
748 IF(det_vce0 <= zero) sign=one
749 DO i=1,3
750 exchange=vce0(i,1)
751 vce0(i,1)=vce0(i,2)
752 vce0(i,2)=exchange*sign
753 ENDDO
754 exchange=dce0(1)
755 dce0(1)=dce0(2)
756 dce0(2)=exchange
757C
758 !COMPUTE UE0
759 d1=sqrt(dce0(1))
760 d2=sqrt(dce0(2))
761 d3=sqrt(dce0(3))
762 de1=log(sqrt(dce0(1)))
763 de2=log(sqrt(dce0(2)))
764 de3=log(sqrt(dce0(3)))
765 DO i=1,3
766 DO j=1,3
767 ue0(i,j)=d1*vce0(i,1)*vce0(j,1)
768 & +d2*vce0(i,2)*vce0(j,2)+d3*vce0(i,3)*vce0(j,3)
769 END DO
770 END DO
771 DO i=1,3
772 DO j=1,3
773 ee0(i,j)=de1*vce0(i,1)*vce0(j,1)
774 & +de2*vce0(i,2)*vce0(j,2)+de3*vce0(i,3)*vce0(j,3)
775 END DO
776 END DO
777 CALL inv_tpu(ue0, ue0_inv, det)
778 CALL axb_tpu(fe0, ue0_inv, rote0, 3)
779 CALL det_tpu(fe0, det_fe)
780C
781C==== DETERMINATION OF MANDEL STRESS M ===
782C
783 CALL tr_tpu(ee0, tr_ee0, 3)
784 coeff1_a = twomu0
785 coeff1_b = bulka0*tr_ee0
786C
787 CALL at_tpu(coeff1_a, ee0, coeff1_b, pmident, xm, 9)
788C
789 !-------------------------------------------------------------72
790 !- CHECK 2
791 !- ELASTIC STRESSES AT THE BEGINNING OF THE STEP
792 !-!------------------------------------------------------------72
793 !-------------------------------------------------------------72
794 !DETERMINATION OF FLOW RULE FOR TIME-DEPENDENT DASHPOT
795 !-------------------------------------------------------------72
796 !DETERMINATION OF EQSTRESS, PRESSURE, DIRECTION OF PLASTIC FLOW
797 CALL st_tpu(dev_alpham, zero, 3*3)
798 CALL st_tpu(betam0, zero, 3*3)
799 CALL st_tpu(dev_xm, zero, 3*3)
800 CALL st_tpu(alpham, zero, 3*3)
801 CALL st_tpu(xmeff, zero, 3*3)
802 CALL st_tpu(dev_meff, zero, 3*3)
803 CALL st_tpu(temp1, zero, 3*3)
804 CALL st_tpu(temp1_inv, zero, 3*3)
805 CALL st_tpu(temp2, zero, 3*3)
806 CALL st_tpu(temp2_inv, zero, 3*3)
807 CALL st_tpu(exp_dvdt, zero, 3*3)
808 CALL st_tpu(fvm_iter, zero, 3*3)
809 CALL st_tpu(dev_meff_vec, zero, 6)
810 DO i = 1, 6!NTENS
811 nv_vec(i)=pident(i)
812 ENDDO
813C
814 CALL vec6x1tomat3x3(beta, betam0, 3)
815 CALL sst_tpu(dmu_b, betam0, alpham, 3*3)
816 CALL dev_tpu(alpham, dev_alpham, 3)
817 CALL dev_tpu(xm, dev_xm, 3)
818 CALL at_tpu(one, dev_xm, -one, dev_alpham, dev_meff, 9)
819 CALL mat3x3tovec6x1(dev_meff,dev_meff_vec,3)
820 CALL sp_tpu(dev_meff_vec,dev_meff_vec, ps, 6)
821 norm_devmeff = sqrt(ps)
822C
823 dkappa1 = ckappa1*des10
824 dkappa2 = ckappa2*des20
825 dkappa = dkappa1 + dkappa2
826 CALL tr_tpu(xm, tr_m, 3)
827 pressure = -third*tr_m
828 eqstress = (norm_devmeff/sqrt(two))
829 & -(dy0+dkappa+alphap*pressure)
830C
831 !DETERMINATION OF GAMVDOT
832C
833 IF (eqstress < small) THEN
834 dgamv_iter = zero
835C WRITE(0 ,*) '!!!NO VISCOSITY!!!! ', DGAMV_ITER
836 ELSE
837C WRITE(0 ,*) '!!!VISCOSITY!!!! ', DGAMV_ITER
838 CALL sst_tpu(one/norm_devmeff, dev_meff_vec, nv_vec,6)
839C
840 gamv0 = gamv_ref*exp(-dq/(dr*theta))
841 dgamv_iter = gamv0*dtime*
842 & (sinh(eqstress*dv1/(two*dkb*theta)))**dm
843 IF (dgamv_iter > ep06) THEN
844 CALL err_tpu('PLASTICITY > 1.0D6!!!')
845 ENDIF
846C
847 ENDIF
848C
849C.......DETERMINATION OF INTERNAL STRAIN
850C
851 destar = destar0
852 & + (destar_s-g0*destar0)*dgamv_iter
853C
854 des1 = des10 + h0*
855 & (one-des10/destar)*dgamv_iter
856C
857 dlambdap_bar = one/ds3*sqrt(tr_bp)
858 des2 = des20 + h1*(dlambdap_bar-one)*
859 & (one-des20/des2_s)*dgamv_iter
860C
861C.......DETERMINATION OF BETA
862C
863 CALL st_tpu(dv_vec, zero, 6)
864 CALL st_tpu(dv, zero, 3*3)
865 CALL st_tpu(tensor1, zero, 3*3)
866 CALL st_tpu(tensor2, zero, 3*3)
867 CALL st_tpu(tensor3, zero, 3*3)
868 CALL st_tpu(tensor4, zero, 3*3)
869 CALL st_tpu(tensor5, zero, 3*3)
870 CALL st_tpu(betam, zero, 3*3)
871C
872 CALL sst_tpu(dgamv_iter/(dtime*sqrt(two)), nv_vec, dv_vec, 6)
873 CALL vec6x1tomat3x3(dv_vec, dv,3)
874C
875C COMPUTE BETA_DOT
876C
877 CALL axb_tpu(dv, betam0, tensor1, 3)
878 CALL axb_tpu(betam0, dv, tensor2, 3)
879 CALL at_tpu(one, tensor1, one, tensor2, tensor3, 9)
880C
881 CALL sst_tpu(drs1,
882 & tensor3, tensor4, 3*3)
883 CALL at_tpu(one, betam0, dtime, tensor4, betam, 9)
884C
885C EIGEN VALUES AND VECTOR OF BETA
886C
887 DO i=1,ndi
888 DO j=1,ndi
889 a(i,j)= betam(i,j)
890 ENDDO
891 ENDDO
892 CALL valprop_tpu(a, 3, 3, dce1, vce1, nrot, ier)
893 CALL eigsrt(dce1, vce1, 3, 3)
894 IF (ier==1) THEN
895 CALL err_tpu('ERROR IN EIG VAL')
896 ENDIF
897 sign=-one
898 CALL det_tpu(vce1, det_vce1)
899 IF(det_vce1<=0.) sign=one
900 DO i=1,3
901 exchange=vce1(i,1)
902 vce1(i,1)=vce1(i,2)
903 vce1(i,2)=exchange*sign
904 ENDDO
905 exchange=dce1(1)
906 dce1(1)=dce1(2)
907 dce1(2)=exchange
908C
909 DO i=1,3
910 DO j=1,3
911 betam0(i,j)=dce1(1)*vce1(i,1)*vce1(j,1)
912 & +dce1(2)*vce1(i,2)*vce1(j,2)+dce1(3)*vce1(i,3)*vce1(j,3)
913 END DO
914 END DO
915C
916C...... DETERMINATION OF MATERIAL FUNCTION MU_B
917 dmu_b = dmu_r*(one-(dce1(1)+dce1(2)+dce1(3)-three)/dlambda_l)**(-1.0)
918C...... DETERMINATION OF ALPHA
919C
920C
921 CALL at_tpu(dtime, dv, one, pmident, temp1, 9)
922 CALL axb_tpu(temp1, fvm, fvm_iter, 3)
923 CALL mat3x3tovec6x1(fvm_iter,fv_iter,3)
924C--------------------------------------------------------------------72
925C-------- ELASTIC STRESSES AT THE END OF THE STEP
926C--------------------------------------------------------------------72
927C
928C---- COMPUTE: FE->CE->UE->RE->EE
929 CALL st_tpu(fvm_inv, zero, 3*3)
930 CALL st_tpu(ftheta, zero, 3*3)
931 CALL st_tpu(ftheta_inv, zero, 3*3)
932 CALL st_tpu(tens0, zero, 3*3)
933 CALL st_tpu(tens1, zero, 3*3)
934 CALL st_tpu(fe1, zero, 3*3)
935 CALL st_tpu(ce1, zero, 3*3)
936 CALL st_tpu(ue1, zero, 3*3)
937 CALL st_tpu(ue1_inv, zero, 3*3)
938 CALL st_tpu(rote1, zero, 3*3)
939 CALL st_tpu(ee1, zero, 3*3)
940C
941C DETERMINATION OF FE AND CE
942C
943 CALL sst_tpu(one+beta0*(theta-thetai), pmident, ftheta, 3*3)
944 CALL inv_tpu(ftheta, ftheta_inv, det)
945 CALL inv_tpu(fvm_iter, fvm_inv, det)
946 CALL axb_tpu(f1, ftheta_inv, tens0, 3)
947 CALL axb_tpu(tens0, fvm_inv, fe1, 3)
948 CALL atxb_tpu(fe1, fe1, ce1, 3)
949 DO i=1,ndi
950 DO j=1,ndi
951 a(i,j)= ce1(i,j)
952 ENDDO
953 ENDDO
954 CALL valprop_tpu(a, 3, 3, dce1, vce1, nrot, ier)
955 CALL eigsrt(dce1, vce1, 3, 3)
956 IF (ier==1) THEN
957 CALL err_tpu('ERROR IN EIG VAL')
958 ENDIF
959 !*** EIGENVALUES (AND ASSOC EIGENVECTORS) ARE ORDERED FROM LARGER TO SMALLER. !
960 !*** REDEFINE AXIS(2) TO BE THE LARGEST IN ORDER TO IMPROVE ACCURACY IN THE !
961 ! CALCULATION OF THE ESHELBY TENSOR. !
962 !*** IF DET(B)<0 MEANS THAT THE SYSTEM IS LEFT HANDED. IT IS MADE RIGHT !
963 ! HANDED BY EXCHANGING 1 AND 2. !
964 sign=-one
965 CALL det_tpu(vce1, det_vce1)
966 IF(det_vce1<=zero) sign=one
967 DO i=1,3
968 exchange=vce1(i,1)
969 vce1(i,1)=vce1(i,2)
970 vce1(i,2)=exchange*sign
971 ENDDO
972 exchange=dce1(1)
973 dce1(1)=dce1(2)
974 dce1(2)=exchange
975C
976C------ COMPUTE UE1
977C
978 d1=sqrt(dce1(1))
979 d2=sqrt(dce1(2))
980 d3=sqrt(dce1(3))
981 de1=log(sqrt(dce1(1)))
982 de2=log(sqrt(dce1(2)))
983 de3=log(sqrt(dce1(3)))
984C
985 DO i=1,3
986 DO j=1,3
987 ue1(i,j)=d1*vce1(i,1)*vce1(j,1)
988 & +d2*vce1(i,2)*vce1(j,2)+d3*vce1(i,3)*vce1(j,3)
989 END DO
990 END DO
991C
992 DO i=1,3
993 DO j=1,3
994 ee1(i,j)=de1*vce1(i,1)*vce1(j,1)
995 & +de2*vce1(i,2)*vce1(j,2)+de3*vce1(i,3)*vce1(j,3)
996 END DO
997 END DO
998C
999 CALL inv_tpu(ue1, ue1_inv, det)
1000 CALL axb_tpu(fe1, ue1_inv, rote1, 3)
1001 CALL det_tpu(fe1, det_fe)
1002C
1003C=======DETERMINATION OF TAU
1004C
1005 CALL st_tpu(xm, zero, 3*3)
1006 CALL st_tpu(taum, zero, 3*3)
1007C
1008 CALL tr_tpu(ee1, tr_ee1, 3)
1009 coeff1_a = twomu0
1010 coeff1_b = bulka0*tr_ee1
1011 CALL at_tpu(coeff1_a,ee1, coeff1_b, pmident, xm, 9)
1012 CALL qaqt_tpu(rote1, xm, taum, 3)
1013C
1014C-------DETERMINATION OF SIGMA AND INTERNAL STATE VARIABLES
1015C
1016 CALL st_tpu(sigmam, zero, 3*3)
1017 CALL st_tpu(sigma_q, zero, 3*3)
1018C
1019 CALL sst_tpu(one/det_fe, taum, sigmam, 3*3)
1020C
1021C-------DETERMINATION OF SIGMA (IN THE ROTATED CONFIGURATION)
1022C
1023 CALL qtaq_tpu(rot1, sigmam, sigma_q, 3)
1024C
1025C====================================================================72
1026C============ UPDATE VARIABLES
1027C====================================================================72
1028C
1029 CALL sp_tpu(dv_vec, dv_vec,ps,6)
1030 gamv_iter = gamv + dgamv_iter
1031 gamveq = gamveq +sqrt((two/three)*ps)*dtime
1032C
1033C-------TEMPERATURE, ASSUMING ADIABATIC CONDITIONS ----
1034C
1035 IF (abs(d_theta_opt-zero) < small .OR.
1036 & abs(d_theta_opt-one) < small) factor = zero
1037C
1038 CALL st_tpu(xm_vec, zero, 6)
1039C
1040 rho_p = rho_p*(1.42d0*thetag+44.7d0)
1041 & / (1.42d0*thetag+0.15d0*theta)
1042 cv = cv*(0.106d0+three*em03*theta)
1043C
1044 CALL mat3x3tovec6x1(xm, xm_vec, 3)
1045 CALL sp_tpu(xm_vec,dv_vec,ps,6)
1046 dtheta_iter = (dtime*factor*ps) /(cv*rho_p)
1047 thetaiter = theta + dtheta_iter
1048 thetadotiter = dtheta_iter/dtime
1049C=======DETERMINATION OF SIGMA AND INTERNAL STATE VARIABLES
1050C
1051C-------SIGMA_Q(I,J)
1052C
1053 signxx(ie) = sigma_q(1,1)
1054 signyy(ie) = sigma_q(2,2)
1055 signzz(ie) = sigma_q(3,3)
1056 signxy(ie) = sigma_q(1,2)
1057 signyz(ie) = sigma_q(2,3)
1058 signzx(ie) = sigma_q(1,3)
1059
1060C
1061 uvar(ie, 1) = fv_iter(1)
1062 uvar(ie, 2) = fv_iter(2)
1063 uvar(ie, 3) = fv_iter(3)
1064 uvar(ie, 4) = fv_iter(4)
1065 uvar(ie, 5) = fv_iter(5)
1066 uvar(ie, 6) = fv_iter(6)
1067 uvar(ie, 7) = des1
1068 uvar(ie, 8) = destar
1069 uvar(ie, 9) = des2
1070 uvar(ie, 10) = betam0(1,1)
1071 uvar(ie, 11) = betam0(2,2)
1072 uvar(ie, 12) = betam0(3,3)
1073 uvar(ie, 13) = betam0(1,2)
1074 uvar(ie, 14) = betam0(2,3)
1075 uvar(ie, 15) = betam0(3,1)
1076 uvar(ie, 16) = dmu_b
1077 uvar(ie, 17) = gamv_iter
1078 uvar(ie, 18) = dgamv_iter
1079 uvar(ie, 19) = gamveq
1080 uvar(ie, 20) = thetaiter
1081 uvar(ie, 21) = thetadotiter
1082 uvar(ie,22) = f1(1,1) !F011
1083 uvar(ie,23) = f1(2,2) !F022
1084 uvar(ie,24) = f1(3,3) !F033
1085 uvar(ie,25) = f1(1,2) !F012
1086 uvar(ie,26) = f1(2,1) !F021
1087 uvar(ie,27) = f1(1,3) !F013
1088 uvar(ie,28) = f1(3,1) !F031
1089 uvar(ie,29) = f1(2,3) !F023
1090 uvar(ie,30) = f1(3,2) !F032
1091 uvar(ie,31) = u1(1,1) !U011
1092 uvar(ie,32) = u1(2,2) !U022
1093 uvar(ie,33) = u1(3,3) !U033
1094 uvar(ie,34) = u1(1,2) !U012
1095 uvar(ie,35) = u1(2,1) !U021
1096 uvar(ie,36) = u1(1,3) !U013
1097 uvar(ie,37) = u1(3,1) !U031
1098 uvar(ie,38) = u1(2,3) !U023
1099 uvar(ie,39) = u1(3,2) !U032
1100 uvar(ie,40) = log(u1(1,1)) !E111
1101 uvar(ie,41) = log(u1(2,2)) !E122
1102 uvar(ie,42) = log(u1(3,3)) !E133
1103
1104C
1105
1106C sound velocity
1107 soundsp(ie) = sqrt((dk+4*dmu/3)/rho0(ie))
1108 viscmax(ie) = zero
1109
1110C................................... MODEL IMPLEMENTATION CHOICE
1111C
111220 ENDDO
1113C
1114C................................... END ELEMENT LOOP
1115
1116
1117C Temperature
1118 temp(1:nel) = tempel(1:nel)
1119
1120
1121 RETURN
1122 END
1123C
1124
1125
1126
1127C=====================================================================72
1128C=====================================================================72
1129C---- SUPPORTING ROUTINES FOR TPU MATERIAL ROUTINE
1130C=====================================================================72
1131C=====================================================================72
1132C
1133!||====================================================================
1134!|| st_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1135!||--- called by ------------------------------------------------------
1136!|| atxb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1137!|| axb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1138!|| axbt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1139!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1140!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1141!|| skw_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1142!||====================================================================
1143 SUBROUTINE st_tpu(A, VALUE, N)
1144#include "implicit_f.inc"
1145
1146 INTEGER N, I
1147 my_real VALUE
1148 my_real A(N)
1149C
1150C---------------------------------------------------------------------72
1151C
1152 DO I = 1, n
1153 a(i) = VALUE
1154 ENDDO
1155
1156 RETURN
1157 END
1158C
1159C=====================================================================72
1160C
1161C
1162C=====================================================================72
1163C
1164!||====================================================================
1165!|| axb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1166!||--- called by ------------------------------------------------------
1167!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1168!|| qaqt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1169!|| qtaq_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1170!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1171!||--- calls -----------------------------------------------------
1172!|| st_tpu ../engine/source/materials/mat/mat101/sigeps101.f
1173!||====================================================================
1174 SUBROUTINE axb_tpu(A, B, C, N)
1175#include "implicit_f.inc"
1176
1177 INTEGER N, I, J, K
1178 my_real A(N, N), B(N, N), C(N, N)
1179C
1180C---------------------------------------------------------------------72
1181C
1182 CALL ST_TPU(C,ZERO, N*N)
1183C
1184 DO I = 1, n
1185 DO j = 1, n
1186 DO k = 1, n
1187 c(i,j) = c(i,j) + a(i,k) * b(k,j)
1188 ENDDO
1189 ENDDO
1190 ENDDO
1191C
1192 RETURN
1193 END
1194C
1195C=====================================================================72
1196C
1197C
1198C=====================================================================72
1199C
1200!||====================================================================
1201!|| atxb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1202!||--- called by ------------------------------------------------------
1203!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1204!|| qtaq_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1205!|| qtaqmt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1206!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1207!||--- calls -----------------------------------------------------
1208!|| st_tpu ../engine/source/materials/mat/mat101/sigeps101.f
1209!||====================================================================
1210 SUBROUTINE atxb_tpu(A, B, C, N)
1211#include "implicit_f.inc"
1212
1213 INTEGER N, I, J, K
1214 my_real A(N, N), B(N, N), C(N, N)
1215C
1216C---------------------------------------------------------------------72
1217C
1218 CALL ST_TPU(C,ZERO, N*N)
1219C
1220 DO I = 1, n
1221 DO j = 1, n
1222 DO k = 1, n
1223 c(i,j) = c(i,j) + a(k,i) * b(k,j)
1224 ENDDO
1225 ENDDO
1226 ENDDO
1227C
1228 RETURN
1229 END
1230C=====================================================================72
1231C
1232C
1233C=====================================================================72
1234C
1235!||====================================================================
1236!|| sst_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1237!||--- called by ------------------------------------------------------
1238!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1239!||====================================================================
1240 SUBROUTINE sst_tpu(COEFA, A, B, N)
1241#include "implicit_f.inc"
1242
1243 INTEGER N, I
1244 my_real COEFA
1245 my_real A(N), B(N)
1246C
1247C---------------------------------------------------------------------72
1248C
1249 DO I = 1, n
1250 b(i) = coefa*a(i)
1251 ENDDO
1252C
1253 RETURN
1254 END
1255C
1256C=====================================================================72
1257C
1258C
1259C=====================================================================72
1260C
1261!||====================================================================
1262!|| axbt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1263!||--- called by ------------------------------------------------------
1264!|| qaqt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1265!|| qtaqmt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1266!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1267!||--- calls -----------------------------------------------------
1268!|| st_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1269!||====================================================================
1270 SUBROUTINE axbt_tpu(A, B, C, N)
1271#include "implicit_f.inc"
1272
1273 INTEGER N, I, J, K
1274 my_real A(N, N), B(N, N), C(N, N)
1275C
1276C---------------------------------------------------------------------72
1277C
1278 CALL ST_TPU(C,ZERO, N*N)
1279C
1280 DO I = 1, n
1281 DO j = 1, n
1282 DO k = 1, n
1283 c(i,j) = c(i,j) + a(i,k) * b(j,k)
1284 ENDDO
1285 ENDDO
1286 ENDDO
1287C
1288 RETURN
1289 END
1290C
1291C=====================================================================72
1292C
1293C
1294C=====================================================================72
1295C
1296!||====================================================================
1297!|| tr_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1298!||--- called by ------------------------------------------------------
1299!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1300!|| dev_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1301!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1302!||--- calls -----------------------------------------------------
1303!|| err_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1304!||====================================================================
1305 SUBROUTINE tr_tpu(TENSOR, TR_TENSOR, N)
1306#include "implicit_f.inc"
1307
1308 INTEGER N
1309 my_real TENSOR(N*N)
1310 my_real TR_TENSOR
1311C
1312C---------------------------------------------------------------------72
1313C
1314 IF (N /= 3)
1315 & call err_tpu('TRACEOFTENSOR: N =! 3')
1316C
1317 tr_tensor = tensor(1) + tensor(5) + tensor(9)
1318C
1319 RETURN
1320 END
1321C
1322C=====================================================================72
1323C
1324C
1325C=====================================================================72
1326C
1327!||====================================================================
1328!|| det_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1329!||--- called by ------------------------------------------------------
1330!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1331!|| inv_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1332!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.f
1333!||====================================================================
1334 SUBROUTINE det_tpu(TENSOR, TENSOR_DET)
1335#include "implicit_f.inc"
1336
1337 my_real tensor(3, 3)
1338 my_real det11, det12, det13, det21, det22, det23
1339 my_real tensor_det
1340C
1341C---------------------------------------------------------------------72
1342C
1343C DETERMINANT OF A SECOND ORDER TENSOR (3X3 MATRIX)
1344C
1345 det11 = tensor(1, 1) * tensor(2, 2) * tensor(3, 3)
1346 det12 = tensor(1, 2) * tensor(2, 3) * tensor(3, 1)
1347 det13 = tensor(1, 3) * tensor(2, 1) * tensor(3, 2)
1348 det21 = tensor(1, 1) * tensor(2, 3) * tensor(3, 2)
1349 det22 = tensor(1, 2) * tensor(2, 1) * tensor(3, 3)
1350 det23 = tensor(1, 3) * tensor(2, 2) * tensor(3, 1)
1351C
1352 tensor_det = det11 + det12 + det13 - det21 - det22 - det23
1353C
1354 RETURN
1355 END
1356C
1357C=====================================================================72
1358C
1359C
1360C=====================================================================72
1361C
1362!||====================================================================
1363!|| inv_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1364!||--- called by ------------------------------------------------------
1365!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1366!|| qtaqmt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1367!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1368!||--- calls -----------------------------------------------------
1369!|| det_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1370!||====================================================================
1371 SUBROUTINE inv_tpu(TENSOR, TENSOR_INV, TENSOR_DET)
1372#include "implicit_f.inc"
1373
1374 my_real tensor_det
1375 my_real tensor(3, 3), tensor_inv(3, 3)
1376 my_real tinv11, tinv12, tinv13
1377 my_real tinv21, tinv22, tinv23
1378 my_real tinv31, tinv32, tinv33
1379C
1380C---------------------------------------------------------------------72
1381C
1382 CALL det_tpu(tensor, tensor_det)
1383C
1384 tinv11 = tensor(2, 2) * tensor(3, 3) - tensor(2, 3) * tensor(3, 2)
1385 tinv12 = tensor(1, 3) * tensor(3, 2) - tensor(1, 2) * tensor(3, 3)
1386 tinv13 = tensor(1, 2) * tensor(2, 3) - tensor(1, 3) * tensor(2, 2)
1387 tinv21 = tensor(2, 3) * tensor(3, 1) - tensor(2, 1) * tensor(3, 3)
1388 tinv22 = tensor(1, 1) * tensor(3, 3) - tensor(1, 3) * tensor(3, 1)
1389 tinv23 = tensor(1, 3) * tensor(2, 1) - tensor(1, 1) * tensor(2, 3)
1390 tinv31 = tensor(2, 1) * tensor(3, 2) - tensor(2, 2) * tensor(3, 1)
1391 tinv32 = tensor(3, 1) * tensor(1, 2) - tensor(1, 1) * tensor(3, 2)
1392 tinv33 = tensor(1, 1) * tensor(2, 2) - tensor(1, 2) * tensor(2, 1)
1393C
1394 IF (tensor_det /= zero) THEN
1395 tensor_inv(1, 1) = tinv11 / tensor_det
1396 tensor_inv(1, 2) = tinv12 / tensor_det
1397 tensor_inv(1, 3) = tinv13 / tensor_det
1398 tensor_inv(2, 1) = tinv21 / tensor_det
1399 tensor_inv(2, 2) = tinv22 / tensor_det
1400 tensor_inv(2, 3) = tinv23 / tensor_det
1401 tensor_inv(3, 1) = tinv31 / tensor_det
1402 tensor_inv(3, 2) = tinv32 / tensor_det
1403 tensor_inv(3, 3) = tinv33 / tensor_det
1404 ENDIF
1405C
1406 RETURN
1407 END
1408C
1409C=====================================================================72
1410C
1411C
1412C=====================================================================72
1413C
1414!||====================================================================
1415!|| at_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1416!||--- called by ------------------------------------------------------
1417!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1418!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1419!||====================================================================
1420 SUBROUTINE at_tpu(COEFA, A, COEFB, B, C, N)
1421#include "implicit_f.inc"
1422
1423 INTEGER N, I
1424 my_real COEFA, COEFB
1425 my_real A(N), B(N), C(N)
1426C
1427C---------------------------------------------------------------------72
1428C
1429 DO I = 1, n
1430 c(i) = coefa*a(i) + coefb*b(i)
1431C WRITE(0 , *) 'CHECK AT_TPU '
1432C WRITE(0 , *) 'COEFA ',COEFA,' ','COEFB ',COEFB
1433C WRITE(0 , *) 'A ',I,' ',A(I)
1434C WRITE(0 , *) 'B ',I,' ',B(I)
1435C WRITE(0 , *) 'C ',I,' ',C(I)
1436 ENDDO
1437C
1438 RETURN
1439 END
1440C
1441C=====================================================================72
1442C
1443C=====================================================================72
1444C
1445!||====================================================================
1446!|| dev_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1447!||--- called by ------------------------------------------------------
1448!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1449!||--- calls -----------------------------------------------------
1450!|| err_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1451!|| tr_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1452!||====================================================================
1453 SUBROUTINE dev_tpu(TENSOR, TENSORDEV, N)
1454#include "implicit_f.inc"
1455
1456 INTEGER N
1457 my_real TENSKK, TR_TENSOR
1458 my_real TENSOR(N*N), TENSORDEV(N*N)
1459C
1460C---------------------------------------------------------------------72
1461C
1462 IF (N /= 3)
1463 & call err_tpu('DEVIATORICTENSOR: N =! 3')
1464C
1465 CALL tr_tpu(tensor, tr_tensor, n)
1466 tenskk = tr_tensor / three
1467C
1468 tensordev(1) = tensor(1) - tenskk
1469 tensordev(5) = tensor(5) - tenskk
1470 tensordev(9) = tensor(9) - tenskk
1471C
1472 tensordev(2) = tensor(2)
1473 tensordev(3) = tensor(3)
1474 tensordev(6) = tensor(6)
1475C
1476 tensordev(4) = tensor(4)
1477 tensordev(7) = tensor(7)
1478 tensordev(8) = tensor(8)
1479C
1480 RETURN
1481 END
1482C
1483C=====================================================================72
1484C
1485C
1486C=====================================================================72
1487C
1488!||====================================================================
1489!|| qaqt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1490!||--- called by ------------------------------------------------------
1491!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1492!||--- calls -----------------------------------------------------
1493!|| axb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1494!|| axbt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1495!||====================================================================
1496 SUBROUTINE qaqt_tpu(Q, A, B, N)
1497#include "implicit_f.inc"
1498
1499 INTEGER N
1500 my_real Q(N,N), A(N,N), B(N,N)
1501 my_real TMP(N,N)
1502C
1503C---------------------------------------------------------------------72
1504C
1505 CALL AXB_TPU(Q, A, TMP, N) ! TMP = Q*A
1506 CALL AXBT_TPU(TMP, Q, B, N) ! B = TMP*QT
1507C
1508 RETURN
1509 END
1510C
1511C=====================================================================72
1512C
1513C
1514C=====================================================================72
1515C
1516!||====================================================================
1517!|| qtaq_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1518!||--- called by ------------------------------------------------------
1519!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1520!||--- calls -----------------------------------------------------
1521!|| atxb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1522!|| axb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1523!||====================================================================
1524 SUBROUTINE qtaq_tpu(Q, A, B, N)
1525#include "implicit_f.inc"
1526
1527 INTEGER N
1528 my_real Q(N,N), A(N,N), B(N,N)
1529 my_real TMP(N,N)
1530C
1531C---------------------------------------------------------------------72
1532C
1533 CALL ATXB_TPU(Q, A, TMP, N) ! TMP = QT*A
1534 CALL AXB_TPU(TMP, Q, B, N) ! B = TMP*Q
1535C
1536 RETURN
1537 END
1538C
1539C=====================================================================72
1540C
1541C
1542C=====================================================================72
1543C
1544!||====================================================================
1545!|| sym_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1546!||--- called by ------------------------------------------------------
1547!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1548!||--- calls -----------------------------------------------------
1549!|| err_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1550!||====================================================================
1551 SUBROUTINE sym_tpu(TENSOR, SYMMTENSOR, N)
1552#include "implicit_f.inc"
1553
1554 INTEGER N
1555 my_real TENSOR(N*N), SYMMTENSOR(N*N)
1556C
1557C---------------------------------------------------------------------72
1558C
1559 IF ((N /= 2) .AND. (n /= 3))
1560 & CALL err_tpu('SYM_TPU: N =! 2 OR 3')
1561C
1562 IF (n == 2) THEN
1563 symmtensor(1) = tensor(1)
1564 symmtensor(4) = tensor(4)
1565 symmtensor(3) = half*(tensor(3) + tensor(2))
1566 symmtensor(2) = symmtensor(3)
1567 ELSE ! N = 3
1568 symmtensor(1) = tensor(1) ! 11
1569 symmtensor(5) = tensor(5) ! 22
1570 symmtensor(9) = tensor(9) ! 33
1571 symmtensor(4) = half*(tensor(4) + tensor(2)) ! 12
1572 symmtensor(7) = half*(tensor(7) + tensor(3)) ! 13
1573 symmtensor(8) = half*(tensor(8) + tensor(6)) ! 23
1574 symmtensor(2) = symmtensor(4) ! 21
1575 symmtensor(3) = symmtensor(7) ! 31
1576 symmtensor(6) = symmtensor(8) ! 32
1577 ENDIF
1578C
1579 RETURN
1580 END
1581C
1582C=====================================================================72
1583C
1584C
1585C=====================================================================72
1586C
1587!||====================================================================
1588!|| skw_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1589!||--- called by ------------------------------------------------------
1590!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.f
1591!||--- calls -----------------------------------------------------
1592!|| err_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1593!|| st_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1594!||====================================================================
1595 SUBROUTINE skw_tpu(TENSOR, SKEWTENSOR, N)
1596#include "implicit_f.inc"
1597
1598 INTEGER N
1599 my_real TENSOR(N*N), SKEWTENSOR(N*N)
1600C
1601C---------------------------------------------------------------------72
1602C
1603 IF ((N /= 2) .AND. (n /= 3))
1604 & CALL err_tpu('SKW_TPU: N =! 2 OR 3')
1605C
1606 CALL st_tpu(skewtensor,zero, n*n)
1607C
1608 IF (n == 2) THEN
1609 skewtensor(3) = half*(tensor(3) - tensor(2))
1610 skewtensor(2) = - skewtensor(3)
1611 ELSE ! N = 3
1612 skewtensor(4) = half*(tensor(4) - tensor(2))
1613 skewtensor(7) = half*(tensor(7) - tensor(3))
1614 skewtensor(8) = half*(tensor(8) - tensor(6))
1615 skewtensor(2) = - skewtensor(4)
1616 skewtensor(3) = - skewtensor(7)
1617 skewtensor(6) = - skewtensor(8)
1618 ENDIF
1619C
1620 RETURN
1621 END
1622C
1623C=====================================================================72
1624C
1625C
1626C=====================================================================72
1627C
1628!||====================================================================
1629!|| mat3x3tovec6x1 ../engine/source/materials/mat/mat101/sigeps101.F
1630!||--- called by ------------------------------------------------------
1631!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1632!||====================================================================
1633 SUBROUTINE mat3x3tovec6x1(
1634 & MATRIX, VECTOR, N
1635 & )
1636C
1637#include "implicit_f.inc"
1638
1639C
1640 INTEGER N
1641 my_real MATRIX(N*N), VECTOR(6)
1642C
1643C---------------------------------------------------------------------72
1644C
1645 VECTOR(1) = matrix(1)
1646 vector(2) = matrix(5)
1647 vector(3) = matrix(9)
1648 vector(4) = matrix(2)
1649 vector(5) = matrix(3)
1650 vector(6) = matrix(6)
1651C
1652 RETURN
1653 END
1654C
1655C=====================================================================72
1656C
1657C
1658C=====================================================================72
1659C
1660!||====================================================================
1661!|| vec6x1tomat3x3 ../engine/source/materials/mat/mat101/sigeps101.F
1662!||--- called by ------------------------------------------------------
1663!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1664!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1665!||====================================================================
1666 SUBROUTINE vec6x1tomat3x3(
1667 & VECTOR, MATRIX, N
1668 & )
1669C
1670#include "implicit_f.inc"
1671
1672C
1673 INTEGER N
1674 my_real MATRIX(N*N), VECTOR(6)
1675C
1676C---------------------------------------------------------------------72
1677C
1678 MATRIX(1) = vector(1)
1679 matrix(5) = vector(2)
1680 matrix(9) = vector(3)
1681 matrix(2) = vector(4)
1682 matrix(3) = vector(5)
1683 matrix(6) = vector(6)
1684 matrix(4) = vector(4)
1685 matrix(7) = vector(5)
1686 matrix(8) = vector(6)
1687C
1688 RETURN
1689 END
1690C
1691C=====================================================================72
1692C
1693C
1694C=====================================================================72
1695C
1696 ! my_real FUNCTION SP_TPU(A, B, N)
1697 ! IMPLICIT NONE
1698 ! INTEGER N, I
1699 ! my_real A(N), B(N)
1700! C
1701! C---------------------------------------------------------------------72
1702! C
1703 ! SP_TPU =ZERO
1704 ! DO I = 1, 3
1705 ! SP_TPU = SP_TPU + A(I)*B(I)
1706 ! ENDDO
1707
1708 ! DO I = 4, N
1709 ! SP_TPU = SP_TPU + TWO*A(I)*B(I)
1710 ! ENDDO
1711
1712 ! RETURN
1713 ! END
1714
1715
1716!||====================================================================
1717!|| sp_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1718!||--- called by ------------------------------------------------------
1719!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1720!||====================================================================
1721 SUBROUTINE sp_tpu(A, B, c, N)
1722#include "implicit_f.inc"
1723
1724 INTEGER N, I
1725 my_real A(N), B(N), c
1726C
1727C---------------------------------------------------------------------72
1728C
1729 c =zero
1730 DO i = 1, 3
1731 c = c + a(i)*b(i)
1732 ENDDO
1733
1734 DO i = 4, n
1735 c = c + two*a(i)*b(i)
1736 ENDDO
1737
1738 RETURN
1739 END
1740
1741
1742C
1743C=====================================================================72
1744C
1745C
1746C=====================================================================72
1747C
1748!||====================================================================
1749!|| qtaqmt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1750!||--- calls -----------------------------------------------------
1751!|| atxb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1752!|| axbt_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1753!|| inv_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1754!||====================================================================
1755 SUBROUTINE qtaqmt_tpu(Q, A, B, N)
1756#include "implicit_f.inc"
1757
1758 INTEGER N
1759 my_real Q(N,N), A(N,N), B(N,N)
1760 my_real TMP(N,N), QM1(N,N)
1761 my_real DET
1762C
1763C---------------------------------------------------------------------72
1764C
1765 CALL ATXB_TPU(Q, A, TMP, N) ! TMP = QT*A
1766 CALL INV_TPU(Q, QM1, DET)
1767 CALL AXBT_TPU(TMP, QM1, B, N) ! B = TMP*QMT
1768C
1769 RETURN
1770 END
1771C
1772C=====================================================================72
1773C
1774C
1775C=====================================================================72
1776C
1777C
1778!||====================================================================
1779!|| emod_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1780!||--- called by ------------------------------------------------------
1781!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1782!||====================================================================
1783 SUBROUTINE emod_tpu(AREF,A0 ,C ,C0 ,DE1 ,DE2, D1,DREF,RES)
1784! CALL EMOD_TPU (DEREF,DE0,THETA,THETA0,DVE1,DVE2,
1785! & EM04,DEDOT_REF,DE)
1786
1787#include "implicit_f.inc"
1788
1789 my_real aref,a0,c,c0,de1,de2,d1,dref,res
1790C---------------------------------------------------------------------72
1791
1792C
1793 res = (aref + a0*(c-c0))*(one+de1/
1794 & (one+ exp(-(log10(max(em20,d1))-log10(max(em20, dref)) ) /max(em20,de2))))
1795C
1796 RETURN
1797 END
1798C
1799C=====================================================================72
1800C
1801C
1802C=====================================================================72
1803C
1804C
1805!||====================================================================
1806!|| valprop_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1807!||--- called by ------------------------------------------------------
1808!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1809!||====================================================================
1810 SUBROUTINE valprop_tpu(A,N,NP,D,V,NROT,IER)
1811#include "implicit_f.inc"
1812
1813 INTEGER N,NP,NROT,NMAX
1814 my_real a(np,np),d(np),v(np,np)
1815 parameter(nmax=500)
1816 INTEGER I,IP,IQ,J,IER
1817 my_real C,G,H,S,SM,T,TAU,THETA,TRESH,B(NMAX),Z(NMAX)
1818 DO 12 IP=1,n
1819 DO 11 iq=1,n
1820 v(ip,iq)=0.
182111 CONTINUE
1822 v(ip,ip)=one
182312 CONTINUE
1824 DO 13 ip=1,n
1825 b(ip)=a(ip,ip)
1826 d(ip)=b(ip)
1827 z(ip)=zero
182813 CONTINUE
1829 nrot=0
1830 DO 24 i=1,50
1831 sm=zero
1832 DO 15 ip=1,n-1
1833 DO 14 iq=ip+1,n
1834 sm=sm+abs(a(ip,iq))
1835
183614 CONTINUE
183715 CONTINUE
1838C
1839 IF(sm==0.)THEN
1840 ier=0
1841 RETURN
1842 ENDIF
1843C
1844 IF(i<4)THEN
1845 tresh=0.2d0*sm/n**2
1846 ELSE
1847 tresh=zero
1848 ENDIF
1849 DO 22 ip=1,n-1
1850 DO 21 iq=ip+1,n
1851 g=hundred*abs(a(ip,iq))
1852 IF((i>4).AND.(abs(d(ip))+
1853 *g==abs(d(ip))).AND.(abs(d(iq))+g==abs(d(iq))))THEN
1854 a(ip,iq)=0.
1855 ELSE IF(abs(a(ip,iq))>tresh)THEN
1856 h=d(iq)-d(ip)
1857 IF(abs(h)+g==abs(h))THEN
1858 t=a(ip,iq)/h
1859C
1860 ELSE
1861 theta=half*h/a(ip,iq)
1862 t=one/(abs(theta)+sqrt(one+theta**2))
1863 IF(theta<zero)t=-t
1864 ENDIF
1865 c=one/sqrt(one+t**2)
1866 s=t*c
1867 tau=s/(one+c)
1868 h=t*a(ip,iq)
1869 z(ip)=z(ip)-h
1870 z(iq)=z(iq)+h
1871 d(ip)=d(ip)-h
1872 d(iq)=d(iq)+h
1873 a(ip,iq)=zero
1874 DO 16 j=1,ip-1
1875 g=a(j,ip)
1876 h=a(j,iq)
1877C
1878 a(j,ip)=g-s*(h+g*tau)
1879 a(j,iq)=h+s*(g-h*tau)
188016 CONTINUE
1881 DO 17 j=ip+1,iq-1
1882 g=a(ip,j)
1883 h=a(j,iq)
1884 a(ip,j)=g-s*(h+g*tau)
1885 a(j,iq)=h+s*(g-h*tau)
188617 CONTINUE
1887 DO 18 j=iq+1,n
1888 g=a(ip,j)
1889 h=a(iq,j)
1890 a(ip,j)=g-s*(h+g*tau)
1891 a(iq,j)=h+s*(g-h*tau)
189218 CONTINUE
1893 DO 19 j=1,n
1894 g=v(j,ip)
1895 h=v(j,iq)
1896C
1897 v(j,ip)=g-s*(h+g*tau)
1898 v(j,iq)=h+s*(g-h*tau)
189919 CONTINUE
1900 nrot=nrot+1
1901 ENDIF
190221 CONTINUE
190322 CONTINUE
1904 DO 23 ip=1,n
1905 b(ip)=b(ip)+z(ip)
1906 d(ip)=b(ip)
1907 z(ip)=zero
190823 CONTINUE
190924 CONTINUE
1910C PAUSE 'TOO MANY ITERATIONS IN JACOBI'
1911C
1912 ier=1
1913C
1914 RETURN
1915 END
1916C
1917C=====================================================================72
1918C
1919C
1920C=====================================================================72
1921C
1922!||====================================================================
1923!|| eigsrt ../engine/source/materials/mat/mat101/sigeps101.F
1924!||--- called by ------------------------------------------------------
1925!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
1926!||====================================================================
1927 SUBROUTINE eigsrt(D,V,N,NP)
1928#include "implicit_f.inc"
1929
1930 INTEGER N,NP
1931 my_real D(NP),V(NP,NP)
1932 INTEGER I,J,K
1933 my_real p
1934 DO 13 i=1,n-1
1935 k=i
1936 p=d(i)
1937 DO 11 j=i+1,n
1938 IF(d(j)>=p)THEN
1939 k=j
1940 p=d(j)
1941 ENDIF
194211 CONTINUE
1943 IF(k/=i)THEN
1944 d(k)=d(i)
1945 d(i)=p
1946 DO 12 j=1,n
1947 p=v(j,i)
1948 v(j,i)=v(j,k)
1949 v(j,k)=p
195012 CONTINUE
1951 ENDIF
195213 CONTINUE
1953 RETURN
1954 END
1955C
1956C=====================================================================72
1957C
1958C
1959C=====================================================================72
1960C
1961!||====================================================================
1962!|| compee_tpu ../engine/source/materials/mat/mat101/sigeps101.f
1963!||--- calls -----------------------------------------------------
1964!|| at_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1965!|| atxb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1966!|| axb_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1967!|| det_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1968!|| inv_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1969!|| st_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1970!|| tr_tpu ../engine/source/materials/mat/mat101/sigeps101.F
1971!|| vec6x1tomat3x3 ../engine/source/materials/mat/mat101/sigeps101.F
1972!||====================================================================
1973 SUBROUTINE compee_tpu(F, FV, ROT, EE, NDI)
1974
1975#include "implicit_f.inc"
1976 INTEGER I, J, NDI
1977 my_real FVM(3,3), FVM_INV(3,3)
1978 my_real FE(3,3), CE(3,3), UE(3,3), F(3,3)
1979 my_real UE_INV(3,3), EE(3,3), XMA(3,3), ROT(3,3)
1980 my_real CE2(3,3), PMIDENT(3,3), TEMP(3,3)
1981 my_real FV(6)
1982 my_real INV1, INV2, INV3, B, C, M, N, T
1983 my_real L1, L2, L3, P1, P2, P3, D1, SMALL
1984 my_real X1, X2, X3, DET
1985 my_real COEFF1, COEFF2, TR_CE2
1986C
1987! DATA PI/3.141592653589793238462643D0/
1988 DATA PMIDENT /1.D0, 0.D0, 0.D0, 0.D0, 1.D0, 0.D0,
1989 & 0.D0, 0.D0, 1.D0/
1990C
1991C*********************************************************************72
1992C
1993C---- NUMERICAL CONSTANTS
1994
1995 small = em06
1996C
1997C--------------------------------------------------------------------72
1998C-------- COMPONENTS USEFUL FOR THE BRANCH A
1999C--------------------------------------------------------------------72
2000C
2001C=======DETERMINATION OF DIFFERENT TENSORS: FE->CE->UE->RE->EE
2002C
2003 CALL st_tpu(fvm, zero, 3*3)
2004 CALL st_tpu(fvm_inv, zero, 3*3)
2005 CALL st_tpu(fe, zero, 3*3)
2006 CALL st_tpu(ce, zero, 3*3)
2007 CALL st_tpu(ce2, zero, 3*3)
2008 CALL st_tpu(ue, zero, 3*3)
2009 CALL st_tpu(ue_inv, zero, 3*3)
2010 CALL st_tpu(rot, zero, 3*3)
2011C
2012C DETERMINATION OF FE AND CE
2013C
2014 CALL vec6x1tomat3x3(fv,fvm,3)
2015 CALL inv_tpu(fvm, fvm_inv, det)
2016 CALL axb_tpu(f, fvm_inv, fe, 3)
2017 CALL atxb_tpu(fe, fe, ce, 3)
2018 CALL axb_tpu(ce, ce, ce2, 3)
2019C
2020C WRITE(0 , *) 'FE IN SUBROUTINE'
2021C WRITE(0 , *) FE(1,1),' ',FE(1,2),' ',FE(1,3)
2022C WRITE(0 , *) FE(2,1),' ',FE(2,2),' ',FE(2,3)
2023C WRITE(0 , *) FE(3,1),' ',FE(3,2),' ',FE(3,3)
2024C
2025C WRITE(0 , *) 'CE IN SUBROUTINE'
2026C WRITE(0 , *) CE(1,1),' ',CE(1,2),' ',CE(1,3)
2027C WRITE(0 , *) CE(2,1),' ',CE(2,2),' ',CE(2,3)
2028C WRITE(0 , *) CE(3,1),' ',CE(3,2),' ',CE(3,3)
2029C
2030C DETERMINATION OF UE USING FRANCA ALGORITHM
2031C (SIMO & HUGHES P243)
2032C
2033 CALL tr_tpu(ce, inv1, 3)
2034 CALL tr_tpu(ce2, tr_ce2, 3)
2035 inv2 = (one/two)*(inv1**two - tr_ce2)
2036 CALL det_tpu(ce, inv3)
2037C
2038 b = inv2 - (inv1**two)/three
2039 c = -(two/27.d0)*(inv1**three) + inv1*inv2/three - inv3
2040C
2041 IF ( abs(b) < em04) THEN
2042 x1 = -abs(c)**(third)
2043 x2 = x1
2044 x3 = x1
2045 ELSE
2046 m = two*sqrt(-b/three)
2047 n = three*c/(m*b)
2048C
2049 IF (n < -one) THEN
2050 n = -one
2051 ELSE IF (n > one) THEN
2052 n = one
2053 ENDIF
2054C
2055
2056 t = atan2(sqrt(one-(n*n)),n)/three
2057 x1 = m * cos(t)
2058 x2 = m * cos(t+two*pi/three)
2059 x3 = m * cos(t+four*pi/three)
2060 ENDIF
2061C
2062 l1 = sqrt(x1 + inv1/three)
2063 l2 = sqrt(x2 + inv1/three)
2064 l3 = sqrt(x3 + inv1/three)
2065C
2066C WRITE(0, *) 'EIG OF CE ', L1,' ',L2,' ',L3
2067C
2068 p1 = l1 + l2 + l3
2069 p2 = l1*l2 + l1*l3 + l2*l3
2070 p3 = l1*l2*l3
2071C
2072 d1 = p1*p2 - p3
2073 coeff1 = p1*p1-p2
2074 coeff2 = p1*p3
2075C
2076 CALL at_tpu(coeff1, ce, coeff2, pmident, temp, 9)
2077 CALL at_tpu(-one/d1, ce2, one/d1, temp, ue, 9)
2078 CALL at_tpu(-p1, ue, p2, pmident, temp, 9)
2079 CALL at_tpu(one/p3, ce, one/p3, temp, ue_inv, 9)
2080 CALL axb_tpu(fe, ue_inv, rot, 3)
2081C
2082C WRITE(0 , *) 'ROT IN SUBR'
2083C WRITE(0 , *) ROT(1,1),' ',ROT(1,2),' ',ROT(1,3)
2084C WRITE(0 , *) ROT(2,1),' ',ROT(2,2),' ',ROT(2,3)
2085C WRITE(0 , *) ROT(3,1),' ',ROT(3,2),' ',ROT(3,3)
2086C
2087C WRITE(0 , *) 'UE IN SUBR'
2088C WRITE(0 , *) UE(1,1),' ',UE(1,2),' ',UE(1,3)
2089C WRITE(0 , *) UE(2,1),' ',UE(2,2),' ',UE(2,3)
2090C WRITE(0 , *) UE(3,1),' ',UE(3,2),' ',UE(3,3)
2091C
2092 CALL st_tpu(ee, zero, 3*3)
2093 ee(1,1) = log(ue(1,1))
2094 ee(2,2) = log(ue(2,2))
2095 ee(3,3) = log(ue(3,3))
2096C
2097C WRITE(0 , *) 'EE IN SUBR'
2098C WRITE(0 , *) EE(1,1),' ',EE(1,2),' ',EE(1,3)
2099C WRITE(0 , *) EE(2,1),' ',EE(2,2),' ',EE(2,3)
2100C WRITE(0 , *) EE(3,1),' ',EE(3,2),' ',EE(3,3)
2101C WRITE(0,*) 'SUB5'
2102C
2103 RETURN
2104 END
2105C
2106C=====================================================================72
2107C
2108C
2109C=====================================================================72
2110C
2111!||====================================================================
2112!|| invl_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2113!||--- calls -----------------------------------------------------
2114!|| err_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2115!|| signum_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2116!||====================================================================
2117 SUBROUTINE invl_tpu(A, B, C)
2118#include "implicit_f.inc"
2119
2120 my_real a, b, c
2121 my_real arg, sarg
2122C
2123C---------------------------------------------------------------------72
2124C
2125 arg = a / b
2126 IF ( abs(arg) < 0.84136d0) THEN
2127 c = 1.31446d0 * tan(1.58986d0*arg) + 0.91209d0 * arg
2128 ELSE IF (abs(arg) >= 0.84136d0 .AND. abs(arg) < one) THEN
2129 CALL signum_tpu(arg, sarg)
2130 c = one / (sarg - arg)
2131 ELSE
2132 WRITE(0 ,*) 'CHECK ERROR LANGEVIN'
2133 WRITE(0 ,*) 'LAMBDA_BAR ', a, ' COEFF ', b
2134 WRITE(0 ,*) 'ARG ', arg
2135 CALL err_tpu('ERROR IN INVERSE LANGEVIN ARGUMENT')
2136 ENDIF
2137C
2138 RETURN
2139 END
2140
2141C=====================================================================72
2142C
2143C
2144C=====================================================================72
2145C
2146!||====================================================================
2147!|| signum_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2148!||--- called by ------------------------------------------------------
2149!|| invl_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2150!||====================================================================
2151 SUBROUTINE signum_tpu(A, B )
2152#include "implicit_f.inc"
2153c
2154 my_real a, b
2155C
2156C---------------------------------------------------------------------72
2157C
2158 IF (a > 0.0) THEN
2159 b = one
2160 ELSE IF (a <zero) THEN
2161 b = -one
2162 ELSE
2163 b =zero
2164 ENDIF
2165C
2166 RETURN
2167 END
2168C
2169C=====================================================================72
2170C
2171C
2172C=====================================================================72
2173C
2174!||====================================================================
2175!|| err_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2176!||--- called by ------------------------------------------------------
2177!|| dev_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2178!|| invl_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2179!|| sigeps101 ../engine/source/materials/mat/mat101/sigeps101.F
2180!|| skw_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2181!|| sym_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2182!|| tr_tpu ../engine/source/materials/mat/mat101/sigeps101.F
2183!||--- calls -----------------------------------------------------
2184!|| arret ../engine/source/system/arret.F
2185!||====================================================================
2186 SUBROUTINE err_tpu(MESSAGE)
2187c
2188#include "implicit_f.inc"
2189 CHARACTER*(*) MESSAGE
2190C
2191C---------------------------------------------------------------------72
2192C
2193 WRITE(*, 1000) MESSAGE
2194 CALL ARRET(2)
2195C
21961000 FORMAT(/,'***ERROR MESSAGE: '/, 3X, A)
2197C
2198 END
2199
2200
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21
qtaqmt_tpu
subroutine qtaqmt_tpu(q, a, b, n)
Definition sigeps101.F:1756
sigeps101
subroutine sigeps101(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, ngl, 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, ismstr, et, ihet, offg, epsth, iexpan, tempel, fpsxx, fpsxy, fpsxz, fpsyx, fpsyy, fpsyz, fpszx, fpszy, fpszz, upsxx, upsyy, upszz, upsxy, upsyz, upsxz)
Definition sigeps101.F:67
atxb_tpu
subroutine atxb_tpu(a, b, c, n)
Definition sigeps101.F:1211
inv_tpu
subroutine inv_tpu(tensor, tensor_inv, tensor_det)
Definition sigeps101.F:1372
sst_tpu
subroutine sst_tpu(coefa, a, b, n)
Definition sigeps101.F:1241
invl_tpu
subroutine invl_tpu(a, b, c)
Definition sigeps101.F:2118
mat3x3tovec6x1
subroutine mat3x3tovec6x1(matrix, vector, n)
Definition sigeps101.F:1636
st_tpu
subroutine st_tpu(a, value, n)
Definition sigeps101.F:1144
det_tpu
subroutine det_tpu(tensor, tensor_det)
Definition sigeps101.F:1335
axb_tpu
subroutine axb_tpu(a, b, c, n)
Definition sigeps101.F:1175
eigsrt
subroutine eigsrt(d, v, n, np)
Definition sigeps101.F:1928
compee_tpu
subroutine compee_tpu(f, fv, rot, ee, ndi)
Definition sigeps101.F:1974
signum_tpu
subroutine signum_tpu(a, b)
Definition sigeps101.F:2152
valprop_tpu
subroutine valprop_tpu(a, n, np, d, v, nrot, ier)
Definition sigeps101.F:1811
tr_tpu
subroutine tr_tpu(tensor, tr_tensor, n)
Definition sigeps101.F:1306
at_tpu
subroutine at_tpu(coefa, a, coefb, b, c, n)
Definition sigeps101.F:1421
dev_tpu
subroutine dev_tpu(tensor, tensordev, n)
Definition sigeps101.F:1454
qaqt_tpu
subroutine qaqt_tpu(q, a, b, n)
Definition sigeps101.F:1497
sp_tpu
subroutine sp_tpu(a, b, c, n)
Definition sigeps101.F:1722
sym_tpu
subroutine sym_tpu(tensor, symmtensor, n)
Definition sigeps101.F:1552
emod_tpu
subroutine emod_tpu(aref, a0, c, c0, de1, de2, d1, dref, res)
Definition sigeps101.F:1784
qtaq_tpu
subroutine qtaq_tpu(q, a, b, n)
Definition sigeps101.F:1525
axbt_tpu
subroutine axbt_tpu(a, b, c, n)
Definition sigeps101.F:1271
err_tpu
subroutine err_tpu(message)
Definition sigeps101.F:2187
vec6x1tomat3x3
subroutine vec6x1tomat3x3(vector, matrix, n)
Definition sigeps101.F:1669
skw_tpu
subroutine skw_tpu(tensor, skewtensor, n)
Definition sigeps101.F:1596