OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dynain_c_strsg.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!|| dynain_c_strsg ../engine/source/output/dynain/dynain_c_strsg.F
25!||--- called by ------------------------------------------------------
26!|| gendynain ../engine/source/output/dynain/gendynain.F
27!||--- calls -----------------------------------------------------
28!|| cortdir3 ../engine/source/elements/shell/coque/cortdir3.F
29!|| layini ../engine/source/elements/shell/coque/layini.F
30!|| shell_local_frame ../engine/source/output/dynain/shell_rota.F
31!|| shell_rota ../engine/source/output/dynain/shell_rota.F
32!|| spmd_rgather9_dp ../engine/source/mpi/interfaces/spmd_outp.F
33!|| spmd_stat_pgather ../engine/source/mpi/output/spmd_stat.F
34!|| strs_txt50 ../engine/source/output/sta/sta_txt.F
35!||--- uses -----------------------------------------------------
36!|| drape_mod ../engine/share/modules/drape_mod.F
37!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
38!|| element_mod ../common_source/modules/elements/element_mod.F90
39!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
40!|| stack_mod ../engine/share/modules/stack_mod.F
41!|| state_mod ../common_source/modules/state_mod.F
42!||====================================================================
43 SUBROUTINE dynain_c_strsg(
44 1 ELBUF_TAB ,IPARG ,IGEO ,IXC ,
45 2 IXTG ,WA ,WAP0 ,IPARTC,IPARTTG,
46 3 DYNAIN_DATA,DYNAIN_INDXC,DYNAIN_INDXTG,SIZP0 ,
47 4 GEO ,STACK ,DRAPE_SH4N ,DRAPE_SH3N,X ,
48 5 THKE ,DRAPEG ,NUMMAT ,MAT_PARAM )
49C-----------------------------------------------
50C M o d u l e s
51C-----------------------------------------------
52 USE elbufdef_mod
53 USE matparam_def_mod
54 USE stack_mod
55 USE drape_mod
56 USE state_mod
57 use element_mod , only : nixc,nixtg
58C-----------------------------------------------
59C I m p l i c i t T y p e s
60C-----------------------------------------------
61#include "implicit_f.inc"
62C-----------------------------------------------
63C C o m m o n B l o c k s
64C-----------------------------------------------
65#include "com01_c.inc"
66#include "mvsiz_p.inc"
67#include "param_c.inc"
68#include "units_c.inc"
69#include "task_c.inc"
70C-----------------------------------------------
71C D u m m y A r g u m e n t s
72C-----------------------------------------------
73 INTEGER NUMMAT
74 INTEGER SIZP0
75 INTEGER IXC(NIXC,*),IXTG(NIXTG,*),
76 . IPARG(NPARG,*),IGEO(NPROPGI,*),
77 . IPARTC(*), IPARTTG(*),
78 . dynain_indxc(*), dynain_indxtg(*)
79 my_real
80 . geo(npropg,*) , x(*) ,thke(*)
81 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP), TARGET :: ELBUF_TAB
82 TYPE (STACK_PLY) :: STACK
83 TYPE (DRAPE_) :: DRAPE_SH4N(NUMELC_DRAPE),DRAPE_SH3N(NUMELTG_DRAPE)
84 TYPE (DRAPEG_) :: DRAPEG
85 double precision WA(*),WAP0(*)
86 TYPE (dynain_database), INTENT(INOUT) :: dynain_data
87 TYPE (MATPARAM_STRUCT_) ,DIMENSION(NUMMAT) ,INTENT(IN) :: MAT_PARAM
88C-----------------------------------------------
89C L o c a l V a r i a b l e s
90C-----------------------------------------------
91 INTEGER I,J,K,N,JJ,LEN,IOFF,IE,NG,NEL,NFT,LFT,NPT,
92 . llt,ity,mlw,ihbe, id, iprt,ir,is,it,
93 . npg,ipg,mpt,ipt,nptr,npts,nptt,nlay,l_pla,ithk,
94 . igtyp,npt_all,il,kk(12),large,irep,ipid,ivisc,
95 . imat,ipmat,ixfem,ixlay,isubstack,iptt,is_written,
96 , laynpt_max,nlay_max,ierr,l_dira,l_dirb,iorth,
97 . jdir,ilay,j1,j2,sedrape,numel_drape,kb
98
99 INTEGER, DIMENSION(:) , ALLOCATABLE :: MATLY !!
100 my_real, DIMENSION(:) , ALLOCATABLE :: THKLY !!
101 my_real, DIMENSION(:,:) , ALLOCATABLE :: POSLY,THK_LY
102 INTEGER , DIMENSION(:),ALLOCATABLE :: PTWA, PTWA_P0
103 INTEGER MAT(MVSIZ),PID(MVSIZ)
104 CHARACTER*80 DELIMIT
105 CHARACTER*100 LINE
106 my_real
107 . sig(6) , mom(3),a1 ,a2
108 my_real
109 . pg,mpg,qpg(2,4),
110 . sk(2),st(2),mk(2),mt(2),shk(2),sht(2),zz
111 my_real
112 . e1x(mvsiz), e1y(mvsiz), e1z(mvsiz),
113 . e2x(mvsiz), e2y(mvsiz), e2z(mvsiz),
114 . e3x(mvsiz), e3y(mvsiz), e3z(mvsiz),
115 . rx(mvsiz), ry(mvsiz), rz(mvsiz),
116 . sx(mvsiz), sy(mvsiz), sz(mvsiz),
117 . thk0(mvsiz)
118 my_real, ALLOCATABLE, DIMENSION(:) :: dira,dirb
119 my_real, DIMENSION(:) ,POINTER :: dir_a, dir_b
120 TARGET :: dira,dirb
121 TYPE(g_bufel_) ,POINTER :: GBUF
122 TYPE(L_BUFEL_) ,POINTER :: LBUF
123 TYPE(buf_lay_) ,POINTER :: BUFLY
124C-----------------------------------------------
125 parameter(pg = .577350269189626)
126 parameter(mpg=-.577350269189626)
127 DATA qpg/mpg,mpg,pg,mpg,pg,pg,mpg,pg/
128 DATA delimit(1:48)
129 ./'$--1---|---2---|---3---|---4---|---5---|---6---|'/
130 DATA delimit(49:80)
131 ./'---7---|---8---|---9---|---10--|'/
132C=======================================================================
133
134C-----------------------
135C Allocation Tabs
136C-----------------------
137 ALLOCATE(ptwa(max(dynain_data%DYNAIN_NUMELC ,
138 . dynain_data%DYNAIN_NUMELTG)),stat=ierr)
139 ALLOCATE(ptwa_p0(0:max(1,dynain_data%DYNAIN_NUMELC_G,
140 . dynain_data%DYNAIN_NUMELTG_G)),stat=ierr)
141C*********************************************
142C 4-NODE SHELLS
143C*********************************************
144 jj = 0
145C
146 ie=0
147 IF (dynain_data%DYNAIN_NUMELC/=0) THEN
148 DO ng=1,ngroup
149 ity = iparg(5,ng)
150 IF (ity == 3) THEN
151 gbuf => elbuf_tab(ng)%GBUF
152 mlw = iparg(1,ng)
153 nel = iparg(2,ng)
154 nft = iparg(3,ng)
155 mpt = iparg(6,ng)
156 ihbe = iparg(23,ng)
157 ithk = iparg(28,ng)
158 igtyp= iparg(38,ng)
159 ixfem = iparg(54,ng)
160 isubstack=iparg(71,ng)
161 ixlay = 0 ! standard element
162 ipid = ixc(6,nft+1)
163 irep = igeo(6,ipid)
164 nptr = elbuf_tab(ng)%NPTR
165 npts = elbuf_tab(ng)%NPTS
166 nptt = elbuf_tab(ng)%NPTT
167 nlay = elbuf_tab(ng)%NLAY
168 npg = nptr*npts
169 npt = nlay*nptt
170 IF (ihbe == 23) npg=4
171 lft=1
172 llt=nel
173!
174 DO i=1,12 ! length max (Hourglass for QEPH)
175 kk(i) = nel*(i-1)
176 ENDDO
177!
178C
179C pre counting of all NPTT (especially for PID_51)
180C
181 ! Npt_max
182 laynpt_max = 1
183 IF (igtyp == 51 .OR. igtyp == 52 ) THEN
184 npt_all = 0
185 DO il=1,nlay
186 npt_all = npt_all + elbuf_tab(ng)%BUFLY(il)%NPTT
187 laynpt_max = max(laynpt_max , elbuf_tab(ng)%BUFLY(il)%NPTT)
188 ENDDO
189 mpt = max(1,npt_all)
190 ENDIF
191 !
192 nlay_max = max(nlay,npt, elbuf_tab(ng)%NLAY)
193 ALLOCATE(matly(mvsiz*nlay_max), thkly(mvsiz*nlay_max*laynpt_max),
194 . posly(mvsiz,nlay_max*laynpt_max),thk_ly(nel,nlay_max*laynpt_max))
195 matly = 0
196 thkly = zero
197 posly = zero
198 thk_ly = zero
199 DO i=lft,llt
200 mat(i)=ixc(1,nft+i)
201 pid(i)=ixc(6,nft+i)
202 ENDDO
203C-------------------------------------------------
204C RELATIVE POSITION OF INTEGRATION POINTS
205C POSLY between -0.5, 0.5 : need to multiply by 2 for LSDYNA
206C------------------------------------------------
207 IF (ithk >0 ) THEN
208 thk0(lft:llt) = gbuf%THK(lft:llt)
209 ELSE
210 thk0(lft:llt) = thke(lft:llt)
211 END IF
212 numel_drape = numelc_drape
213 sedrape = scdrape
214 CALL layini(
215 . elbuf_tab(ng),lft ,llt ,geo ,igeo ,
216 . mat ,pid ,thkly ,matly ,posly ,
217 . igtyp ,ixfem ,ixlay ,nlay ,npt ,
218 . isubstack ,stack ,drape_sh4n ,nft ,thke ,
219 . nel ,thk_ly ,drapeg%INDX_SH4N ,sedrape,numel_drape)
220
221C-------------------------------------------------------
222C ELEMENT LOCAL FRAME : for rotation local -> Global
223C-------------------------------------------------------
224 CALL shell_local_frame(
225 1 lft , llt ,ity ,ihbe ,igtyp ,
226 2 ixc ,ixtg ,nft ,x ,gbuf%OFF,
227 3 rx ,ry ,rz ,sx ,sy ,
228 4 sz ,e1x ,e2x ,e3x ,e1y ,
229 5 e2y ,e3y ,e1z ,e2z ,e3z )
230
231C-------------------------------------------------------
232C ORTHOTROPY DIRECTION : for rotation orthotropic -> local
233C-------------------------------------------------------
234
235 l_dira = elbuf_tab(ng)%BUFLY(1)%LY_DIRA
236 l_dirb = elbuf_tab(ng)%BUFLY(1)%LY_DIRB
237 ALLOCATE(dira(nlay*nel*l_dira))
238 ALLOCATE(dirb(nlay*nel*l_dirb))
239 dira=zero
240 dirb=zero
241 IF (l_dira == 0) THEN
242 CONTINUE
243 ELSEIF (irep == 0) THEN
244 DO j=1,nlay
245 j1 = 1+(j-1)*l_dira*nel
246 j2 = j*l_dira*nel
247 dira(j1:j2) = elbuf_tab(ng)%BUFLY(j)%DIRA(1:nel*l_dira)
248 ENDDO
249 ENDIF
250 dir_a => dira(1:nlay*nel*l_dira)
251 dir_b => dirb(1:nlay*nel*l_dirb)
252
253 CALL cortdir3(elbuf_tab(ng),dir_a ,dir_b ,lft ,llt ,
254 . nlay ,irep ,rx ,ry ,rz ,
255 . sx ,sy ,sz ,e1x ,e1y ,
256 . e1z ,e2x ,e2y ,e2z ,nel )
257
258 iorth = 0
259 IF (igtyp /= 1) THEN
260 iorth = 1
261 ENDIF
262C-------------------------------------------------------
263C Viscosity stress add
264C-------------------------------------------------------
265
266 ivisc = 0
267 IF (igtyp == 11) THEN
268 ipmat = 100
269 DO n=1,mpt
270 DO i=lft,llt
271 imat = matly((n-1)*llt + i)
272 IF (mat_param(imat)%IVISC > 0) ivisc = 1
273 ENDDO
274 ENDDO
275 ELSEIF (igtyp == 9 .OR. igtyp == 10) THEN
276 DO n=1,mpt
277 DO i=lft,llt
278 imat = matly((n-1)*llt + i)
279 IF (mat_param(imat)%IVISC > 0) ivisc = 1
280 ENDDO
281 ENDDO
282 ELSEIF(igtyp == 17 .OR. igtyp == 51 .OR.igtyp == 52)THEN
283 ipmat = 2 + nlay
284 DO n=1,nlay
285 DO i=lft,llt
286 imat = matly((n-1)*llt + i)
287 IF (mat_param(imat)%IVISC > 0 ) ivisc = 1
288 END DO
289 ENDDO
290 ENDIF ! IGTYP
291
292C-------------------------------------------------------
293C- Loop over 4 node shell elements
294C-------------------------------------------------------
295 DO i=lft,llt
296 n = i + nft
297 iprt=ipartc(n)
298 IF (dynain_data%IPART_DYNAIN(iprt)==0) cycle
299 jj = jj + 1
300 IF (mlw /= 0 .AND. mlw /= 13) THEN
301 wa(jj) = gbuf%OFF(i)
302 ELSE
303 wa(jj) = zero
304 ENDIF
305 jj = jj + 1
306 wa(jj) = ixc(nixc,n)
307 jj = jj + 1
308 IF (mpt == 0) THEN ! global integration
309 wa(jj) = 3 ! Membrane - Lower - Upper
310 ELSE
311 wa(jj) = mpt ! Integration points
312 ENDIF
313 jj = jj + 1
314 wa(jj) = npg ! Gauss Points
315 jj = jj + 1
316 wa(jj) = one ! LARGE
317
318c---------
319
320 IF (mpt == 0) THEN ! global integration
321 IF (mlw == 0 .or. mlw == 13) THEN
322 DO ipg=1,npg
323 jj = jj + 1
324 wa(jj) = zero
325 DO j=1,8 ! forces
326 jj = jj + 1
327 wa(jj) = zero
328 ENDDO
329 ENDDO
330 ELSEIF (npg == 1) THEN
331
332! LOWER
333 a1 = one
334 a2 = -six
335
336 sig(1) = a1*gbuf%FOR(kk(1)+i) + a2* gbuf%MOM(kk(1)+i)
337 sig(2) = a1*gbuf%FOR(kk(2)+i) + a2* gbuf%MOM(kk(2)+i)
338 sig(3) = a1*gbuf%FOR(kk(3)+i) + a2* gbuf%MOM(kk(3)+i)
339 sig(4) = gbuf%FOR(kk(4)+i)
340 sig(5) = gbuf%FOR(kk(5)+i)
341 sig(6) = zero
342
343 iorth = 0
344 ilay = 0
345 jdir = 0
346
347 CALL shell_rota(
348 1 i ,ilay ,nel ,iorth ,ity ,
349 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
350 3 rx ,ry ,rz ,sx ,sy ,
351 4 sz ,e1x ,e2x ,e3x ,e1y ,
352 5 e2y ,e3y ,e1z ,e2z ,e3z ,
353 6 dir_a ,dir_b )
354
355 jj = jj + 1
356 wa(jj) = -one
357 jj = jj + 1
358 wa(jj) = sig(1)
359 jj = jj + 1
360 wa(jj) = sig(2)
361 jj = jj + 1
362 wa(jj) = sig(3)
363 jj = jj + 1
364 wa(jj) = sig(4)
365 jj = jj + 1
366 wa(jj) = sig(5)
367 jj = jj + 1
368 wa(jj) = sig(6)
369c
370 jj = jj + 1
371 IF (gbuf%G_PLA > 0) THEN
372 wa(jj) = gbuf%PLA(i)
373 ELSE
374 wa(jj) = zero
375 ENDIF
376! Membrane
377 a1 = one
378 a2 = zero
379
380 sig(1) = a1*gbuf%FOR(kk(1)+i) + a2* gbuf%MOM(kk(1)+i)
381 sig(2) = a1*gbuf%FOR(kk(2)+i) + a2* gbuf%MOM(kk(2)+i)
382 sig(3) = a1*gbuf%FOR(kk(3)+i) + a2* gbuf%MOM(kk(3)+i)
383 sig(4) = gbuf%FOR(kk(4)+i)
384 sig(5) = gbuf%FOR(kk(5)+i)
385 sig(6) = zero
386
387 iorth = 0
388 ilay = 0
389 jdir = 0
390
391 CALL shell_rota(
392 1 i ,ilay ,nel ,iorth ,ity ,
393 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
394 3 rx ,ry ,rz ,sx ,sy ,
395 4 sz ,e1x ,e2x ,e3x ,e1y ,
396 5 e2y ,e3y ,e1z ,e2z ,e3z ,
397 6 dir_a ,dir_b )
398
399 jj = jj + 1
400 wa(jj) = zero
401 jj = jj + 1
402 wa(jj) = sig(1)
403 jj = jj + 1
404 wa(jj) = sig(2)
405 jj = jj + 1
406 wa(jj) = sig(3)
407 jj = jj + 1
408 wa(jj) = sig(4)
409 jj = jj + 1
410 wa(jj) = sig(5)
411 jj = jj + 1
412 wa(jj) = sig(6)
413c
414 jj = jj + 1
415 IF (gbuf%G_PLA > 0) THEN
416 wa(jj) = gbuf%PLA(i)
417 ELSE
418 wa(jj) = zero
419 ENDIF
420! Upper
421
422 a1 = one
423 a2 = six
424
425 sig(1) = a1*gbuf%FOR(kk(1)+i) + a2* gbuf%MOM(kk(1)+i)
426 sig(2) = a1*gbuf%FOR(kk(2)+i) + a2* gbuf%MOM(kk(2)+i)
427 sig(3) = a1*gbuf%FOR(kk(3)+i) + a2* gbuf%MOM(kk(3)+i)
428 sig(4) = gbuf%FOR(kk(4)+i)
429 sig(5) = gbuf%FOR(kk(5)+i)
430 sig(6) = zero
431
432 iorth = 0
433 ilay = 0
434 jdir = 0
435
436 CALL shell_rota(
437 1 i ,ilay ,nel ,iorth ,ity ,
438 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
439 3 rx ,ry ,rz ,sx ,sy ,
440 4 sz ,e1x ,e2x ,e3x ,e1y ,
441 5 e2y ,e3y ,e1z ,e2z ,e3z ,
442 6 dir_a ,dir_b )
443
444 jj = jj + 1
445 wa(jj) = -one
446 jj = jj + 1
447 wa(jj) = sig(1)
448 jj = jj + 1
449 wa(jj) = sig(2)
450 jj = jj + 1
451 wa(jj) = sig(3)
452 jj = jj + 1
453 wa(jj) = sig(4)
454 jj = jj + 1
455 wa(jj) = sig(5)
456 jj = jj + 1
457 wa(jj) = sig(6)
458c
459 jj = jj + 1
460 IF (gbuf%G_PLA > 0) THEN
461 wa(jj) = gbuf%PLA(i)
462 ELSE
463 wa(jj) = zero
464 ENDIF
465
466 ELSE ! NPG > 1
467
468 IF(ihbe ==23) THEN !IHBE = 23 (QEPH)
469
470 st(1) = gbuf%HOURG(kk(1)+i)
471 st(2) =-gbuf%HOURG(kk(2)+i)
472 mt(1) = gbuf%HOURG(kk(3)+i)
473 mt(2) =-gbuf%HOURG(kk(4)+i)
474 sk(1) =-gbuf%HOURG(kk(7)+i)
475 sk(2) = gbuf%HOURG(kk(8)+i)
476 mk(1) =-gbuf%HOURG(kk(9)+i)
477 mk(2) = gbuf%HOURG(kk(10)+i)
478 sht(1)= gbuf%HOURG(kk(5)+i)
479 sht(2)=-gbuf%HOURG(kk(6)+i)
480 shk(1)=-gbuf%HOURG(kk(11)+i)
481 shk(2)= gbuf%HOURG(kk(12)+i)
482
483! LOWER
484 a1 = one
485 a2 = -six
486
487 DO ipg=1,npg
488
489 sig(1) = gbuf%FOR(kk(1)+i)
490 . + st(1)*qpg(2,ipg)+sk(1)*qpg(1,ipg)
491 sig(2) = gbuf%FOR(kk(2)+i)
492 . + st(2)*qpg(2,ipg)+sk(2)*qpg(1,ipg)
493 sig(3) = gbuf%FOR(kk(3)+i)
494 sig(4) = gbuf%FOR(kk(4)+i)
495 . + sht(2)*qpg(2,ipg)+shk(2)*qpg(1,ipg)
496 sig(5) = gbuf%FOR(kk(5)+i)
497 . + sht(1)*qpg(2,ipg)+shk(1)*qpg(1,ipg)
498
499 mom(1) = gbuf%MOM(kk(1)+i)
500 . + mt(1)*qpg(2,ipg)+mk(1)*qpg(1,ipg)
501 mom(2) = gbuf%MOM(kk(2)+i)
502 . + mt(2)*qpg(2,ipg)+mk(2)*qpg(1,ipg)
503 mom(3) = gbuf%MOM(kk(3)+i)
504 sig(1) = a1*sig(1) + a2*mom(1)
505 sig(2) = a1*sig(2) + a2*mom(2)
506 sig(3) = a1*sig(3) + a2*mom(3)
507 sig(6) = zero
508
509 iorth = 0
510 ilay = 0
511 jdir = 0
512
513 CALL shell_rota(
514 1 i ,ilay ,nel ,iorth ,ity ,
515 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
516 3 rx ,ry ,rz ,sx ,sy ,
517 4 sz ,e1x ,e2x ,e3x ,e1y ,
518 5 e2y ,e3y ,e1z ,e2z ,e3z ,
519 6 dir_a ,dir_b )
520
521 jj = jj + 1
522 wa(jj) = -one
523 jj = jj + 1
524 wa(jj) = sig(1)
525 jj = jj + 1
526 wa(jj) = sig(2)
527 jj = jj + 1
528 wa(jj) = sig(3)
529 jj = jj + 1
530 wa(jj) = sig(4)
531 jj = jj + 1
532 wa(jj) = sig(5)
533 jj = jj + 1
534 wa(jj) = sig(6)
535c
536 jj = jj + 1
537 IF (gbuf%G_PLA > 0) THEN
538 wa(jj) = gbuf%PLA(i)
539 ELSE
540 wa(jj) = zero
541 ENDIF
542
543 ENDDO
544
545! MEMBRANE
546 a1 = one
547 a2 = zero
548
549 DO ipg=1,npg
550
551 sig(1) = gbuf%FOR(kk(1)+i)
552 . + st(1)*qpg(2,ipg)+sk(1)*qpg(1,ipg)
553 sig(2) = gbuf%FOR(kk(2)+i)
554 . + st(2)*qpg(2,ipg)+sk(1)*qpg(1,ipg)
555 sig(3) = gbuf%FOR(kk(3)+i)
556 sig(4) = gbuf%FOR(kk(4)+i)
557 . + sht(2)*qpg(2,ipg)+shk(2)*qpg(1,ipg)
558 sig(5) = gbuf%FOR(kk(5)+i)
559 . + sht(1)*qpg(2,ipg)+shk(1)*qpg(1,ipg)
560
561 mom(1) = gbuf%MOM(kk(1)+i)
562 . + mt(1)*qpg(2,ipg)+mk(1)*qpg(1,ipg)
563 mom(2) = gbuf%MOM(kk(2)+i)
564 . + mt(2)*qpg(2,ipg)+mk(2)*qpg(1,ipg)
565 mom(3) = gbuf%MOM(kk(3)+i)
566 sig(1) = a1*sig(1) + a2*mom(1)
567 sig(2) = a1*sig(2) + a2*mom(2)
568 sig(3) = a1*sig(3) + a2*mom(3)
569 sig(6) = zero
570
571 iorth = 0
572 ilay = 0
573 jdir = 0
574
575 CALL shell_rota(
576 1 i ,ilay ,nel ,iorth ,ity ,
577 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
578 3 rx ,ry ,rz ,sx ,sy ,
579 4 sz ,e1x ,e2x ,e3x ,e1y ,
580 5 e2y ,e3y ,e1z ,e2z ,e3z ,
581 6 dir_a ,dir_b )
582
583 jj = jj + 1
584 wa(jj) = zero
585 jj = jj + 1
586 wa(jj) = sig(1)
587 jj = jj + 1
588 wa(jj) = sig(2)
589 jj = jj + 1
590 wa(jj) = sig(3)
591 jj = jj + 1
592 wa(jj) = sig(4)
593 jj = jj + 1
594 wa(jj) = sig(5)
595 jj = jj + 1
596 wa(jj) = sig(6)
597c
598 jj = jj + 1
599 IF (gbuf%G_PLA > 0) THEN
600 wa(jj) = gbuf%PLA(i)
601 ELSE
602 wa(jj) = zero
603 ENDIF
604
605 ENDDO
606! Upper
607 a1 = one
608 a2 = six
609
610 DO ipg=1,npg
611
612 sig(1) = gbuf%FOR(kk(1)+i)
613 . + st(1)*qpg(2,ipg)+sk(1)*qpg(1,ipg)
614 sig(2) = gbuf%FOR(kk(2)+i)
615 . + st(2)*qpg(2,ipg)+sk(2)*qpg(1,ipg)
616 sig(3) = gbuf%FOR(kk(3)+i)
617 sig(4) = gbuf%FOR(kk(4)+i)
618 . + sht(2)*qpg(2,ipg)+shk(2)*qpg(1,ipg)
619 sig(5) = gbuf%FOR(kk(5)+i)
620 . + sht(1)*qpg(2,ipg)+shk(1)*qpg(1,ipg)
621
622 mom(1) = gbuf%MOM(kk(1)+i)
623 . + mt(1)*qpg(2,ipg)+mk(1)*qpg(1,ipg)
624 mom(2) = gbuf%MOM(kk(2)+i)
625 . + mt(2)*qpg(2,ipg)+mk(2)*qpg(1,ipg)
626 mom(3) = gbuf%MOM(kk(3)+i)
627 sig(1) = a1*sig(1) + a2*mom(1)
628 sig(2) = a1*sig(2) + a2*mom(2)
629 sig(3) = a1*sig(3) + a2*mom(3)
630 sig(6) = zero
631
632 iorth = 0
633 ilay = 0
634 jdir = 0
635
636 CALL shell_rota(
637 1 i ,ilay ,nel ,iorth ,ity ,
638 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
639 3 rx ,ry ,rz ,sx ,sy ,
640 4 sz ,e1x ,e2x ,e3x ,e1y ,
641 5 e2y ,e3y ,e1z ,e2z ,e3z ,
642 6 dir_a ,dir_b )
643
644 jj = jj + 1
645 wa(jj) = one
646 jj = jj + 1
647 wa(jj) = sig(1)
648 jj = jj + 1
649 wa(jj) = sig(2)
650 jj = jj + 1
651 wa(jj) = sig(3)
652 jj = jj + 1
653 wa(jj) = sig(4)
654 jj = jj + 1
655 wa(jj) = sig(5)
656 jj = jj + 1
657 wa(jj) = sig(6)
658c
659 jj = jj + 1
660 IF (gbuf%G_PLA > 0) THEN
661 wa(jj) = gbuf%PLA(i)
662 ELSE
663 wa(jj) = zero
664 ENDIF
665
666 ENDDO
667
668 ELSE ! IHBE = 23
669! LOWER
670 a1 = one
671 a2 = -six
672
673 DO ir=1,nptr
674 DO is=1,npts
675 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
676 ipg = nptr*(is-1) + ir
677 k = (ipg-1)*nel*5
678 kb= (ipg-1)*nel*3
679 sig(1) = a1*gbuf%FORPG(k + kk(1) + i) + a2*gbuf%MOMPG(kb + kk(1) + i)
680 sig(2) = a1*gbuf%FORPG(k + kk(2) + i) + a2*gbuf%MOMPG(kb + kk(2) + i)
681 sig(3) = a1*gbuf%FORPG(k + kk(3) + i) + a2*gbuf%MOMPG(kb + kk(3) + i)
682 sig(4) = gbuf%FORPG(k + kk(4) + i)
683 sig(5) = gbuf%FORPG(k + kk(5) + i)
684 sig(6) = zero
685
686 iorth = 0
687 ilay = 0
688 jdir = 0
689
690 CALL shell_rota(
691 1 i ,ilay ,nel ,iorth ,ity ,
692 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
693 3 rx ,ry ,rz ,sx ,sy ,
694 4 sz ,e1x ,e2x ,e3x ,e1y ,
695 5 e2y ,e3y ,e1z ,e2z ,e3z ,
696 6 dir_a ,dir_b )
697
698 jj = jj + 1
699 wa(jj) = -one
700 jj = jj + 1
701 wa(jj) = sig(1)
702 jj = jj + 1
703 wa(jj) = sig(2)
704 jj = jj + 1
705 wa(jj) = sig(3)
706 jj = jj + 1
707 wa(jj) = sig(4)
708 jj = jj + 1
709 wa(jj) = sig(5)
710 jj = jj + 1
711 wa(jj) = sig(6)
712c
713 jj = jj + 1
714 IF (gbuf%G_PLA > 0) THEN
715 wa(jj) = gbuf%PLA(i)
716 ELSE
717 wa(jj) = zero
718 ENDIF
719 ENDDO
720 ENDDO
721
722! MEMBRANE
723 a1 = one
724 a2 = zero
725
726 DO ir=1,nptr
727 DO is=1,npts
728 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
729 ipg = nptr*(is-1) + ir
730 k = (ipg-1)*nel*5
731 kb = (ipg-1)*nel*3
732
733 sig(1) = a1*gbuf%FORPG(k + kk(1) + i) + a2*gbuf%MOMPG(kb + kk(1) + i)
734 sig(2) = a1*gbuf%FORPG(k + kk(2) + i) + a2*gbuf%MOMPG(kb + kk(2) + i)
735 sig(3) = a1*gbuf%FORPG(k + kk(3) + i) + a2*gbuf%MOMPG(kb + kk(3) + i)
736 sig(4) = gbuf%FORPG(k + kk(4) + i)
737 sig(5) = gbuf%FORPG(k + kk(5) + i)
738 sig(6) = zero
739
740 iorth = 0
741 ilay = 0
742 jdir = 0
743
744 CALL shell_rota(
745 1 i ,ilay ,nel ,iorth ,ity ,
746 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
747 3 rx ,ry ,rz ,sx ,sy ,
748 4 sz ,e1x ,e2x ,e3x ,e1y ,
749 5 e2y ,e3y ,e1z ,e2z ,e3z ,
750 6 dir_a ,dir_b )
751
752 jj = jj + 1
753 wa(jj) = zero
754 jj = jj + 1
755 wa(jj) = sig(1)
756 jj = jj + 1
757 wa(jj) = sig(2)
758 jj = jj + 1
759 wa(jj) = sig(3)
760 jj = jj + 1
761 wa(jj) = sig(4)
762 jj = jj + 1
763 wa(jj) = sig(5)
764 jj = jj + 1
765 wa(jj) = sig(6)
766c
767 jj = jj + 1
768 IF (gbuf%G_PLA > 0) THEN
769 wa(jj) = gbuf%PLA(i)
770 ELSE
771 wa(jj) = zero
772 ENDIF
773 ENDDO
774 ENDDO
775! Upper
776 a1 = one
777 a2 = six
778
779 DO ir=1,nptr
780 DO is=1,npts
781 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
782 ipg = nptr*(is-1) + ir
783 k = (ipg-1)*nel*5
784 kb =(ipg-1)*nel*3
785
786 sig(1) = a1*gbuf%FORPG(k + kk(1) + i) + a2*gbuf%MOMPG(kb + kk(1) + i)
787 sig(2) = a1*gbuf%FORPG(k + kk(2) + i) + a2*gbuf%MOMPG(kb + kk(2) + i)
788 sig(3) = a1*gbuf%FORPG(k + kk(3) + i) + a2*gbuf%MOMPG(kb + kk(3) + i)
789 sig(4) = gbuf%FORPG(k + kk(4) + i)
790 sig(5) = gbuf%FORPG(k + kk(5) + i)
791 sig(6) = zero
792
793 iorth = 0
794 ilay = 0
795 jdir = 0
796
797 CALL shell_rota(
798 1 i ,ilay ,nel ,iorth ,ity ,
799 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
800 3 rx ,ry ,rz ,sx ,sy ,
801 4 sz ,e1x ,e2x ,e3x ,e1y ,
802 5 e2y ,e3y ,e1z ,e2z ,e3z ,
803 6 dir_a ,dir_b )
804
805 jj = jj + 1
806 wa(jj) = one
807 jj = jj + 1
808 wa(jj) = sig(1)
809 jj = jj + 1
810 wa(jj) = sig(2)
811 jj = jj + 1
812 wa(jj) = sig(3)
813 jj = jj + 1
814 wa(jj) = sig(4)
815 jj = jj + 1
816 wa(jj) = sig(5)
817 jj = jj + 1
818 wa(jj) = sig(6)
819c
820 jj = jj + 1
821 IF (gbuf%G_PLA > 0) THEN
822 wa(jj) = gbuf%PLA(i)
823 ELSE
824 wa(jj) = zero
825 ENDIF
826 ENDDO
827 ENDDO
828
829 ENDIF ! IHBE = 23
830
831 ENDIF ! IF (MLW == 0 .or. MLW == 13)
832C-------------------------
833C (MPT /=0 ):
834C-------------------------
835 ELSEIF (mlw == 0 .or. mlw == 13) THEN
836
837 DO k=1,mpt
838 DO ipg=1,npg
839 DO j=1,8 ! Stress + pos + plastic strain
840 jj = jj + 1
841 wa(jj) = zero
842 ENDDO
843 ENDDO
844 ENDDO
845
846 ELSEIF(ihbe == 23) THEN
847
848 st(1) = gbuf%HOURG(kk(1)+i)
849 st(2) =-gbuf%HOURG(kk(2)+i)
850 mt(1) = gbuf%HOURG(kk(3)+i)
851 mt(2) =-gbuf%HOURG(kk(4)+i)
852 sk(1) =-gbuf%HOURG(kk(7)+i)
853 sk(2) = gbuf%HOURG(kk(8)+i)
854 mk(1) =-gbuf%HOURG(kk(9)+i)
855 mk(2) = gbuf%HOURG(kk(10)+i)
856 sht(1)= gbuf%HOURG(kk(5)+i)
857 sht(2)=-gbuf%HOURG(kk(6)+i)
858 shk(1)=-gbuf%HOURG(kk(11)+i)
859 shk(2)= gbuf%HOURG(kk(12)+i)
860
861 iptt = 0
862 DO il = 1,nlay
863 bufly => elbuf_tab(ng)%BUFLY(il)
864 nptt = bufly%NPTT
865
866 jdir = 1 + (il-1)*nel*2
867
868 DO it=1,nptt
869 ipt = iptt + it
870 lbuf => bufly%LBUF(1,1,it)
871 l_pla = bufly%L_PLA
872 zz = posly(i,ipt)*thk0(i)
873 DO ipg=1,npg
874
875 sig(1) = lbuf%SIG(kk(1)+i)
876 . + (st(1)+zz*mt(1))*qpg(2,ipg)
877 . + (sk(1)+zz*mk(1))*qpg(1,ipg)
878 sig(2) = lbuf%SIG(kk(2)+i)
879 . + (st(2)+zz*mt(2))*qpg(2,ipg)
880 . + (sk(2)+zz*mk(2))*qpg(1,ipg)
881 sig(3) = lbuf%SIG(kk(3)+i)
882 sig(4) = lbuf%SIG(kk(4)+i)
883 . + sht(2)*qpg(2,ipg)+shk(2)*qpg(1,ipg)
884 sig(5) = lbuf%SIG(kk(5)+i)
885 . + sht(1)*qpg(2,ipg)+shk(1)*qpg(1,ipg)
886 sig(6) = zero
887
888 CALL shell_rota(
889 1 i ,il ,nel ,iorth ,ity ,
890 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
891 3 rx ,ry ,rz ,sx ,sy ,
892 4 sz ,e1x ,e2x ,e3x ,e1y ,
893 5 e2y ,e3y ,e1z ,e2z ,e3z ,
894 6 dir_a ,dir_b )
895
896 jj = jj + 1
897 wa(jj) = two * posly(i,ipt)
898
899 jj = jj + 1
900 wa(jj) = sig(1)
901 jj = jj + 1
902 wa(jj) = sig(2)
903 jj = jj + 1
904 wa(jj) = sig(3)
905 jj = jj + 1
906 wa(jj) = sig(4)
907 jj = jj + 1
908 wa(jj) = sig(5)
909 jj = jj + 1
910 wa(jj) = sig(6)
911 jj = jj + 1
912 IF (bufly%L_PLA > 0) THEN
913 wa(jj) = lbuf%PLA(i)
914 ELSE
915 wa(jj) = zero
916 ENDIF
917
918 ENDDO ! DO IPG=1,NPG
919 ENDDO ! DO IT=1,NPTT
920 iptt= iptt + nptt
921
922 ENDDO ! DO IL=1,NLAY
923
924
925 ELSEIF (nlay == 1) THEN ! PID1
926
927 bufly => elbuf_tab(ng)%BUFLY(1)
928 nptt = bufly%NPTT
929
930 DO it=1,nptt
931 DO is=1,npts
932 DO ir=1,nptr
933 lbuf => bufly%LBUF(ir,is,it)
934 ipg = nptr*(is-1) + ir
935 sig(1) = lbuf%SIG(kk(1)+i)
936 sig(2) = lbuf%SIG(kk(2)+i)
937 sig(3) = lbuf%SIG(kk(3)+i)
938 sig(4) = lbuf%SIG(kk(4)+i)
939 sig(5) = lbuf%SIG(kk(5)+i)
940 sig(6) = zero
941C
942C viscous stress added
943C
944 IF (ivisc > 0) THEN
945
946 sig(1) = sig(1) + lbuf%VISC(kk(1)+i)
947 sig(2) = sig(2) + lbuf%VISC(kk(2)+i)
948 sig(3) = sig(3) + lbuf%VISC(kk(3)+i)
949 sig(4) = sig(4) + lbuf%VISC(kk(4)+i)
950 sig(5) = sig(5) + lbuf%VISC(kk(5)+i)
951
952 ENDIF ! IF (IVISC > 0)
953
954 ilay = 1
955 jdir = 1 + (ilay-1)*llt*2
956
957 CALL shell_rota(
958 1 i ,ilay ,nel ,iorth ,ity ,
959 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
960 3 rx ,ry ,rz ,sx ,sy ,
961 4 sz ,e1x ,e2x ,e3x ,e1y ,
962 5 e2y ,e3y ,e1z ,e2z ,e3z ,
963 6 dir_a ,dir_b )
964
965
966
967 jj = jj + 1
968 wa(jj) = two * posly(i,it)
969
970 jj = jj + 1
971 wa(jj) = sig(1)
972 jj = jj + 1
973 wa(jj) = sig(2)
974 jj = jj + 1
975 wa(jj) = sig(3)
976 jj = jj + 1
977 wa(jj) = sig(4)
978 jj = jj + 1
979 wa(jj) = sig(5)
980 jj = jj + 1
981 wa(jj) = sig(6)
982 jj = jj + 1
983 IF (bufly%L_PLA > 0) THEN
984 wa(jj) = lbuf%PLA(i)
985 ELSE
986 wa(jj) = zero
987 ENDIF
988
989 ENDDO
990 ENDDO
991 ENDDO ! DO IPT = 1,NPTT
992
993
994 ELSE ! NLAY > 1, PID10,PID11,PID16,PID17,PID51,PID52
995
996 iptt = 0
997 DO il = 1,nlay
998 bufly => elbuf_tab(ng)%BUFLY(il)
999 nptt = bufly%NPTT
1000
1001 jdir = 1 + (il-1)*llt*2
1002
1003 DO it=1,nptt
1004 ipt = iptt + it
1005 DO is=1,npts
1006 DO ir=1,nptr
1007 lbuf => bufly%LBUF(ir,is,it)
1008
1009 sig(1) = lbuf%SIG(kk(1)+i)
1010 sig(2) = lbuf%SIG(kk(2)+i)
1011 sig(3) = lbuf%SIG(kk(3)+i)
1012 sig(4) = lbuf%SIG(kk(4)+i)
1013 sig(5) = lbuf%SIG(kk(5)+i)
1014 sig(6) = zero
1015C
1016C viscous stress added
1017C
1018 IF (ivisc > 0) THEN
1019 sig(1) = sig(1) + lbuf%VISC(kk(1)+i)
1020 sig(2) = sig(2) + lbuf%VISC(kk(2)+i)
1021 sig(3) = sig(3) + lbuf%VISC(kk(3)+i)
1022 sig(4) = sig(4) + lbuf%VISC(kk(4)+i)
1023 sig(5) = sig(5) + lbuf%VISC(kk(5)+i)
1024 ENDIF ! IF (IVISC > 0)
1025
1026 CALL shell_rota(
1027 1 i ,il ,nel ,iorth ,ity ,
1028 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1029 3 rx ,ry ,rz ,sx ,sy ,
1030 4 sz ,e1x ,e2x ,e3x ,e1y ,
1031 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1032 6 dir_a ,dir_b )
1033
1034 jj = jj + 1
1035 wa(jj) = two * posly(i,ipt)
1036
1037 jj = jj + 1
1038 wa(jj) = sig(1)
1039 jj = jj + 1
1040 wa(jj) = sig(2)
1041 jj = jj + 1
1042 wa(jj) = sig(3)
1043 jj = jj + 1
1044 wa(jj) = sig(4)
1045 jj = jj + 1
1046 wa(jj) = sig(5)
1047 jj = jj + 1
1048 wa(jj) = sig(6)
1049 jj = jj + 1
1050 IF (bufly%L_PLA > 0) THEN
1051 wa(jj) = lbuf%PLA(i)
1052 ELSE
1053 wa(jj) = zero
1054 ENDIF
1055 ENDDO
1056 ENDDO
1057 ENDDO
1058 iptt = iptt + nptt
1059 ENDDO
1060 ENDIF ! MPT, NLAY
1061C
1062 ie=ie+1
1063C Pointer last position of shell IE in WA
1064 ptwa(ie)=jj
1065 ENDDO ! DO I=LFT,LLT
1066
1067 IF (ALLOCATED(dirb)) DEALLOCATE(dirb)
1068 IF (ALLOCATED(dira)) DEALLOCATE(dira)
1069c------- end loop over 4 node shell elements
1070 DEALLOCATE(matly, thkly, posly, thk_ly)
1071 ENDIF ! ITY == 3
1072 ENDDO ! NG = 1, NGROUP
1073 ENDIF ! DYNAIN_NUMEL /= 0
1074
1075
1076c-----------------------------------------------------------------------
1077c 4N SHELLS - WRITE
1078c-----------------------------------------------------------------------
1079 IF (nspmd == 1) THEN
1080C recopying for code simplification
1081 ptwa_p0(0)=0
1082 DO n=1,dynain_data%DYNAIN_NUMELC
1083 ptwa_p0(n)=ptwa(n)
1084 ENDDO
1085 len=jj
1086 DO j=1,len
1087 wap0(j)=wa(j)
1088 ENDDO
1089 ELSE
1090C Global pointers WAP0
1091 CALL spmd_stat_pgather(ptwa,dynain_data%DYNAIN_NUMELC,ptwa_p0,dynain_data%DYNAIN_NUMELC_G)
1092 len = 0
1093 CALL spmd_rgather9_dp(wa,jj,wap0,sizp0,len)
1094 ENDIF
1095c-------------------------------------
1096 is_written = 0
1097 IF (ispmd == 0.AND.len > 0) THEN
1098 IF(dynain_data%ZIPDYNAIN==0) THEN
1099 WRITE(iudynain,'(A)') delimit
1100 WRITE(iudynain,'(A)')'*INITIAL_STRESS_SHELL'
1101 WRITE(iudynain,'(A)')
1102 . '$ SHELLID NPG NBINT LARGE '
1103 WRITE(iudynain,'(A)')
1104 . '$ IF(NPT == 0), REPEAT I=1,NPG :'
1105 WRITE(iudynain,'(A)')
1106 . '$ IF(NPT /= 0) REPEAT K=1,NPT : REPEAT I=1,NPG :'
1107 WRITE(iudynain,'(A)')
1108 . '$ T SIGXX SIGYY SIGZZ SIGXY SIGYZ SIGZX EPSP '
1109 WRITE(iudynain,'(A)') delimit
1110 ELSE
1111 WRITE(line,'(A)') delimit
1112 CALL strs_txt50(line,100)
1113 WRITE(line,'(A)')'*INITIAL_STRESS_SHELL'
1114 CALL strs_txt50(line,100)
1115 WRITE(line,'(A)')
1116 . '$ SHELLID NPG NBINT LARGE '
1117 CALL strs_txt50(line,100)
1118 WRITE(line,'(A)')
1119 . '$ IF(NPT == 0), REPEAT I=1,NPG :'
1120 CALL strs_txt50(line,100)
1121 WRITE(line,'(A)')
1122 . '$ IF(NPT /= 0) REPEAT K=1,NPT : REPEAT I=1,NPG :'
1123 CALL strs_txt50(line,100)
1124 WRITE(line,'(A)')
1125 . '$ T SIGXX SIGYY SIGZZ SIGXY SIGYZ SIGZX EPSP '
1126 CALL strs_txt50(line,100)
1127 WRITE(line,'(A)') delimit
1128 CALL strs_txt50(line,100)
1129
1130 ENDIF
1131 is_written = 1
1132 DO n=1,dynain_data%DYNAIN_NUMELC_G
1133C Retrieving shell ID in increasing order
1134 k=dynain_indxc(n)
1135C Adress in WAP0
1136 j=ptwa_p0(k-1)
1137
1138 ioff = nint(wap0(j + 1))
1139 IF (ioff >= 1) THEN
1140c
1141 id = nint(wap0(j + 2))
1142 npt = nint(wap0(j + 3))
1143 npg = nint(wap0(j + 4))
1144 large = nint(wap0(j + 5))
1145
1146 j = j + 5
1147 IF(dynain_data%ZIPDYNAIN==0) THEN
1148 WRITE(iudynain,'(3I8,16X,I8)')id,npg,npt,large
1149 ELSE
1150 WRITE(line,'(3I8,16X,I8)')id,npg,npt,large
1151 CALL strs_txt50(line,100)
1152 ENDIF
1153 IF (npt == 0) THEN
1154 DO ipg=1,npg
1155 IF(dynain_data%ZIPDYNAIN==0) THEN
1156 WRITE(iudynain,'(1P5G16.9)')(wap0(jj + k),k=1,5)
1157 WRITE(iudynain,'(1P3G16.9)')(wap0(jj + k),k=6,8)
1158 ELSE
1159 WRITE(line,'(1P5G16.9)')(wap0(jj + k),k=1,5)
1160 CALL strs_txt50(line,100)
1161 WRITE(line,'(1P3G16.9)')(wap0(jj + k),k=6,8)
1162 CALL strs_txt50(line,100)
1163 ENDIF
1164 j = j + 8
1165 ENDDO
1166 ELSE
1167 DO ipt=1,npt
1168 DO ipg=1,npg
1169 IF(dynain_data%ZIPDYNAIN==0) THEN
1170 WRITE(iudynain,'(1P5G16.9)')(wap0(j + k),k=1,5)
1171 WRITE(iudynain,'(1P3G16.9)')(wap0(j + k),k=6,8)
1172 ELSE
1173 WRITE(line,'(1P5G16.9)')(wap0(j + k),k=1,5)
1174 CALL strs_txt50(line,100)
1175 WRITE(line,'(1P3G16.9)')(wap0(j + k),k=6,8)
1176 CALL strs_txt50(line,100)
1177 ENDIF
1178 j = j + 8
1179 ENDDO
1180 ENDDO
1181
1182 ENDIF ! IF (NPT == 0)
1183 ENDIF ! IF (IOFF >= 1)
1184 ENDDO ! DO N=1,DYNAIN_NUMELC_G
1185 ENDIF ! IF (ISPMD == 0.AND.LEN > 0)
1186
1187C***********************************************
1188C 3-NODE SHELLS
1189C***********************************************
1190 jj = 0
1191 ie=0
1192C
1193 IF(dynain_data%DYNAIN_NUMELTG/=0) THEN
1194 DO ng=1,ngroup
1195 ity = iparg(5,ng)
1196 IF (ity == 7) THEN
1197 gbuf => elbuf_tab(ng)%GBUF
1198 mlw = iparg(1,ng)
1199 nel = iparg(2,ng)
1200 nft = iparg(3,ng)
1201 mpt = iparg(6,ng)
1202 ihbe = iparg(23,ng)
1203 ithk = iparg(28,ng)
1204 igtyp= iparg(38,ng)
1205 ipid = ixtg(5,nft+1)
1206 irep = igeo(6,ipid)
1207 nptr = elbuf_tab(ng)%NPTR
1208 npts = elbuf_tab(ng)%NPTS
1209 nptt = elbuf_tab(ng)%NPTT
1210 nlay = elbuf_tab(ng)%NLAY
1211 npg = nptr*npts
1212 npt = nlay*nptt
1213 lft=1
1214 llt=nel
1215!
1216 DO i=1,6
1217 kk(i) = nel*(i-1)
1218 ENDDO
1219!
1220C
1221C pre counting of all NPTT (especially for PID_51)
1222C
1223 ! Npt_max
1224 laynpt_max = 1
1225 IF (igtyp == 51 .OR. igtyp == 52 ) THEN
1226 npt_all = 0
1227 DO k=1,nlay
1228 npt_all = npt_all + elbuf_tab(ng)%BUFLY(k)%NPTT
1229 laynpt_max = max(laynpt_max , elbuf_tab(ng)%BUFLY(k)%NPTT)
1230 ENDDO
1231 mpt = max(1,npt_all)
1232 ENDIF
1233 !
1234 nlay_max = max(nlay,npt, elbuf_tab(ng)%NLAY)
1235 ALLOCATE(matly(mvsiz*nlay_max), thkly(mvsiz*nlay_max*laynpt_max),
1236 . posly(mvsiz,nlay_max*laynpt_max),thk_ly(nel,nlay_max*laynpt_max))
1237 matly = 0
1238 thkly = zero
1239 posly = zero
1240 thk_ly = zero
1241 DO i=lft,llt
1242 mat(i)=ixtg(1,nft+i)
1243 pid(i)=ixtg(5,nft+i)
1244 ENDDO
1245C-------------------------------------------------
1246C RELATIVE POSITION OF INTEGRATION POINTS
1247C POSLY between -0.5, 0.5 : need to multiply by 2 for LSDYNA
1248C------------------------------------------------
1249 IF (ithk >0 ) THEN
1250 thk0(lft:llt) = gbuf%THK(lft:llt)
1251 ELSE
1252 thk0(lft:llt) = thke(lft:llt)
1253 END IF
1254 numel_drape = numeltg_drape
1255 sedrape = stdrape
1256 CALL layini(
1257 . elbuf_tab(ng),lft ,llt ,geo ,igeo ,
1258 . mat ,pid ,thkly ,matly ,posly ,
1259 . igtyp ,ixfem ,ixlay ,nlay ,npt ,
1260 . isubstack ,stack ,drape_sh3n ,nft ,thke ,
1261 . nel ,thk_ly ,drapeg%INDX_SH3N,sedrape,numel_drape)
1262
1263C-------------------------------------------------------
1264C ELEMENT LOCAL FRAME : for rotation local -> Global
1265C-------------------------------------------------------
1266 CALL shell_local_frame(
1267 1 lft , llt ,ity ,ihbe ,igtyp ,
1268 2 ixc ,ixtg ,nft ,x ,gbuf%OFF,
1269 3 rx ,ry ,rz ,sx ,sy ,
1270 4 sz ,e1x ,e2x ,e3x ,e1y ,
1271 5 e2y ,e3y ,e1z ,e2z ,e3z )
1272
1273C-------------------------------------------------------
1274C ORTHOTROPY DIRECTION : for rotation orthotropic -> local
1275C-------------------------------------------------------
1276
1277 l_dira = elbuf_tab(ng)%BUFLY(1)%LY_DIRA
1278 l_dirb = elbuf_tab(ng)%BUFLY(1)%LY_DIRB
1279 ALLOCATE(dira(nlay*nel*l_dira))
1280 ALLOCATE(dirb(nlay*nel*l_dirb))
1281 dira=zero
1282 dirb=zero
1283 IF (l_dira == 0) THEN
1284 CONTINUE
1285 ELSEIF (irep == 0) THEN
1286 DO j=1,nlay
1287 j1 = 1+(j-1)*l_dira*nel
1288 j2 = j*l_dira*nel
1289 dira(j1:j2) = elbuf_tab(ng)%BUFLY(j)%DIRA(1:nel*l_dira)
1290 ENDDO
1291 ENDIF
1292 dir_a => dira(1:nlay*nel*l_dira)
1293 dir_b => dirb(1:nlay*nel*l_dirb)
1294
1295 CALL cortdir3(elbuf_tab(ng),dir_a ,dir_b ,lft ,llt ,
1296 . nlay ,irep ,rx ,ry ,rz ,
1297 . sx ,sy ,sz ,e1x ,e1y ,
1298 . e1z ,e2x ,e2y ,e2z ,nel )
1299
1300 iorth = 1
1301C-------------------------------------------------------
1302C Viscosity stress add
1303C-------------------------------------------------------
1304
1305 ivisc = 0
1306 IF (igtyp == 11) THEN
1307 ipmat = 100
1308 DO n=1,mpt
1309 DO i=lft,llt
1310 imat = matly((n-1)*llt + i)
1311 IF (mat_param(imat)%IVISC > 0 ) ivisc = 1
1312 ENDDO
1313 ENDDO
1314 ELSEIF (igtyp == 9 .OR. igtyp == 10) THEN
1315 DO n=1,mpt
1316 DO i=lft,llt
1317 imat = matly((n-1)*llt + i)
1318 IF (mat_param(imat)%IVISC > 0 ) ivisc = 1
1319 ENDDO
1320 ENDDO
1321 ELSEIF(igtyp == 17 .OR. igtyp == 51 .OR.igtyp == 52)THEN
1322 ipmat = 2 + nlay
1323 DO n=1,nlay
1324 DO i=lft,llt
1325 imat = matly((n-1)*llt + i)
1326 IF (mat_param(imat)%IVISC > 0 ) ivisc = 1
1327 END DO
1328 ENDDO
1329 ENDIF ! IGTYP
1330
1331C-------------------------------------------------------
1332C- Loop over 3 node shell elements
1333C-------------------------------------------------------
1334
1335 DO i=lft,llt
1336 n = i + nft
1337 iprt=iparttg(n)
1338 IF (dynain_data%IPART_DYNAIN(iprt) == 0) cycle
1339 jj = jj + 1
1340 IF (mlw /= 0 .AND. mlw /= 13) THEN
1341 wa(jj) = gbuf%OFF(i)
1342 ELSE
1343 wa(jj) = zero
1344 ENDIF
1345 jj = jj + 1
1346 wa(jj) = ixtg(nixtg,n)
1347 jj = jj + 1
1348 IF (mpt == 0) THEN ! global integration
1349 wa(jj) = 3 ! Membrane - Lower - Upper
1350 ELSE
1351 wa(jj) = mpt ! Integration points
1352 ENDIF
1353 jj = jj + 1
1354 wa(jj) = npg
1355 jj = jj + 1
1356 wa(jj) = one ! LARGE
1357c----
1358 IF (mpt == 0) THEN ! global integration
1359 IF (mlw == 0 .or. mlw == 13) THEN
1360 DO ipg=1,npg
1361 DO j=1,9
1362 jj = jj + 1
1363 wa(jj) = zero
1364 ENDDO
1365 ENDDO
1366 ELSEIF (npg == 1) THEN
1367
1368! LOWER
1369 a1 = one
1370 a2 = -six
1371
1372 sig(1) = a1*gbuf%FOR(kk(1)+i) + a2* gbuf%MOM(kk(1)+i)
1373 sig(2) = a1*gbuf%FOR(kk(2)+i) + a2* gbuf%MOM(kk(2)+i)
1374 sig(3) = a1*gbuf%FOR(kk(3)+i) + a2* gbuf%MOM(kk(3)+i)
1375 sig(4) = gbuf%FOR(kk(4)+i)
1376 sig(5) = gbuf%FOR(kk(5)+i)
1377 sig(6) = zero
1378
1379 iorth = 0
1380 ilay = 0
1381 jdir = 0
1382
1383 CALL shell_rota(
1384 1 i ,ilay ,nel ,iorth ,ity ,
1385 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1386 3 rx ,ry ,rz ,sx ,sy ,
1387 4 sz ,e1x ,e2x ,e3x ,e1y ,
1388 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1389 6 dir_a ,dir_b )
1390
1391 jj = jj + 1
1392 wa(jj) = -one
1393 jj = jj + 1
1394 wa(jj) = sig(1)
1395 jj = jj + 1
1396 wa(jj) = sig(2)
1397 jj = jj + 1
1398 wa(jj) = sig(3)
1399 jj = jj + 1
1400 wa(jj) = sig(4)
1401 jj = jj + 1
1402 wa(jj) = sig(5)
1403 jj = jj + 1
1404 wa(jj) = sig(6)
1405c
1406 jj = jj + 1
1407 IF (gbuf%G_PLA > 0) THEN
1408 wa(jj) = gbuf%PLA(i)
1409 ELSE
1410 wa(jj) = zero
1411 ENDIF
1412! Membrane
1413 a1 = one
1414 a2 = zero
1415
1416 sig(1) = a1*gbuf%FOR(kk(1)+i) + a2* gbuf%MOM(kk(1)+i)
1417 sig(2) = a1*gbuf%FOR(kk(2)+i) + a2* gbuf%MOM(kk(2)+i)
1418 sig(3) = a1*gbuf%FOR(kk(3)+i) + a2* gbuf%MOM(kk(3)+i)
1419 sig(4) = gbuf%FOR(kk(4)+i)
1420 sig(5) = gbuf%FOR(kk(5)+i)
1421 sig(6) = zero
1422
1423 iorth = 0
1424 ilay = 0
1425 jdir = 0
1426
1427 CALL shell_rota(
1428 1 i ,ilay ,nel ,iorth ,ity ,
1429 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1430 3 rx ,ry ,rz ,sx ,sy ,
1431 4 sz ,e1x ,e2x ,e3x ,e1y ,
1432 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1433 6 dir_a ,dir_b )
1434
1435 jj = jj + 1
1436 wa(jj) = zero
1437 jj = jj + 1
1438 wa(jj) = sig(1)
1439 jj = jj + 1
1440 wa(jj) = sig(2)
1441 jj = jj + 1
1442 wa(jj) = sig(3)
1443 jj = jj + 1
1444 wa(jj) = sig(4)
1445 jj = jj + 1
1446 wa(jj) = sig(5)
1447 jj = jj + 1
1448 wa(jj) = sig(6)
1449c
1450 jj = jj + 1
1451 IF (gbuf%G_PLA > 0) THEN
1452 wa(jj) = gbuf%PLA(i)
1453 ELSE
1454 wa(jj) = zero
1455 ENDIF
1456! Upper
1457
1458 a1 = one
1459 a2 = six
1460
1461 sig(1) = a1*gbuf%FOR(kk(1)+i) + a2* gbuf%MOM(kk(1)+i)
1462 sig(2) = a1*gbuf%FOR(kk(2)+i) + a2* gbuf%MOM(kk(2)+i)
1463 sig(3) = a1*gbuf%FOR(kk(3)+i) + a2* gbuf%MOM(kk(3)+i)
1464 sig(4) = gbuf%FOR(kk(4)+i)
1465 sig(5) = gbuf%FOR(kk(5)+i)
1466 sig(6) = zero
1467
1468 iorth = 0
1469 ilay = 0
1470 jdir = 0
1471
1472 CALL shell_rota(
1473 1 i ,ilay ,nel ,iorth ,ity ,
1474 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1475 3 rx ,ry ,rz ,sx ,sy ,
1476 4 sz ,e1x ,e2x ,e3x ,e1y ,
1477 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1478 6 dir_a ,dir_b )
1479
1480 jj = jj + 1
1481 wa(jj) = one
1482 jj = jj + 1
1483 wa(jj) = sig(1)
1484 jj = jj + 1
1485 wa(jj) = sig(2)
1486 jj = jj + 1
1487 wa(jj) = sig(3)
1488 jj = jj + 1
1489 wa(jj) = sig(4)
1490 jj = jj + 1
1491 wa(jj) = sig(5)
1492 jj = jj + 1
1493 wa(jj) = sig(6)
1494c
1495 jj = jj + 1
1496 IF (gbuf%G_PLA > 0) THEN
1497 wa(jj) = gbuf%PLA(i)
1498 ELSE
1499 wa(jj) = zero
1500 ENDIF
1501 ELSE
1502! LOWER
1503 a1 = one
1504 a2 = -six
1505
1506 DO ir=1,nptr
1507 DO is=1,npts
1508 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
1509 ipg = nptr*(is-1) + ir
1510 k = (ipg-1)*nel*5
1511 kb = (ipg-1)*nel*3
1512
1513 sig(1) = a1*gbuf%FORPG(k + kk(1) + i) + a2*gbuf%MOMPG(kb + kk(1) + i)
1514 sig(2) = a1*gbuf%FORPG(k + kk(2) + i) + a2*gbuf%MOMPG(kb + kk(2) + i)
1515 sig(3) = a1*gbuf%FORPG(k + kk(3) + i) + a2*gbuf%MOMPG(kb + kk(3) + i)
1516 sig(4) = gbuf%FORPG(k + kk(4) + i)
1517 sig(5) = gbuf%FORPG(k + kk(5) + i)
1518 sig(6) = zero
1519
1520 iorth = 0
1521 ilay = 0
1522 jdir = 0
1523
1524 CALL shell_rota(
1525 1 i ,ilay ,nel ,iorth ,ity ,
1526 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1527 3 rx ,ry ,rz ,sx ,sy ,
1528 4 sz ,e1x ,e2x ,e3x ,e1y ,
1529 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1530 6 dir_a ,dir_b )
1531
1532 jj = jj + 1
1533 wa(jj) = -one
1534 jj = jj + 1
1535 wa(jj) = sig(1)
1536 jj = jj + 1
1537 wa(jj) = sig(2)
1538 jj = jj + 1
1539 wa(jj) = sig(3)
1540 jj = jj + 1
1541 wa(jj) = sig(4)
1542 jj = jj + 1
1543 wa(jj) = sig(5)
1544 jj = jj + 1
1545 wa(jj) = sig(6)
1546c
1547 jj = jj + 1
1548 IF (gbuf%G_PLA > 0) THEN
1549 wa(jj) = lbuf%PLA(i)
1550 ELSE
1551 wa(jj) = zero
1552 ENDIF
1553 ENDDO
1554 ENDDO
1555
1556! MEMBRANE
1557 a1 = one
1558 a2 = zero
1559
1560 DO ir=1,nptr
1561 DO is=1,npts
1562 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
1563 ipg = nptr*(is-1) + ir
1564 k = (ipg-1)*nel*5
1565 kb = (ipg-1)*nel*3
1566 sig(1) = a1*gbuf%FORPG(k + kk(1) + i) + a2*gbuf%MOMPG(kb + kk(1) + i)
1567 sig(2) = a1*gbuf%FORPG(k + kk(2) + i) + a2*gbuf%MOMPG(kb + kk(2) + i)
1568 sig(3) = a1*gbuf%FORPG(k + kk(3) + i) + a2*gbuf%MOMPG(kb + kk(3) + i)
1569 sig(4) = gbuf%FORPG(k + kk(4) + i)
1570 sig(5) = gbuf%FORPG(k + kk(5) + i)
1571 sig(6) = zero
1572
1573 iorth = 0
1574 ilay = 0
1575 jdir = 0
1576
1577 CALL shell_rota(
1578 1 i ,ilay ,nel ,iorth ,ity ,
1579 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1580 3 rx ,ry ,rz ,sx ,sy ,
1581 4 sz ,e1x ,e2x ,e3x ,e1y ,
1582 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1583 6 dir_a ,dir_b )
1584
1585 jj = jj + 1
1586 wa(jj) = zero
1587 jj = jj + 1
1588 wa(jj) = sig(1)
1589 jj = jj + 1
1590 wa(jj) = sig(2)
1591 jj = jj + 1
1592 wa(jj) = sig(3)
1593 jj = jj + 1
1594 wa(jj) = sig(4)
1595 jj = jj + 1
1596 wa(jj) = sig(5)
1597 jj = jj + 1
1598 wa(jj) = sig(6)
1599c
1600 jj = jj + 1
1601 IF (gbuf%G_PLA > 0) THEN
1602 wa(jj) = lbuf%PLA(i)
1603 ELSE
1604 wa(jj) = zero
1605 ENDIF
1606 ENDDO
1607 ENDDO
1608
1609! Upper
1610 a1 = one
1611 a2 = six
1612 DO ir=1,nptr
1613 DO is=1,npts
1614 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,1)
1615 ipg = nptr*(is-1) + ir
1616 k = (ipg-1)*nel*5
1617 kb = (ipg-1)*nel*3
1618 sig(1) = a1*gbuf%FORPG(k + kk(1) + i) + a2*gbuf%MOMPG(kb + kk(1) + i)
1619 sig(2) = a1*gbuf%FORPG(k + kk(2) + i) + a2*gbuf%MOMPG(kb + kk(2) + i)
1620 sig(3) = a1*gbuf%FORPG(k + kk(3) + i) + a2*gbuf%MOMPG(kb + kk(3) + i)
1621 sig(4) = gbuf%FORPG(k + kk(4) + i)
1622 sig(5) = gbuf%FORPG(k + kk(5) + i)
1623 sig(6) = zero
1624
1625 iorth = 0
1626 ilay = 0
1627 jdir = 0
1628
1629 CALL shell_rota(
1630 1 i ,ilay ,nel ,iorth ,ity ,
1631 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1632 3 rx ,ry ,rz ,sx ,sy ,
1633 4 sz ,e1x ,e2x ,e3x ,e1y ,
1634 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1635 6 dir_a ,dir_b )
1636
1637 jj = jj + 1
1638 wa(jj) = one
1639 jj = jj + 1
1640 wa(jj) = sig(1)
1641 jj = jj + 1
1642 wa(jj) = sig(2)
1643 jj = jj + 1
1644 wa(jj) = sig(3)
1645 jj = jj + 1
1646 wa(jj) = sig(4)
1647 jj = jj + 1
1648 wa(jj) = sig(5)
1649 jj = jj + 1
1650 wa(jj) = sig(6)
1651c
1652 jj = jj + 1
1653 IF (gbuf%G_PLA > 0) THEN
1654 wa(jj) = lbuf%PLA(i)
1655 ELSE
1656 wa(jj) = zero
1657 ENDIF
1658 ENDDO
1659 ENDDO
1660
1661 ENDIF ! IF (MLW == 0 .or. MLW == 13)
1662 ELSE ! MPT > 0
1663 IF (mlw == 0 .or. mlw == 13) THEN
1664 DO ipg=1,npg
1665 DO j=1,8
1666 jj = jj + 1
1667 wa(jj) = zero
1668 ENDDO
1669 ENDDO
1670 ELSEIF (nlay == 1) THEN ! PID1
1671
1672 bufly => elbuf_tab(ng)%BUFLY(1)
1673 nptt = bufly%NPTT
1674
1675 DO it=1,nptt
1676 DO is=1,npts
1677 DO ir=1,nptr
1678 lbuf => bufly%LBUF(ir,is,it)
1679 ipg = nptr*(is-1) + ir
1680 sig(1) = lbuf%SIG(kk(1)+i)
1681 sig(2) = lbuf%SIG(kk(2)+i)
1682 sig(3) = lbuf%SIG(kk(3)+i)
1683 sig(4) = lbuf%SIG(kk(4)+i)
1684 sig(5) = lbuf%SIG(kk(5)+i)
1685 sig(6) = zero
1686C
1687C viscous stress added for animation
1688C
1689 IF (ivisc > 0) THEN
1690
1691 sig(1) = sig(1) + lbuf%VISC(kk(1)+i)
1692 sig(2) = sig(2) + lbuf%VISC(kk(2)+i)
1693 sig(3) = sig(3) + lbuf%VISC(kk(3)+i)
1694 sig(4) = sig(4) + lbuf%VISC(kk(4)+i)
1695 sig(5) = sig(5) + lbuf%VISC(kk(5)+i)
1696
1697 ENDIF ! IF (IVISC > 0)
1698
1699 ilay = 1
1700 jdir = 1 + (ilay-1)*llt*2
1701
1702 CALL shell_rota(
1703 1 i ,ilay ,nel ,iorth ,ity ,
1704 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1705 3 rx ,ry ,rz ,sx ,sy ,
1706 4 sz ,e1x ,e2x ,e3x ,e1y ,
1707 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1708 6 dir_a ,dir_b )
1709
1710 jj = jj + 1
1711 wa(jj) = two * posly(i,it)
1712 jj = jj + 1
1713 wa(jj) = sig(1)
1714 jj = jj + 1
1715 wa(jj) = sig(2)
1716 jj = jj + 1
1717 wa(jj) = sig(3)
1718 jj = jj + 1
1719 wa(jj) = sig(4)
1720 jj = jj + 1
1721 wa(jj) = sig(5)
1722 jj = jj + 1
1723 wa(jj) = sig(6)
1724 jj = jj + 1
1725 IF (bufly%L_PLA > 0) THEN
1726 wa(jj) = lbuf%PLA(i)
1727 ELSE
1728 wa(jj) = zero
1729 ENDIF
1730 ENDDO
1731 ENDDO
1732 ENDDO ! DO IPT = 1,NPTT
1733
1734 ELSE
1735 iptt = 0
1736 DO il=1,nlay
1737 bufly => elbuf_tab(ng)%BUFLY(il)
1738 nptt = bufly%NPTT
1739 DO it=1,nptt
1740 DO ipg=1,npg
1741 lbuf => bufly%LBUF(ipg,1,it)
1742 l_pla = bufly%L_PLA
1743
1744 sig(1) = lbuf%SIG(kk(1)+i)
1745 sig(2) = lbuf%SIG(kk(2)+i)
1746 sig(3) = lbuf%SIG(kk(3)+i)
1747 sig(4) = lbuf%SIG(kk(4)+i)
1748 sig(5) = lbuf%SIG(kk(5)+i)
1749 sig(6) = zero
1750C
1751C viscous stress added for animation
1752C
1753 IF (ivisc > 0) THEN
1754
1755 sig(1) = sig(1) + lbuf%VISC(kk(1)+i)
1756 sig(2) = sig(2) + lbuf%VISC(kk(2)+i)
1757 sig(3) = sig(3) + lbuf%VISC(kk(3)+i)
1758 sig(4) = sig(4) + lbuf%VISC(kk(4)+i)
1759 sig(5) = sig(5) + lbuf%VISC(kk(5)+i)
1760
1761 ENDIF ! IF (IVISC > 0)
1762
1763 ilay = 0
1764 jdir = 1 + (il-1)*llt*2
1765
1766 CALL shell_rota(
1767 1 i ,ilay ,nel ,iorth ,ity ,
1768 2 igtyp ,mlw ,jdir ,sig ,elbuf_tab(ng),
1769 3 rx ,ry ,rz ,sx ,sy ,
1770 4 sz ,e1x ,e2x ,e3x ,e1y ,
1771 5 e2y ,e3y ,e1z ,e2z ,e3z ,
1772 6 dir_a ,dir_b )
1773 jj = jj + 1
1774 wa(jj) = two * posly(i,iptt)
1775 jj = jj + 1
1776 wa(jj) = sig(1)
1777 jj = jj + 1
1778 wa(jj) = sig(2)
1779 jj = jj + 1
1780 wa(jj) = sig(3)
1781 jj = jj + 1
1782 wa(jj) = sig(4)
1783 jj = jj + 1
1784 wa(jj) = sig(5)
1785 jj = jj + 1
1786 wa(jj) = sig(6)
1787 jj = jj + 1
1788 IF (l_pla > 0) THEN
1789 wa(jj) = lbuf%PLA(i)
1790 ELSE
1791 wa(jj) = zero
1792 ENDIF
1793 ENDDO ! DO IPG=1,NPG
1794 ENDDO ! DO IT=1,NPTT
1795 iptt = iptt + nptt
1796 ENDDO ! DO IL=1,NLAY
1797 ENDIF ! IF (MLW == 0 .or. MLW == 13)
1798 ENDIF ! IF (MPT == 0)
1799C
1800 ie=ie+1
1801C pointer for last position for element IE
1802 ptwa(ie)=jj
1803 ENDDO ! DO I=LFT,LLT
1804
1805 IF (ALLOCATED(dirb)) DEALLOCATE(dirb)
1806 IF (ALLOCATED(dira)) DEALLOCATE(dira)
1807 DEALLOCATE(matly, thkly, posly, thk_ly)
1808 ENDIF ! IF (ITY == 7)
1809 ENDDO ! DO NG=1,NGROUP
1810 ENDIF
1811C
1812c-----------------------------------------------------------------------
1813 IF (nspmd == 1) THEN
1814C recopying for code simplification
1815 len=jj
1816 DO j=1,len
1817 wap0(j)=wa(j)
1818 ENDDO
1819 ptwa_p0(0)=0
1820 DO n=1,dynain_data%DYNAIN_NUMELTG
1821 ptwa_p0(n)=ptwa(n)
1822 ENDDO
1823 ELSE
1824C Global pointers WAP0
1825 CALL spmd_stat_pgather(ptwa,dynain_data%DYNAIN_NUMELTG,ptwa_p0,dynain_data%DYNAIN_NUMELTG_G)
1826 len = 0
1827 CALL spmd_rgather9_dp(wa,jj,wap0,sizp0,len)
1828 ENDIF
1829
1830 IF (ispmd == 0.AND.len > 0) THEN
1831 IF(is_written == 0 ) THEN
1832 IF(dynain_data%ZIPDYNAIN==0) THEN
1833 WRITE(iudynain,'(A)') delimit
1834 WRITE(iudynain,'(A)')'*INITIAL_STRESS_SHELL'
1835 WRITE(iudynain,'(A)')
1836 . '$ SHELLID NPG NBINT LARGE '
1837 WRITE(iudynain,'(A)')
1838 . '$ IF(NPT == 0), REPEAT I=1,NPG :'
1839 WRITE(iudynain,'(A)')
1840 . '$ IF(NPT /= 0) REPEAT K=1,NPT : REPEAT I=1,NPG :'
1841 WRITE(iudynain,'(A)')
1842 . '$ T SIGXX SIGYY SIGZZ SIGXY SIGYZ SIGZX EPSP '
1843 WRITE(iudynain,'(A)') delimit
1844 ELSE
1845 WRITE(line,'(A)')'*INITIAL_STRESS_SHELL'
1846 CALL strs_txt50(line,100)
1847 WRITE(line,'(A)')
1848 . '$ SHELLID NPG NBINT LARGE '
1849 CALL strs_txt50(line,100)
1850 WRITE(line,'(A)')
1851 . '$ IF(NPT == 0), REPEAT I=1,NPG :'
1852 CALL strs_txt50(line,100)
1853 WRITE(line,'(A)')
1854 . '$ IF(NPT /= 0) REPEAT K=1,NPT : REPEAT I=1,NPG :'
1855 CALL strs_txt50(line,100)
1856 WRITE(line,'(A)')
1857 . '$ T SIGXX SIGYY SIGZZ SIGXY SIGYZ SIGZX EPSP '
1858 CALL strs_txt50(line,100)
1859 WRITE(line,'(A)') delimit
1860 CALL strs_txt50(line,100)
1861 ENDIF
1862 is_written = 1
1863 ENDIF
1864
1865 DO n=1,dynain_data%DYNAIN_NUMELTG_G
1866C Retrieving shell ID in increasing order
1867 k=dynain_indxtg(n)
1868C Adress in WAP0
1869 j=ptwa_p0(k-1)
1870C
1871 ioff = nint(wap0(j + 1))
1872 IF (ioff >= 1) THEN
1873 id = nint(wap0(j + 2))
1874 npt = nint(wap0(j + 3))
1875 npg = nint(wap0(j + 4))
1876 large = nint(wap0(j + 5))
1877 j = j + 5
1878 IF(dynain_data%ZIPDYNAIN==0) THEN
1879 WRITE(iudynain,'(3I8,16X,I8)')id,npg,npt,large
1880 ELSE
1881 WRITE(line,'(3I8,16X,I8)')id,npg,npt,large
1882 CALL strs_txt50(line,100)
1883 ENDIF
1884 IF (npt == 0) THEN
1885 DO ipg=1,npg
1886 IF(dynain_data%ZIPDYNAIN==0) THEN
1887 WRITE(iudynain,'(1P5G16.9)')(wap0(jj + k),k=1,5)
1888 WRITE(iudynain,'(1P3G16.9)')(wap0(jj + k),k=6,8)
1889 ELSE
1890 WRITE(line,'(1P5G16.9)')(wap0(jj + k),k=1,5)
1891 CALL strs_txt50(line,100)
1892 WRITE(line,'(1P3G16.9)')(wap0(jj + k),k=6,8)
1893 CALL strs_txt50(line,100)
1894 ENDIF
1895 j = j + 8
1896 ENDDO
1897 ELSE
1898 DO ipt=1,npt
1899 DO ipg=1,npg
1900 IF(dynain_data%ZIPDYNAIN==0) THEN
1901 WRITE(iudynain,'(1P5G16.9)')(wap0(j + k),k=1,5)
1902 WRITE(iudynain,'(1P3G16.9)')(wap0(j + k),k=6,8)
1903 ELSE
1904 WRITE(line,'(1P5G16.9)')(wap0(j + k),k=1,5)
1905 CALL strs_txt50(line,100)
1906 WRITE(line,'(1P3G16.9)')(wap0(j + k),k=6,8)
1907 CALL strs_txt50(line,100)
1908 ENDIF
1909 j = j + 8
1910 ENDDO
1911 ENDDO
1912
1913 ENDIF ! IF (NPT == 0)
1914 ENDIF ! IF (IOFF >= 1)
1915 ENDDO ! DO N=1,DYNAIN_NUMELTG_G
1916 ENDIF ! IF (ISPMD == 0.AND.LEN > 0)
1917C
1918 DEALLOCATE(ptwa,ptwa_p0)
1919C
1920 RETURN
1921 END
cortdir3
subroutine cortdir3(elbuf_str, dir_a, dir_b, jft, jlt, nlay, irep, rx, ry, rz, sx, sy, sz, e1x, e1y, e1z, e2x, e2y, e2z, nel)
Definition cortdir3.F:45
my_real
#define my_real
Definition cppsort.cpp:32
dynain_c_strsg
subroutine dynain_c_strsg(elbuf_tab, iparg, igeo, ixc, ixtg, wa, wap0, ipartc, iparttg, dynain_data, dynain_indxc, dynain_indxtg, sizp0, geo, stack, drape_sh4n, drape_sh3n, x, thke, drapeg, nummat, mat_param)
Definition dynain_c_strsg.F:49
layini
subroutine layini(elbuf_str, jft, jlt, geo, igeo, mat, pid, thkly, matly, posly, igtyp, ixfem, ixlay, nlay, npt, isubstack, stack, drape, nft, thk, nel, ratio_thkly, indx_drape, sedrape, numel_drape)
Definition layini.F:47
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
drape_mod
Definition drape_mod.F:54
drape_mod::scdrape
integer scdrape
Definition drape_mod.F:92
drape_mod::stdrape
integer stdrape
Definition drape_mod.F:92
stack_mod
Definition stack_mod.F:75
state_mod
Definition state_mod.F:54
shell_local_frame
subroutine shell_local_frame(jft, jlt, ity, ihbe, igtyp, ixc, ixtg, nft, x, offg, rx, ry, rz, sx, sy, sz, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z)
Definition shell_rota.F:227
shell_rota
subroutine shell_rota(i, ilay, nel, iorth, ity, igtyp, mlw, jdir, sig, elbuf_str, rx, ry, rz, sx, sy, sz, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, dir_a, dir_b)
Definition shell_rota.F:36
spmd_rgather9_dp
subroutine spmd_rgather9_dp(v, len, vp0, lenp0, iad)
Definition spmd_outp.F:1019
spmd_stat_pgather
subroutine spmd_stat_pgather(ptv, ptlen, ptv_p0, ptlen_p0)
Definition spmd_stat.F:53
strs_txt50
subroutine strs_txt50(text, length)
Definition sta_txt.F:87
state_mod::dynain_database
Definition state_mod.F:64