OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
thcoq.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!|| thcoq ../engine/source/output/th/thcoq.F
25!||--- called by ------------------------------------------------------
26!|| hist2 ../engine/source/output/th/hist2.F
27!||--- uses -----------------------------------------------------
28!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
29!|| element_mod ../common_source/modules/elements/element_mod.f90
30!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
31!|| pinchtype_mod ../common_source/modules/pinchtype_mod.F
32!|| stack_mod ../engine/share/modules/stack_mod.F
33!||====================================================================
34 SUBROUTINE thcoq(ELBUF_TAB,MATPARAM_TAB,NTHGRP2 , ITHGRP ,
35 . IPARG ,ITHBUF ,WA ,
36 . IPM ,IGEO ,IXC ,IXTG ,PM ,
37 . RTHBUF ,THKE ,STACK)
38C-----------------------------------------------
39C M o d u l e s
40C-----------------------------------------------
41 USE elbufdef_mod
42 USE stack_mod
43 USE pinchtype_mod
44 USE matparam_def_mod
45 use element_mod , only : nixc,nixtg
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50C-----------------------------------------------
51C C o m m o n B l o c k s
52C-----------------------------------------------
53#include "mvsiz_p.inc"
54#include "com01_c.inc"
55#include "com04_c.inc"
56#include "task_c.inc"
57#include "param_c.inc"
58C-----------------------------------------------
59C D u m m y A r g u m e n t s
60C-----------------------------------------------
61 INTEGER IPARG(NPARG,*),ITHBUF(*),IXC(NIXC,*),
62 . IXTG(NIXTG,*),IPM(NPROPMI,*),IGEO(NPROPGI,*)
63 INTEGER, INTENT(in) :: NTHGRP2
64 INTEGER, DIMENSION(NITHGR,*), INTENT(in) :: ITHGRP
65 my_real
66 . wa(*),pm(npropm,*),rthbuf(*),thke(*)
67 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP), TARGET :: ELBUF_TAB
68 TYPE (MATPARAM_STRUCT_) ,DIMENSION(NUMMAT) ,INTENT(IN) :: MATPARAM_TAB
69C-----------------------------------------------
70C L o c a l V a r i a b l e s
71C-----------------------------------------------
72 INTEGER I,J,K,L,II,N, IH, NG, ITY, MTE, M5,M8,
73 . NPT,MPT,NPG,NPTR,NPTS,NPTT,NLAY,IP,IR,IS,IT,IL,IPT,
74 . NEL,NFT,I1,IUV,IAA,IADR,
75 . istrain,nu,nuvar,nuvarv,nuvard,igtyp,ihbe,
76 . ifailure,ivisc,ipmat,ishplyxfem,ipmat_iply,
77 . mat_iply,isubstack,ithk,npt_all,
78 . matly,kk(8),ipg,imat,mat_orth, idrape
79 INTEGER PID(MVSIZ),MAT(MVSIZ)
80 INTEGER :: NITER,IAD,NN,IADV,NVAR,ITYP,IJK
81 my_real :: WWA(50000),FUNC(6),SIG(5),SIGG(5)
82 my_real ,DIMENSION(MVSIZ) :: dam1,dam2,wpla,
83 . fail,fail1,fail2,fail3
84 my_real :: f1,f2,f3,f4,f5,f11,f22,f33,f44,f55,cp,sp,mm1,mm2,mm3,
85 . mm11,mm22,mm33,d1,d2,d11,d12,d22,val_ly_ip,val_ly_average
86 TYPE(g_bufel_) ,POINTER :: GBUF
87 TYPE(l_bufel_) ,POINTER :: LBUF
88 TYPE(BUF_LAY_) ,POINTER :: BUFLY
89 my_real ,DIMENSION(:), POINTER :: uvar,dir_a
90 my_real ,DIMENSION(:,:), ALLOCATABLE :: var
91 TYPE (STACK_PLY) :: STACK
92C-------------------------
93C ELEMENTS COQUES
94C=======================================================================
95 ijk = 0
96 DO niter=1,nthgrp2
97 ityp=ithgrp(2,niter)
98 nn =ithgrp(4,niter)
99 iad =ithgrp(5,niter)
100 nvar=ithgrp(6,niter)
101 iadv=ithgrp(7,niter)
102 ii=0
103 IF(ityp==3.OR.ityp==7)THEN
104! -------------------------------
105 ii=0
106 ih=iad
107C SPECIFIC SPMD
108C decalage IH
109 DO WHILE((ithbuf(ih+nn)/=ispmd).AND.(ih<iad+nn))
110 ih = ih + 1
111 ENDDO
112 IF (ih>=iad+nn) GOTO 666
113C-------------------
114 DO ng=1,ngroup
115 ity=iparg(5,ng)
116 IF (ity == ityp) THEN
117 mte=iparg(1,ng)
118 nel=iparg(2,ng)
119 nft=iparg(3,ng)
120 npt = iparg(6,ng)
121 igtyp = iparg(38,ng)
122 istrain = iparg(44,ng)
123 ihbe = iparg(23,ng)
124 ifailure = iparg(43,ng)
125 ishplyxfem = iparg(50,ng)
126 isubstack = iparg(71,ng)
127 ithk =iparg(28,ng)
128 gbuf => elbuf_tab(ng)%GBUF
129 nptr = elbuf_tab(ng)%NPTR
130 npts = elbuf_tab(ng)%NPTS
131 nptt = elbuf_tab(ng)%NPTT
132 nlay = elbuf_tab(ng)%NLAY
133 idrape = elbuf_tab(ng)%IDRAPE
134 npg = nptr*npts
135cc NPT = NLAY*NPTT ! not compatible with PID51 (shell)
136 mpt = max(1,npt)
137!
138 DO i=1,8 ! length max of GBUF%G_STRA = 8
139 kk(i) = nel*(i-1)
140 ENDDO
141!
142C
143 IF (igtyp == 51 .OR. igtyp == 52) THEN
144 npt_all = 0
145 DO ipt=1,nlay
146 npt_all = npt_all + elbuf_tab(ng)%BUFLY(ipt)%NPTT
147 ENDDO
148 IF (nlay == 1) mpt = max(1,npt_all)
149 ENDIF
150C
151 ivisc = 0
152 nuvar = 0
153 nuvarv = 0
154 nuvard = 0
155c
156 IF (mte /= 13 .and. mte /= 0) THEN
157c
158 IF ((mte>=29.AND.mte<=31).OR.
159 . mte == 35.OR.mte == 36.OR.mte == 43.OR.
160 . mte == 44.OR.mte == 45.OR.mte == 48.OR.mte>=50) THEN
161 CONTINUE
162c
163 ELSEIF (mte == 25) THEN
164C
165 DO i=1,nel
166 dam1(i)=zero
167 dam2(i)=zero
168 wpla(i)=zero
169 fail(i)=zero
170 fail1(i)=zero
171 fail2(i)=zero
172 fail3(i)=zero
173 ENDDO
174c
175 IF (ity == 3) THEN
176 DO i=1,nel
177 mat(i)=ixc(1,nft+i)
178 pid(i)=ixc(6,nft+i)
179 ENDDO
180 ELSE
181 DO i=1,nel
182 mat(i)=ixtg(1,nft+i)
183 pid(i)=ixtg(5,nft+i)
184 ENDDO
185 ENDIF
186c---
187 IF (igtyp == 11) THEN
188 ipmat = 100
189 DO n=1,mpt
190 DO i=1,nel
191 matly = igeo(ipmat+n,pid(i))
192 IF (matparam_tab(matly)%IVISC > 0) THEN
193 ivisc = 1
194 nuvarv = max(nuvarv, matparam_tab(matly)%VISC%NUVAR)
195 END IF
196 ENDDO
197 ENDDO
198 ELSEIF (igtyp == 9 .OR. igtyp == 10) THEN
199 DO n=1,mpt
200 DO i=1,nel
201 matly=mat(i)
202 IF (matparam_tab(matly)%IVISC > 0) THEN
203 ivisc = 1
204 nuvarv = max(nuvarv, matparam_tab(matly)%VISC%NUVAR)
205 END IF
206 ENDDO
207 ENDDO
208 ELSEIF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
209 ipmat = 2 + nlay
210 DO n=1,nlay
211 DO i=1,nel
212 matly = stack%IGEO(ipmat+n,isubstack)
213 IF (matparam_tab(matly)%IVISC > 0) THEN
214 ivisc = 1
215 nuvarv = max(nuvarv, matparam_tab(matly)%VISC%NUVAR)
216 END IF
217 ENDDO
218 ENDDO
219c
220 IF (ishplyxfem > 0) THEN
221 ipmat_iply = ipmat + mpt
222 DO j=1,npt -1
223 DO i=1,nel
224 mat_iply = stack%IGEO(ipmat_iply + j ,isubstack)
225 nuvard = max(nuvard, ipm(221,mat_iply))
226 ENDDO
227 ENDDO
228 ENDIF
229 ENDIF
230 ENDIF ! MTE
231c---------
232c
233c---------
234 DO i=1,nel
235 n=i+nft
236 k=ithbuf(ih)
237 ip=ithbuf(ih+nn)
238 iadr=ithbuf(ih+3*nn) ! Adress of cos sin related to skew
239C
240 IF (k == n) THEN
241 ih=ih+1
242C SPECT SPMD treatment
243C search for the correct ii
244 ii = ((ih-1) - iad)*nvar
245 DO WHILE((ithbuf(ih+nn) /= ispmd) .AND. (ih < iad+nn))
246 ih = ih + 1
247 ENDDO
248C
249 IF (ih > iad+nn) GOTO 666
250C
251 m5=5*(i-1)
252 m8=8*(i-1)
253c
254 IF (iadr /= 0) THEN ! output with respect to a (non global) SKEW
255 cp=rthbuf(iadr)
256 sp=rthbuf(iadr+1)
257c
258 f11 = gbuf%FOR(kk(1)+i)
259 f22 = gbuf%FOR(kk(2)+i)
260 f33 = gbuf%FOR(kk(3)+i)
261 f44 = gbuf%FOR(kk(4)+i)
262 f55 = gbuf%FOR(kk(5)+i)
263c
264 mm11 = gbuf%MOM(kk(1)+i)
265 mm22 = gbuf%MOM(kk(2)+i)
266 mm33 = gbuf%MOM(kk(3)+i)
267c
268 f1 = cp*cp*f11
269 . + sp*sp*f22
270 . + two*cp*sp*f33
271c
272 f2 = sp*sp*f11
273 . + cp*cp*f22
274 . - two*cp*sp*f33
275c
276 f3 =-cp*sp*f11
277 . + cp*sp*f22
278 . + (cp*cp-sp*sp )*f33
279c
280 f4 =-sp*f55+cp*f44
281 f5 = cp*f55+sp*f44
282c
283 mm1 = cp*cp*mm11
284 . + sp*sp*mm22
285 . + two*cp*sp*mm33
286c
287 mm2 = sp*sp*mm11
288 . + cp*cp*mm22
289 . - two*cp*sp*mm33
290c
291 mm3 =-cp*sp*mm11
292 . + cp*sp*mm22
293 . + (cp*cp-sp*sp )*mm33
294 ELSE !output with respect to the global SKEW.
295 f1 = gbuf%FOR(kk(1)+i)
296 f2 = gbuf%FOR(kk(2)+i)
297 f3 = gbuf%FOR(kk(3)+i)
298 f4 = gbuf%FOR(kk(4)+i)
299 f5 = gbuf%FOR(kk(5)+i)
300c
301 mm1 = gbuf%MOM(kk(1)+i)
302 mm2 = gbuf%MOM(kk(2)+i)
303 mm3 = gbuf%MOM(kk(3)+i)
304 ENDIF
305 wwa(1) = f1
306 wwa(2) = f2
307 wwa(3) = f3
308 wwa(4) = f4
309 wwa(5) = f5
310 wwa(6) = mm1
311 wwa(7) = mm2
312 wwa(8) = mm3
313 wwa(9) = gbuf%EINT(i)
314 wwa(10)= gbuf%EINT(i+nel)
315 wwa(11)= gbuf%OFF(i)
316 IF (ithk > 0) THEN
317 wwa(12)= gbuf%THK(i)
318 ELSE
319 wwa(12)= thke(n)
320 ENDIF
321 wwa(13)=zero
322 wwa(14)=zero
323 wwa(15)=zero
324 wwa(16)=zero
325 wwa(17)=zero
326 wwa(18)=zero
327 wwa(19)=zero
328 wwa(20)=zero
329 wwa(21)=zero
330 wwa(22)=zero
331 IF (gbuf%G_EPSD == 0) THEN
332 wwa(23)=zero
333 ELSE
334 wwa(23)=gbuf%EPSD(i)
335 ENDIF
336 DO j = 24,50000
337 wwa(j)=zero
338 ENDDO
339c----------------------
340c Stress tensor
341c----------------------
342c---- mean stress over Gauss points in each layer
343 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52) )THEN
344 DO il = 1,nlay
345 bufly => elbuf_tab(ng)%BUFLY(il)
346 imat = bufly%IMAT
347 nptt = bufly%NPTT
348 k = 183 + (il-1)*5
349 sigg(1:5) = zero
350 DO it=1,nptt
351 dir_a => elbuf_tab(ng)%BUFLY(il)%LBUF_DIR(it)%DIRA
352 sig(1:5) = zero
353 DO ir=1,nptr
354 DO is=1,npts
355 lbuf => bufly%LBUF(ir,is,it)
356 DO j = 1,5
357 sig(j) = sig(j) + lbuf%SIG(kk(j) + i) / npg
358 ENDDO
359 ENDDO
360 ENDDO
361 d1 = dir_a(i)
362 d2 = dir_a(i+nel)
363 d11 = d1*d1
364 d22 = d2*d2
365 d12 = d1*d2
366 sigg(1) = sigg(1) + (d11*sig(1) + d22*sig(2) + two*d12 *sig(3)) /nptt
367 sigg(2) = sigg(2) + (d22*sig(1) + d11*sig(2) - two*d12 *sig(3)) /nptt
368 sigg(3) = sigg(3) + (d12*sig(2) + (d11-d22)*sig(3) -d12*sig(1)) /nptt
369 sigg(4) = sigg(4) + (d1 *sig(4) - d2 *sig(5)) / nptt
370 sigg(5) = sigg(5) + (d1 *sig(5) + d2 *sig(4)) / nptt
371 ENDDO
372 wwa(k + 1) =sigg(1)
373 wwa(k + 2) =sigg(2)
374 wwa(k + 3) =sigg(3)
375 wwa(k + 4) =sigg(4)
376 wwa(k + 5) =sigg(5)
377 ENDDO ! DO IL=1,NLAY
378 ELSE
379 DO il = 1,nlay
380 bufly => elbuf_tab(ng)%BUFLY(il)
381 imat = bufly%IMAT
382 nptt = bufly%NPTT
383 sig(1:5) = zero
384 k = 183 + (il-1)*5
385 DO ir=1,nptr
386 DO is=1,npts
387 DO it=1,nptt
388 lbuf => bufly%LBUF(ir,is,it)
389 DO j = 1,5
390 sig(j) = sig(j) + lbuf%SIG(kk(j) + i) / (nptt*npg)
391 ENDDO
392 ENDDO
393 ENDDO
394 ENDDO
395 mat_orth = matparam_tab(imat)%ORTHOTROPY
396 IF (mat_orth == 1) THEN
397 DO j = 1,5
398 wwa(k + j) = sig(j)
399 ENDDO
400 ELSE IF (mat_orth == 2) THEN ! rotate sig to global coord
401 dir_a => elbuf_tab(ng)%BUFLY(il)%DIRA
402 d1 = dir_a(i)
403 d2 = dir_a(i+nel)
404 d11 = d1*d1
405 d22 = d2*d2
406 d12 = d1*d2
407 wwa(k + 1) = d11*sig(1) + d22*sig(2) + two*d12 *sig(3)
408 wwa(k + 2) = d22*sig(1) + d11*sig(2) - two*d12 *sig(3)
409 wwa(k + 3) =-d12*sig(1) + d12*sig(2) +(d11-d22)*sig(3)
410 wwa(k + 4) =-d2 *sig(5) + d1 *sig(4)
411 wwa(k + 5) = d1 *sig(5) + d2 *sig(4)
412 END IF
413 ENDDO ! DO IL=1,NLAY
414 ENDIF ! idrape
415c------------ Viscous stress
416c
417 DO il = 1,nlay
418 bufly => elbuf_tab(ng)%BUFLY(il)
419 imat = bufly%IMAT
420 ivisc = matparam_tab(imat)%IVISC
421 nptt = bufly%NPTT
422 IF (ivisc > 0) THEN
423 k = 30382+(il-1)*5
424 func(1:5) = zero
425 DO ir=1,nptr
426 DO is=1,npts
427 DO it=1,nptt
428 lbuf => bufly%LBUF(ir,is,it)
429 DO j = 1,5
430 func(j) = func(j) + lbuf%VISC(kk(j) + i) / nptt
431 ENDDO
432 ENDDO
433 DO j = 1,5
434 wwa(k+j) = func(j) / npg
435 ENDDO
436 ENDDO ! DO IS=1,NPTS
437 ENDDO ! DO IR=1,NPTR
438c
439 ENDIF ! IVISC > 0
440 ENDDO ! DO IL=1,NLAY
441c
442c------------ Viscous stress
443c
444c
445c------------ Max/Min Plastic strain
446 IF (gbuf%G_PLA > 0) THEN
447 wwa(13) = ep30
448 wwa(14) = zero
449c
450c IF (NPG == 1 .and. NLAY == 1) THEN
451c DO IPT = 1, NPT
452c LBUF => ELBUF_TAB(NG)%BUFLY(1)%LBUF(1,1,IPT)
453c WWA(13) = MIN(WWA(13),LBUF%PLA(I))
454c WWA(14) = MAX(WWA(14),LBUF%PLA(I))
455c ENDDO
456c ELSEIF (NLAY > 1) THEN
457c DO IPT=1,NLAY
458c BUFLY => ELBUF_TAB(NG)%BUFLY(IPT)
459c IF (BUFLY%L_PLA > 0) THEN
460c WWA(13) = MIN(WWA(13),ABS(BUFLY%PLAPT(I)))
461c WWA(14) = MAX(WWA(14),ABS(BUFLY%PLAPT(I)))
462c ENDIF
463c ENDDO
464c ELSEIF (ELBUF_TAB(NG)%BUFLY(1)%LY_PLAPT > 0) THEN
465c DO IPT = 1, NPT
466c BUFLY => ELBUF_TAB(NG)%BUFLY(1)
467c JJ = (IPT-1)*NEL
468c WWA(13) = MIN(WWA(13),ABS(BUFLY%PLAPT(JJ+I)))
469c WWA(14) = MAX(WWA(14),ABS(BUFLY%PLAPT(JJ+I)))
470c ENDDO
471c ENDIF
472 IF (nlay > 1) THEN
473 IF (npg > 1) THEN
474 DO ipt = 1,nlay
475 bufly => elbuf_tab(ng)%BUFLY(ipt)
476 IF (bufly%L_PLA > 0) THEN
477 wwa(13) = min(wwa(13),abs(bufly%PLAPT(i)))
478 wwa(14) = max(wwa(14),abs(bufly%PLAPT(i)))
479 ENDIF
480 ENDDO
481 ELSE ! (NPG = 1)
482 DO ipt = 1,nlay
483 bufly => elbuf_tab(ng)%BUFLY(ipt)
484 nptt = bufly%NPTT
485 IF (bufly%L_PLA > 0) THEN
486 func(6) = zero
487 DO it=1,nptt
488 lbuf => bufly%LBUF(1,1,it)
489 func(6) = func(6) + abs(lbuf%PLA(i))/nptt
490 ENDDO
491 wwa(13) = min(wwa(13),func(6))
492 wwa(14) = max(wwa(14),func(6))
493 ENDIF
494 ENDDO
495 ENDIF ! IF (NPG > 1) THEN
496 ELSE ! (NLAY = 1)
497 IF (npg > 1) THEN
498 bufly => elbuf_tab(ng)%BUFLY(1)
499 nptt = bufly%NPTT
500 DO it=1,nptt
501 i1 = (it-1)*nel
502 IF (bufly%L_PLA > 0) THEN
503 wwa(13) = min(wwa(13),abs(bufly%PLAPT(i1+i)))
504 wwa(14) = max(wwa(14),abs(bufly%PLAPT(i1+i)))
505 ENDIF
506 ENDDO
507 ELSE ! (NPG = 1)
508 bufly => elbuf_tab(ng)%BUFLY(1)
509 nptt = bufly%NPTT
510 DO it=1,nptt
511 lbuf => bufly%LBUF(1,1,it)
512 IF (bufly%L_PLA > 0) THEN
513 wwa(13) = min(wwa(13),abs(lbuf%PLA(i)))
514 wwa(14) = max(wwa(14),abs(lbuf%PLA(i)))
515 ENDIF
516 ENDDO
517 ENDIF
518 ENDIF ! IF (NLAY > 1) THEN
519c----
520c------------ Plastic strain per layer
521c----
522 IF (mte == 25) THEN
523 IF (ifailure == 0)THEN
524 wwa(30279) = fail(i)
525 wwa(30280) = 100*fail(i)/npt
526 wwa(30281) = fail1(i)
527 wwa(30282) = fail2(i)
528 wwa(30283) = fail3(i)
529 ENDIF ! IFAILURE == 0
530C
531 DO ipt=1,nlay
532 IF(ipt > 99) EXIT
533 bufly => elbuf_tab(ng)%BUFLY(ipt)
534 nptt = bufly%NPTT
535 val_ly_average = zero
536 DO ir=1,nptr
537 DO is=1,npts
538 val_ly_ip = zero
539 DO it=1,nptt
540 lbuf => bufly%LBUF(ir,is,it)
541 val_ly_ip = val_ly_ip + lbuf%PLA(i)/nptt
542 ENDDO
543 val_ly_average = val_ly_average + val_ly_ip/npg
544 ENDDO ! DO IS=1,NPTS
545 ENDDO ! DO IR=1,NPTR
546 wwa(30283 + ipt ) = val_ly_average
547 ENDDO ! DO IPT=1,NLAY
548 ENDIF ! MTE == 25
549 ENDIF ! GBUF%G_PLA > 0
550c----
551c------------ Non-local plastic strain and non-local plastic strain rate
552c----
553 IF (gbuf%G_PLANL > 0) THEN
554 bufly => elbuf_tab(ng)%BUFLY(1)
555 wwa(37855) = zero
556 nptt = bufly%NPTT
557 DO ir = 1,nptr
558 DO is = 1,npts
559 DO it = 1,nptt
560 wwa(37855) = wwa(37855) +
561 . bufly%LBUF(ir,is,it)%PLANL(i)/(nptr*npts*nptt)
562 ENDDO
563 ENDDO
564 ENDDO
565 ENDIF
566 IF (gbuf%G_EPSDNL > 0) THEN
567 bufly => elbuf_tab(ng)%BUFLY(1)
568 wwa(37856) = zero
569 nptt = bufly%NPTT
570 DO ir = 1,nptr
571 DO is = 1,npts
572 DO it = 1,nptt
573 wwa(37856) = wwa(37856) +
574 . bufly%LBUF(ir,is,it)%EPSDNL(i)/(nptr*npts*nptt)
575 ENDDO
576 ENDDO
577 ENDDO
578 ENDIF
579c----
580c------------ User variables
581c----
582 IF ((mte>=29.AND.mte<=31).OR.
583 . mte==35.OR.mte==36.OR.mte==43.OR.
584 . mte==44.OR.mte==45.OR.mte==48.OR.mte>=50) THEN
585c
586 nuvar = elbuf_tab(ng)%BUFLY(1)%NVAR_MAT
587 ALLOCATE (var(nuvar,max(1,mpt)))
588 var = zero
589C---
590 IF (mte == 58 .or. mte == 158) THEN
591 IF (nlay > 1) THEN
592 DO il=1,nlay
593 nuvar = elbuf_tab(ng)%BUFLY(il)%NVAR_MAT
594 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
595 DO ir=1,nptr
596 DO is=1,npts
597 k = nptr*(is-1) + ir
598 DO it=1,nptt
599 uvar=>elbuf_tab(ng)%BUFLY(il)%MAT(ir,is,it)%VAR
600 DO j = 1, nuvar
601 i1 = (j-1)*nel
602 IF (j==4 .OR. j==5) THEN ! convert to eng strain
603 var(j,il) = var(j,il) + (exp(uvar(i1+i))-one)/npg
604 ELSE
605 var(j,il) = var(j,il) + uvar(i1+i)/npg
606 ENDIF
607 wwa(6518 + (il-1)*60*4 + (k-1)*60 + j) =
608 . uvar(i1 + i)
609 ENDDO
610 ENDDO
611 ENDDO
612 ENDDO
613 ENDDO
614 ELSE ! NLAY = 1
615 nuvar = elbuf_tab(ng)%BUFLY(1)%NVAR_MAT
616 nptt = elbuf_tab(ng)%BUFLY(1)%NPTT
617 DO ipt=1,nptt
618 DO ir=1,nptr
619 DO is=1,npts
620 uvar=>elbuf_tab(ng)%BUFLY(1)%MAT(ir,is,ipt)%VAR
621 k = nptr*(is-1) + ir
622 DO j = 1, nuvar
623 i1 = (j-1)*nel
624 IF (j==4 .OR. j==5) THEN ! convert to eng strain
625 var(j,ipt) = var(j, ipt) + (exp(uvar(i1+i))-one)/npg
626 ELSE
627 var(j,ipt) = var(j, ipt) + uvar(i1 + i)/npg
628 ENDIF
629 wwa(6518 + (ipt-1)*60*4 + (k-1)*60 + j) =
630 . uvar(i1 + i)
631 ENDDO
632 ENDDO
633 ENDDO
634 ENDDO
635 ENDIF ! NLAY
636 ELSE ! (MTE /= 58)
637 IF (nlay > 1) THEN
638 DO il=1,nlay
639 nuvar = elbuf_tab(ng)%BUFLY(il)%NVAR_MAT
640 nptt = elbuf_tab(ng)%BUFLY(il)%NPTT
641 DO ir=1,nptr
642 DO is=1,npts
643 k = nptr*(is-1) + ir
644 DO it=1,nptt
645 uvar=>elbuf_tab(ng)%BUFLY(il)%MAT(ir,is,it)%VAR
646 DO j = 1, nuvar
647 i1 = (j-1)*nel
648 var(j,il) = var(j,il) + uvar(i1+i)/npg
649 wwa(6518 + (il-1)*60*4 + (k-1)*60 + j) =
650 . uvar(i1 + i)
651 ENDDO
652 ENDDO
653 ENDDO
654 ENDDO
655 ENDDO
656 ELSE ! NLAY = 1
657 nuvar = elbuf_tab(ng)%BUFLY(1)%NVAR_MAT
658 nptt = elbuf_tab(ng)%BUFLY(1)%NPTT
659 DO ipt=1,nptt
660 DO ir=1,nptr
661 DO is=1,npts
662 uvar=>elbuf_tab(ng)%BUFLY(1)%MAT(ir,is,ipt)%VAR
663 k = nptr*(is-1) + ir
664 DO j = 1, nuvar
665 i1 = (j-1)*nel
666 var(j,ipt) = var(j, ipt) + uvar(i1 + i)/npg
667 wwa(6518 + (ipt-1)*60*4 + (k-1)*60 + j) =
668 . uvar(i1 + i)
669 ENDDO
670 ENDDO
671 ENDDO
672 ENDDO
673 ENDIF ! NLAY
674 ENDIF ! MTE=58
675C---
676 nu = min(60,nuvar)
677 DO j = 1, nu
678 wwa(23+j)=var(j,iabs(mpt)/2 + 1) ! UVAR(NU) - membrane
679 DO ipt = 1, mpt
680 IF (j <= 20) THEN
681 IF (ipt <= 5) THEN
682 iuv = 83
683 iaa = 5
684 ELSEIF(ipt > 5)THEN
685 iuv = 678
686 iaa = 94
687 ENDIF
688 ELSE
689 iuv = 2558
690 iaa = 99
691 ENDIF
692 wwa(iuv + (j - 1)*iaa + ipt) = var(j, ipt) ! UVAR(NU, IPT)
693 ENDDO ! IPT =1,MPT
694 ENDDO ! J = 1, NU
695c
696 DEALLOCATE (var)
697 ENDIF ! MTE user
698c----
699c------------ Strain
700c----
701 IF (istrain /= 0) THEN
702 wwa(15)=gbuf%STRA(kk(1)+i)
703 wwa(16)=gbuf%STRA(kk(2)+i)
704 wwa(17)=gbuf%STRA(kk(3)+i)
705 wwa(18)=gbuf%STRA(kk(4)+i)
706 wwa(19)=gbuf%STRA(kk(5)+i)
707 wwa(20)=gbuf%STRA(kk(6)+i)
708 wwa(21)=gbuf%STRA(kk(7)+i)
709 wwa(22)=gbuf%STRA(kk(8)+i)
710 ENDIF
711C pinching data
712 IF(ihbe ==11.AND.npinch > 0) THEN
713 wwa(37848:37853) = zero
714 DO ipg=1,4
715 wwa(37847+1) = wwa(37847+1) + fourth*gbuf%EPGPINCHXZ(4*(i-1)+ipg)
716 wwa(37847+2) = wwa(37847+2) + fourth*gbuf%EPGPINCHYZ(4*(i-1)+ipg)
717 wwa(37847+3) = wwa(37847+3) + fourth*gbuf%EPGPINCHZZ(4*(i-1)+ipg)
718 wwa(37847+4) = wwa(37847+4) + fourth*gbuf%FORPGPINCH(4*(i-1)+ipg)
719 wwa(37847+5) = wwa(37847+5) + fourth*gbuf%MOMPGPINCH(8*(i-1)+2*(ipg-1)+1)
720 wwa(37847+6) = wwa(37847+6) + fourth*gbuf%MOMPGPINCH(8*(i-1)+2*ipg)
721 ENDDO
722 wwa(37847+7) = gbuf%THK(i)
723 ENDIF
724C end of pinching data
725c---------------------------------
726 DO l=iadv,iadv+nvar-1
727 k = ithbuf(l)
728 ijk = ijk+1
729 wa(ijk)=wwa(k)
730 ENDDO
731 ijk = ijk+1
732 wa(ijk) = ii
733c--------
734 ENDIF ! K==N
735 ENDDO ! I=1,NEL
736c--------
737 ENDIF ! MTE /= 13
738 ENDIF ! ITY == ITYP
739 ENDDO ! NG=1,NGROUP
740! -------------------------------
741 ENDIF
742 666 continue
743 ENDDO
744
745C-----------
746 RETURN
747 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
pinchtype_mod
Definition pinchtype_mod.F:69
pinchtype_mod::npinch
integer npinch
Definition pinchtype_mod.F:73
stack_mod
Definition stack_mod.F:75
thcoq
subroutine thcoq(elbuf_tab, matparam_tab, nthgrp2, ithgrp, iparg, ithbuf, wa, ipm, igeo, ixc, ixtg, pm, rthbuf, thke, stack)
Definition thcoq.F:38