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