38 SUBROUTINE c3inmas(X,XREFTG ,IXTG ,GEO,PM ,MS ,TINER,THKE ,
39 . PARTSAV,V ,IPART ,MSTG ,INTG ,
40 . PTG ,IGEO ,IMAT ,IPROP ,AREA ,
41 . ETNOD ,NSHNOD ,STTG ,SH3TREE ,MCP ,
42 . MCPTG , TEMP ,SH3TRIM,ISUBSTACK,NLAY ,
43 . ELBUF_STR, STACK,THKI ,RNOISE ,DRAPE ,
44 . PERTURB,IX1 ,IX2 ,IX3 ,NINTEMP,
45 . X2 ,X3 ,Y3 ,IDRAPE ,INDX)
54 use element_mod ,
only : nixtg
60#include "implicit_f.inc"
70#include "remesh_c.inc"
74#include "vect01_c.inc"
79 INTEGER IXTG(NIXTG,*),IPART(*),IMAT,IPROP,
80 . IGEO(NPROPGI,*),IGTYP, NSHNOD(*), SH3TREE(KSH3TREE,*),
81 . SH3TRIM(*),ISUBSTACK,NLAY,IDRAPE,
83 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ)
84 INTEGER,
DIMENSION (STDRAPE) :: INDX
85 INTEGER ,
INTENT(IN) :: NINTEMP
87 . x(3,*), geo(npropg,*), pm(npropm,*), ms(*), tiner(*),thke(*),
88 . v(3,*),partsav(20,*),mstg(*),intg(*),ptg(3,*),
89 . etnod(*), sttg(*),mcp(*),mcptg(*),temp(*),xreftg(3,3,*),thki(*),
90 . rnoise(nperturb,*),
area(mvsiz)
91 TYPE(elbuf_struct_),
TARGET :: ELBUF_STR
92 TYPE (STACK_PLY) :: STACK
93 TYPE (DRAPE_),
DIMENSION(NUMELC_DRAPE + NUMELTG_DRAPE),
TARGET :: DRAPE
97 INTEGER I,J,K,N,IP,I1,I2,I3,I4,IPTHK,IPMAT,IPPOS,MATLY,IDPROP,
98 . ipid,ippid,iigeo,iadi,iadr,ipang,nthk,igmat,ilay,nptt,
99 . it,max_nptt,iint,ipid_ly,ipgmat,ipos,ipt,ipt_all,nslice,ie,
102 CHARACTER(LEN=NCHARTITLE) :: TITR
103 my_real XX,YY,ZZ,XY,YZ,ZX,ZSHIFT
104 my_real x2(mvsiz), x3(mvsiz), y3(mvsiz),
105 . rho(mvsiz), thk(mvsiz), em(mvsiz), xi(mvsiz),
106 . p1, p2, p3, aa, bb, cc, a2, b2, c2, xin, ems, msl, xil,
107 . mst,dms,ins,posly,thkly,thickt, et,emscp(mvsiz),rhocp(mvsiz),
108 . ems_nptt,ms_nptt,pos_0,thinning,scale,thick_drape
109 my_real,
DIMENSION(:) ,
ALLOCATABLE :: thk_it, pos_it
110 my_real,
DIMENSION(:,:),
ALLOCATABLE :: laypos, ratio_thkly
112 TYPE (DRAPE_PLY_) ,
POINTER :: DRAPE_PLY
115 . A_GAUSS(9,9),W_GAUSS(9,9)
121 2 -.577350269189626,0.577350269189626,0. ,
124 3 -.774596669241483,0. ,0.774596669241483,
127 4 -.861136311594053,-.339981043584856,0.339981043584856,
128 4 0.861136311594053,0. ,0. ,
130 5 -.906179845938664,-.538469310105683,0. ,
131 5 0.538469310105683,0.906179845938664,0. ,
133 6 -.932469514203152,-.661209386466265,-.238619186083197,
134 6 0.238619186083197,0.661209386466265,0.932469514203152,
136 7 -.949107912342759,-.741531185599394,-.405845151377397,
137 7 0. ,0.405845151377397,0.741531185599394,
138 7 0.949107912342759,0. ,0. ,
139 8 -.960289856497536,-.796666477413627,-.525532409916329,
140 8 -.183434642495650,0.183434642495650,0.525532409916329,
141 8 0.796666477413627,0.960289856497536,0. ,
142 9 -.968160239507626,-.836031107326636,-.613371432700590,
143 9 -.324253423403809,0. ,0.324253423403809,
144 9 0.613371432700590,0.836031107326636,0.968160239507626/
152 3 0.555555555555556,0.888888888888889,0.555555555555556,
155 4 0.347854845137454,0.652145154862546,0.652145154862546,
156 4 0.347854845137454,0. ,0. ,
158 5 0.236926885056189,0.478628670499366,0.568888888888889,
159 5 0.478628670499366,0.236926885056189,0. ,
161 6 0.171324492379170,0.360761573048139,0.467913934572691,
162 6 0.467913934572691,0.360761573048139,0.171324492379170,
164 7 0.129484966168870,0.279705391489277,0.381830050505119,
165 7 0.417959183673469,0.381830050505119,0.279705391489277,
166 7 0.129484966168870,0. ,0. ,
167 8 0.101228536290376,0.222381034453374,0.313706645877887,
168 8 0.362683783378362,0.362683783378362,0.313706645877887,
169 8 0.222381034453374,0.101228536290376,0. ,
170 9 0.081274388361574,0.180648160694857,0.260610696402935,
171 9 0.312347077040003,0.330239355001260,0.312347077040003,
172 9 0.260610696402935,0.180648160694857,0.081274388361574/
174 igtyp = nint(geo(12,iprop))
175 igmat = igeo(98,iprop)
177 IF(igtyp == 51 .OR. igtyp == 52)
THEN
179 max_nptt =
max(max_nptt ,elbuf_str%BUFLY
181 ALLOCATE(thk_it(max_nptt),pos_it(max_nptt))
183 ALLOCATE(thk_it(0),pos_it(0))
185 nlay_max =
max(nlay, npt)
186 ALLOCATE( laypos(nlay_max*max_nptt,mvsiz),ratio_thkly(nlay_max*max_nptt,mvsiz) )
187 IF ((igtyp == 17 .OR. igtyp == 51) .AND. igmat < 0)
THEN
189 IF (ndrape == 0)
THEN
190 IF (thke(i)== zero)
THEN
191 thk(i) = stack%GEO(1,isubstack)
196 thki(i) = stack%GEO(1,isubstack)
197 ELSEIF (ndrape > 0)
THEN
202 IF (nperturb /= 0)
THEN
204 IF (perturb(j) == 1 .AND.
205 . rnoise(j,numelc+i+nft) /= zero)
THEN
206 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
207 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
208 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
214 rhocp(i) = pm(69,imat)
215 IF(thk(i)==zero.AND.igtyp/=0)
THEN
218 . anmode=aninfo_blind_1,
222 ELSEIF (igtyp == 52 .OR.
223 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))
THEN
225 IF (ndrape == 0)
THEN
226 IF (thke(i)== zero)
THEN
227 thk(i) = stack%GEO(1,isubstack)
232 thki(i) = stack%GEO(1,isubstack)
233 ELSEIF (ndrape > 0)
THEN
238 IF (nperturb /= 0)
THEN
240 IF (perturb(j) == 1 .AND.
241 . rnoise(j,numelc+i+nft) /= zero)
THEN
242 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
243 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
244 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
249 rho(i) = stack%PM(1,isubstack)
250 rhocp(i) = stack%PM(15,isubstack)* rho(i)/stack%PM(14,isubstack)
251 IF (thk(i) == zero .AND. igtyp /= 0)
THEN
254 . anmode=aninfo_blind_1,
260 IF(thke(i)==zero)
THEN
268 IF (nperturb /= 0)
THEN
270 IF (perturb(j) == 1 .AND.
271 . rnoise(j,numelc+i+nft) /= zero)
THEN
272 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
273 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
274 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
280 rhocp(i) = pm(69,imat)
281 IF(thk(i)==zero.AND.igtyp/=0)
THEN
284 . anmode=aninfo_blind_1,
292 IF (igtyp==11 .OR. igtyp==16)
THEN
301 thickt = geo(200,iprop)
302 thkly = geo(i1,iprop)*thk(i)
303 matly = igeo(i2,iprop)
304 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
307 mst = rho(i)*thk(i)*
area(i)
308 dms = abs(em(i)-mst)/mst
314 idprop = nint(geo(40,iprop))
316 . igeo(npropgi-ltitr+1,iprop),ltitr)
318 . msgtype=msgwarning,
319 . anmode=aninfo_blind_2,
321 . c1=titr,i2=ixtg(nixtg,nft+i),
326 ELSEIF (igtyp == 17)
THEN
341 thkly = stack%GEO(i1,isubstack)*thk(i)
342 matly = stack%IGEO(i3,isubstack)
343 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
346 mst = rho(i)*thk(i)*
area(i)
347 dms = abs(em(i)-mst)/mst
350 idprop = nint(geo(40,iprop))
352 . igeo(npropgi-ltitr+1,iprop),ltitr)
354 . msgtype=msgwarning,
355 . anmode=aninfo_blind_2,
357 . c1=titr,i2=ixtg(nixtg,nft+i),
371 thickt = stack%GEO(1,isubstack)
372 thkly = stack%GEO(i1,isubstack)*thk(i)
373 matly = stack%IGEO(i3,isubstack)
374 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
377 thick_drape = drape(ie)%THICK
378 scale = thk(i)/thick_drape
383 ip = drape(ie)%INDX_PLY(n)
385 drape_ply => drape(ie)%DRAPE_PLY(ip)
386 nslice = drape_ply%NSLICE
387 thinning = drape_ply%RDRAPE(1,1)
388 thickt = stack%GEO(1,isubstack)
389 thkly = stack%GEO(i1,isubstack)*thickt
390 matly = stack%IGEO(i3
391 thkly = thkly*thinning*scale
392 matly = stack%IGEO(i3,isubstack)
393 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
395 thickt = stack%GEO(1,isubstack)
396 thkly = stack%GEO(i1,isubstack)*thickt*scale
397 matly = stack%IGEO(i3,isubstack)
398 matly = stack%IGEO(i3,isubstack)
399 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
404 mst = rho(i)*thk(i)*
area(i)
405 dms = abs(em(i)-mst)/mst
408 idprop = nint(geo(40,iprop))
410 . igeo(npropgi-ltitr+1,iprop),ltitr)
412 . msgtype=msgwarning,
413 . anmode=aninfo_blind_2,
415 . c1=titr,i2=ixtg(nixtg,nft+i),
421 ELSEIF (igtyp == 51 .OR. igtyp == 52)
THEN
429 ipos = igeo(99,iprop)
430 thickt = stack%GEO(1,isubstack)
431 zshift = geo(199,iprop)
432 IF(ipos == 2 ) zshift = - zshift /
max(thickt,em20)
435 scale = thk(i)/thickt
439 nptt = elbuf_str%BUFLY(ilay)%NPTT
443 thickt = stack%GEO(1,isubstack)
444 thkly = stack%GEO(i1,isubstack)*thickt
445 ipid_ly= stack%IGEO(i2,isubstack)
446 matly = stack%IGEO(i3,isubstack)
447 ipid = stack%IGEO(ippid,isubstack) ! or ipid = ix(nixc-1,nft+i)
448 iint = igeo(47,iprop)
451 thk_it(it) = thkly/nptt
453 ELSEIF (iint == 2)
THEN
455 thk_it(it) = half*thkly*w_gauss(it,nptt)
460 thk_it(it) = thk_it(it)*scale
461 ems_nptt = pm(89,matly)*thk_it(it)*
area
462 em(i) = em(i) + ems_nptt
467 mst = rho(i)*thk(i)*
area(i)
468 dms = abs(em(i)-mst)/mst
469 IF(igtyp == 51 .AND.igmat < 0)
THEN
471 idprop = nint(geo(40,iprop))
473 . igeo(npropgi-ltitr+1,iprop),ltitr)
475 . msgtype=msgwarning,
476 . anmode=aninfo_blind_2,
478 . c1=titr,i2=ixtg(nixtg,nft+i),
484 thickt = stack%GEO(1,isubstack)
489 scale = thk(i)/thickt
492 nptt = elbuf_str%BUFLY(ilay)%NPTT
496 ipid = stack%IGEO(i2,isubstack)
497 iint = igeo(47,iprop)
498 thickt = stack%GEO(1,isubstack)
500 thickt = stack%GEO(1,isubstack)
501 thkly = stack%GEO(i1,isubstack)*thickt
502 ipid_ly= stack%IGEO(i2,isubstack)
503 matly = stack%IGEO(i3,isubstack)
505 thk_it(it) = thkly/nptt
507 ELSEIF(iint == 2)
THEN
508 thickt = stack%GEO(1,isubstack)
509 thkly = stack%GEO(i1,isubstack)*thickt
510 ipid_ly= stack%IGEO(i2,isubstack)
511 matly = stack%IGEO(i3,isubstack)
513 thk_it(it) = half*thkly*w_gauss(it,nptt)
518 thk_it(it) = scale*thk_it(it)
519 ems_nptt = pm(89,matly)*thk_it(it)*
area
520 em(i) = em(i) + ems_nptt
525 thick_drape = drape(ie)%THICK
526 scale = thk(i)/thick_drape
528 nptt = elbuf_str%BUFLY(ilay)%NPTT
532 ipid = stack%IGEO(i2,isubstack)
533 iint = igeo(47,iprop)
534 thickt = stack%GEO(1,isubstack)
535 ip = drape(ie)%INDX_PLY(ilay)
537 thickt = stack%GEO(1,isubstack)
538 thkly = stack%GEO(i1,isubstack)*thickt
539 ipid_ly= stack%IGEO(i2,isubstack)
540 matly = stack%IGEO(i3,isubstack)
543 thk_it(it) = thkly/nptt
545 ELSEIF(iint == 2)
THEN
547 thk_it(it) = half*thkly*w_gauss(it,nptt)
551 drape_ply => drape(ie)%DRAPE_PLY(ip)
552 nslice = drape_ply%NSLICE
553 thickt = stack%GEO(1,isubstack)
554 thkly = stack%GEO(i1,isubstack)*thickt
555 ipid_ly= stack%IGEO(i2,isubstack)
556 matly = stack%IGEO(i3,isubstack)
559 thinning = drape_ply%RDRAPE(it,1)
560 thk_it(it) = thkly*thinning/nptt
562 ELSEIF(iint == 2)
THEN
564 thinning = drape_ply%RDRAPE(it,1)
565 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
570 thk_it(it) = scale*thk_it(it)
571 ems_nptt = pm(89,matly)*thk_it(it)*
area(i)
572 em(i) = em(i) + ems_nptt
577 mst = rho(i)*thk(i)*
area(i)
578 dms = abs(em(i)-mst)/mst
579 IF(igtyp == 51 .AND.igmat < 0)
THEN
581 idprop = nint(geo(40,iprop))
583 . igeo(npropgi-ltitr+1,iprop),ltitr)
585 . msgtype=msgwarning,
586 . anmode=aninfo_blind_2,
588 . c1=titr,i2=ixtg(nixtg,nft+i),
596 em(i) = rho(i)*thk(i)*
area(i)
599 emscp(i)=rhocp(i)*thk(i)*
area(i)
605 IF (igtyp==11 .OR. igtyp==16)
THEN
612 thickt= geo(200,iprop)
613 thkly = geo(i1,iprop)*thk(i)
614 posly = geo(i3,iprop)*thk(i)
615 matly = igeo(i2,iprop)
616 msl = pm(89,matly)*thkly*
area(i)
617 xil = msl*(
area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
622 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
623 dms = abs(xi(i)-ins)/ins
625 idprop = nint(geo(40,iprop))
627 . igeo(npropgi-ltitr+1,iprop),ltitr)
629 . msgtype=msgwarning,
630 . anmode=aninfo_blind_2,
632 . c1=titr,i2=ixtg(nixtg,nft+i))
636 ELSEIF (igtyp == 17)
THEN
637 ipos = igeo(99,iprop)
638 thickt = stack%GEO(1,isubstack)
639 zshift = geo(199,iprop)
640 IF(ipos == 2 ) zshift = - zshift /
max(thickt,em20)
641 IF(idrape == 0 )
THEN
649 thickt= stack%GEO(1,isubstack)
650 thkly = stack%GEO(i3,isubstack)*thk(i)
651 posly = stack%GEO(i4,isubstack)*thk(i)
652 matly = stack%IGEO(i2,isubstack)
653 msl = pm(89,matly)*thkly*
area(i)
654 xil = msl*(
area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
658 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
659 dms = abs(xi(i)-ins)/ins
663 idprop = nint(geo(40,iprop))
665 . igeo(npropgi-ltitr+1,iprop),ltitr)
667 . msgtype=msgwarning,
668 . anmode=aninfo_blind_2,
670 . c1=titr,i2=ixtg(nixtg,nft+i))
684 thickt= stack%GEO(1,isubstack)
685 thkly = stack%GEO(i3,isubstack)*thk(i)
686 posly = stack%GEO(i4,isubstack)*thk(i)
687 matly = stack%IGEO(i2,isubstack)
688 msl = pm(89,matly)*thkly*
area(i)
689 xil = msl*(
area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
693 thick_drape = drape(ie)%THICK
694 scale = thk(i)/thick_drape
696 ip = drape(ie)%INDX_PLY(n)
702 thkly = stack%GEO(i3,isubstack)*thickt
703 matly = stack%IGEO(i2,isubstack)
704 ratio_thkly(n,i) = thkly/ thick_drape
706 drape_ply => drape(ie)%DRAPE_PLY(ip)
707 thinning = drape_ply%RDRAPE(1,1)
708 thkly = stack%GEO(i3,isubstack)*thickt
709 thkly = thkly*thinning
710 matly = stack%IGEO(i2,isubstack)
711 ratio_thkly(n,i) = thkly/thick_drape
714 laypos(n,i) = zshift + half*ratio_thkly(n,i)
716 laypos(n,i) = laypos(n-1,i) + half*(ratio_thkly(n,i)+ratio_thkly(n-1,i))
718 posly = laypos(n,i)*thk(i)
720 msl = pm(89,matly)*thkly*
area(i)
721 xil = msl*(
area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
726 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
727 dms = abs(xi(i)-ins)/ins
731 idprop = nint(geo(40,iprop))
733 . igeo(npropgi-ltitr+1,iprop),ltitr)
735 . msgtype=msgwarning,
736 . anmode=aninfo_blind_2,
738 . c1=titr,i2=ixtg(nixtg,nft+i))
743 ELSEIF (igtyp == 51 .OR. igtyp == 52)
THEN
752 ipos = igeo(99,iprop)
753 thickt = stack%GEO(1,isubstack)
754 zshift = geo(199,iprop)
755 IF(ipos == 2 ) zshift = - zshift /
max(thickt,em20)
757 thickt = stack%GEO(1,isubstack)
761 scale = thk(i)/thickt
763 nptt = elbuf_str%BUFLY(ilay)%NPTT
768 ipid_ly = stack%IGEO(i1,isubstack)
769 matly = stack%IGEO(i2,isubstack)
770 iint = igeo(47,iprop)
772 thickt = stack%GEO(1,isubstack)
773 thkly = stack%GEO(i3,isubstack)*thickt
776 thk_it(it) = thkly/nptt
777 ratio_thkly(ipt,i) = thk_it(it)/thickt
779 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
781 laypos(ipt,i) = laypos(ipt-1,i) +
782 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
784 pos_it(it) = laypos(ipt,i)*thk(i)
786 ELSEIF(iint == 2)
THEN
787 thickt = stack%GEO(1,isubstack)
788 thkly = stack%GEO(i3,isubstack)*thickt
791 thk_it(it) = half*thkly*w_gauss(it,nptt)
792 ratio_thkly(ipt,i) = thk_it(it)/thickt
794 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
796 laypos(ipt,i) = laypos(ipt-1,i) +
797 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
799 pos_it(it) = laypos(ipt,i)*thk(i)
804 ms_nptt = pm(89,matly)*thk_it(it)
805 thk_it(it) = thk_it(it)*scale
806 xil = ms_nptt*(
area(i)/nine_over_2+thk_it(it)*thk_it(it
807 . pos_it(it)*pos_it(it))
810 ipt_all = ipt_all + nptt
813 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
814 dms = abs(xi(i)-ins)/ins
815 IF(igtyp == 51 .AND. igmat < 0)
THEN
817 idprop = nint(geo(40,iprop))
819 . igeo(npropgi-ltitr+1,iprop),ltitr)
821 . msgtype=msgwarning,
822 . anmode=aninfo_blind_2,
824 . c1=titr,i2=ixtg(nixtg,nft+i))
829 thickt = stack%GEO(1,isubstack)
834 scale = thk(i)/thickt
837 nptt = elbuf_str%BUFLY(ilay)%NPTT
842 ipid_ly = stack%IGEO(i1,isubstack)
843 matly = stack%IGEO(i2,isubstack)
844 iint = igeo(47,iprop)
846 thickt = stack%GEO(1,isubstack)
847 thkly = stack%GEO(i3,isubstack)*thickt
850 thk_it(it) = thkly/nptt
851 ratio_thkly(ipt,i) = thk_it(it)/thickt
853 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
855 laypos(ipt,i) = laypos(ipt-1,i) +
856 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
858 pos_it(it) = laypos(ipt,i)*thk(i)
860 ELSEIF(iint == 2)
THEN
861 thickt = stack%GEO(1,isubstack)
862 thkly = stack%GEO(i3,isubstack)*thickt
865 thk_it(it) = half*thkly*w_gauss(it,nptt)
866 ratio_thkly(ipt,i) = thk_it(it)/thickt
868 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
870 laypos(ipt,i) = laypos(ipt-1,i) +
871 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
873 pos_it(it) = laypos(ipt,i)*thk(i)
877 thk_it(it) = scale*thk_it(it)
878 ms_nptt = pm(89,matly)*thk_it(it)*
area(i)
879 xil = ms_nptt*(
area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
880 . pos_it(it)*pos_it(it))
883 ipt_all = ipt_all + nptt
887 thick_drape = drape(ie)%THICK
888 scale = thk(i)/thick_drape
890 ip = drape(ie)%INDX_PLY(ilay)
891 nptt = elbuf_str%BUFLY(ilay)%NPTT
896 ipid_ly = stack%IGEO(i1,isubstack)
897 matly = stack%IGEO(i2,isubstack)
898 iint = igeo(47,iprop)
900 drape_ply => drape(ie)%DRAPE_PLY(ip)
901 nslice = drape_ply%NSLICE
903 thickt = stack%GEO(1,isubstack)
904 thkly = stack%GEO(i3,isubstack)*thickt
907 thinning = drape_ply%RDRAPE(it,1)
908 thk_it(it) = thkly*thinning/nptt
909 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
911 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
913 laypos(ipt,i) = laypos(ipt-1,i) +
914 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
916 pos_it(it) = laypos(ipt,i)*thk(i)
918 ELSEIF(iint == 2)
THEN
919 thickt = stack%GEO(1,isubstack)
920 thkly = stack%GEO(i3,isubstack)*thickt
923 thinning = drape_ply%RDRAPE(it,1)
924 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
925 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
927 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
929 laypos(ipt,i) = laypos(ipt-1,i) +
930 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
932 pos_it(it) = laypos(ipt,i)*thk(i)
937 thickt = stack%GEO(1,isubstack)
938 thkly = stack%GEO(i3,isubstack)*thickt
941 thk_it(it) = thkly/nptt
942 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
944 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
946 laypos(ipt,i) = laypos(ipt-1,i) +
947 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
949 pos_it(it) = laypos(ipt,i)*thk(i)
951 ELSEIF(iint == 2)
THEN
952 thickt = stack%GEO(1,isubstack)
953 thkly = stack%GEO(i3,isubstack)*thickt
956 thk_it(it) = half*thkly*w_gauss(it,nptt)
957 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
959 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
961 laypos(ipt,i) = laypos(ipt-1,i) +
962 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
964 pos_it(it) = laypos(ipt,i)*thk(i)
969 thk_it(it) = scale*thk_it(it)
970 ms_nptt = pm(89,matly)*thk_it(it)*
area(i)
971 xil = ms_nptt*(
area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
972 . pos_it(it)*pos_it(it))
975 ipt_all = ipt_all + nptt
979 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
980 dms = abs(xi(i)-ins)/ins
981 IF(igtyp == 51 .AND. igmat < 0)
THEN
983 idprop = nint(geo(40,iprop))
985 . igeo(npropgi-ltitr+1,iprop),ltitr)
987 . msgtype=msgwarning,
988 . anmode=aninfo_blind_2,
990 . c1=titr,i2=ixtg(nixtg,nft+i))
999 xi(i)=em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
1022 b2 = (x2(i)-x3(i))**2 + y3(i)**2
1024 c2 = x3(i)**2 + y3(i)**2
1026 p1 = acos((a2 + c2 - b2)/(two * aa * cc)) / pi
1027 p2 = acos((a2 + b2 - c2)/(two * aa * bb)) / pi
1028 p3 = acos((b2 + c2 - a2)/(two * bb * cc)) / pi
1032 IF ( ( ((a2 + c2 - b2)/(two * aa * cc) <= -one) .OR.
1033 . ( (a2 + c2 - b2)/(two * aa * cc) >= one) .OR.
1034 . ( (a2 + b2 - c2)/(two * aa * bb) <= -one) .OR.
1035 . ( (a2 + b2 - c2)/(two * aa * bb) >= one) .OR.
1036 . ( (b2 + c2 - a2)/(two * bb * cc) <= -one) .OR.
1037 . ( (b2 + c2 - a2)/(two * bb * cc) >= one)) .AND.
1038 . lsh3trim /= 0 )
THEN
1039 IF (sh3trim(nft+i) /= -1)
CALL ancmsg(msgid=750,
1043 ELSEIF ( ( ((a2 + c2 - b2)/(two * aa * cc) <= -one) .OR.
1044 . ( (a2 + c2 - b2)/(two * aa * cc) >= one) .OR.
1045 . ( (a2 + b2 - c2)/(two * aa * bb) <= -one) .OR.
1046 . ( (a2 + b2 - c2)/(two * aa * bb) >= one) .OR.
1047 . ( (b2 + c2 - a2)/(two * bb * cc) <= -one) .OR.
1048 . ( (b2 + c2 - a2)/(two * bb * cc) >= one)) .AND.
1049 . lsh3trim == 0 )
THEN
1059 IF (nxref == 0)
THEN
1072 partsav(1,ip)=partsav(1,ip) + ems
1073 partsav(2,ip)=partsav(2,ip) + ems*
1074 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1075 partsav(3,ip)=partsav(3,ip) + ems*
1076 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1077 partsav(4,ip)=partsav(4,ip) + ems*
1078 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1079 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1080 . +x(1,i3)*x(1,i3)*p3)
1081 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1082 . +x(1,i3)*x(2,i3)*p3)
1083 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1084 . +x(2,i3)*x(2,i3)*p3)
1085 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1086 . +x(2,i3)*x(3,i3)*p3)
1087 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1088 . +x(3,i3)*x(3,i3)*p3)
1089 zx = (x(3,i1)*x(1,i1)*p1
1090 . +x(3,i3)*x(1,i3)*p3)
1091 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1092 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1093 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1094 partsav(8,ip) =partsav(8,ip) - ems * xy
1095 partsav(9,ip) =partsav(9,ip) - ems * yz
1096 partsav(10,ip)=partsav(10,ip) - ems * zx
1098 partsav(11,ip)=partsav(11,ip) + ems*
1099 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1100 partsav(12,ip)=partsav(12,ip) + ems*
1101 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1102 partsav(13,ip)=partsav(13,ip) + ems*
1103 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1104 partsav(14,ip)=partsav(14,ip
1105 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1106 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1107 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1122 partsav(1,ip)=partsav(1,ip) + ems
1123 partsav(2,ip)=partsav(2,ip)
1124 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)+p3*xreftg(3,1,i))
1125 partsav(3,ip)=partsav(3,ip) + ems*
1126 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1127 partsav(4,ip)=partsav(4,ip) + ems*
1128 . (p1*xreftg(1,3,i)+p2*xreftg(2,3,i)+p3*xreftg(3,3,i))
1129 xx = p1*xreftg(1,1,i)**2
1130 . + p2*xreftg(2,1,i)**2
1131 . + p3*xreftg(3,1,i)**2
1132 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1133 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1134 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1135 yy = p1*xreftg(1,2,i)**2
1136 . + p2*xreftg(2,2,i)**2
1137 . + p3*xreftg(3,2,i)**2
1138 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1139 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1140 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1141 zz = p1*xreftg(1,3,i)**2
1142 . + p2*xreftg(2,3,i)**2
1143 . + p3*xreftg(3,3,i)**2
1144 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1145 . + p2*xreftg(2,3,i)*xreftg(2,1,i)
1146 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1148 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1149 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1150 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1151 partsav(8,ip) =partsav(8,ip) - ems * xy
1152 partsav(9,ip) =partsav(9,ip) - ems * yz
1153 partsav(10,ip)=partsav(10,ip) - ems * zx
1155 partsav(11,ip)=partsav(11,ip) + ems*
1156 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1157 partsav(12,ip)=partsav(12,ip) + ems*
1158 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1159 partsav(13,ip)=partsav(13,ip) + ems*
1160 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1161 partsav(14,ip)=partsav(14,ip) + half * ems *
1162 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1163 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1164 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v
1167 ELSEIF (istatcnd==0)
THEN
1168 IF (nxref == 0)
THEN
1171 IF(sh3tree(3,n) >= 0)
THEN
1183 partsav(1,ip)=partsav(1,ip) + ems
1184 partsav(2,ip)=partsav(2,ip) + ems*
1185 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1186 partsav(3,ip)=partsav(3,ip) + ems*
1187 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1188 partsav(4,ip)=partsav(4,ip) + ems*
1189 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1190 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1191 . +x(1,i3)*x(1,i3)*p3)
1192 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1193 . +x(1,i3)*x(2,i3)*p3)
1194 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1195 . +x(2,i3)*x(2,i3)*p3)
1196 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1197 . +x(2,i3)*x(3,i3)*p3)
1198 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1199 . +x(3,i3)*x(3,i3)*p3)
1200 zx = (x(3,i1)*x(1,i1)*p1+x(3,i2)*x(1,i2)*p2
1201 . +x(3,i3)*x(1,i3)*p3)
1203 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1204 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1205 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1206 partsav(8,ip) =partsav(8,ip) - ems * xy
1207 partsav(9,ip) =partsav(9,ip) - ems * yz
1208 partsav(10,ip)=partsav(10,ip) - ems * zx
1210 partsav(11,ip)=partsav(11,ip) + ems*
1211 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1212 partsav(12,ip)=partsav(12,ip) + ems*
1213 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1214 partsav(13,ip)=partsav(13,ip) + ems*
1215 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1216 partsav(14,ip)=partsav(14,ip) + half * ems *
1217 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1218 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1219 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1225 IF (sh3tree(3,n) >= 0)
THEN
1237 partsav(1,ip)=partsav(1,ip) + ems
1238 partsav(2,ip)=partsav(2,ip) + ems*
1239 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)+p3*xreftg(3,1,i))
1240 partsav(3,ip)=partsav(3,ip) + ems*
1241 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1242 partsav(4,ip)=partsav(4,ip) + ems*
1243 . (p1*xreftg(1,3,i)+p2*xreftg(2,3,i)+p3*xreftg(3,3,i))
1244 xx = p1*xreftg(1,1,i)**2
1245 . + p2*xreftg(2,1,i)**2
1246 . + p3*xreftg(3,1,i)**2
1247 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1248 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1249 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1250 yy = p1*xreftg(1,2,i)**2
1251 . + p2*xreftg(2,2,i)**2
1252 . + p3*xreftg(3,2,i)**2
1253 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1254 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1255 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1256 zz = p1*xreftg(1,3,i)**2
1257 . + p2*xreftg(2,3,i)**2
1258 . + p3*xreftg(3,3,i)**2
1259 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1260 . + p2*xreftg(2,3,i)*xreftg(2,1,i)
1261 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1263 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1264 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1265 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1266 partsav(8,ip) =partsav(8,ip) - ems * xy
1267 partsav(9,ip) =partsav(9,ip) - ems * yz
1268 partsav(10,ip)=partsav(10,ip) - ems * zx
1270 partsav(11,ip)=partsav(11,ip) + ems*
1271 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1272 partsav(12,ip)=partsav(12,ip) + ems*
1273 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1274 partsav(13,ip)=partsav(13,ip) + ems*
1275 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1276 partsav(14,ip)=partsav(14,ip) + half * ems *
1277 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1278 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1279 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1284 IF (nxref == 0)
THEN
1287 IF(sh3tree(3,n) == 0 .OR. sh3tree(3,n) == -1 )
THEN
1299 partsav(1,ip)=partsav(1,ip) + ems
1300 partsav(2,ip)=partsav(2,ip) + ems*
1301 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1302 partsav(3,ip)=partsav(3,ip) + ems*
1303 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1304 partsav(4,ip)=partsav(4,ip) + ems*
1305 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1306 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1307 . +x(1,i3)*x(1,i3)*p3)
1308 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1309 . +x(1,i3)*x(2,i3)*p3)
1310 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1311 . +x(2,i3)*x(2,i3)*p3)
1312 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1313 . +x(2,i3)*x(3,i3)*p3)
1314 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1315 . +x(3,i3)*x(3,i3)*p3)
1316 zx = (x(3,i1)*x(1,i1)*p1+x(3,i2)*x(1,i2)*p2
1317 . +x(3,i3)*x(1,i3)*p3)
1319 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1320 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1321 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1322 partsav(8,ip) =partsav(8,ip) - ems * xy
1323 partsav(9,ip) =partsav(9,ip) - ems * yz
1324 partsav(10,ip)=partsav(10,ip) - ems * zx
1326 partsav(11,ip)=partsav(11,ip) + ems*
1327 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1328 partsav(12,ip)=partsav(12,ip) + ems*
1329 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1330 partsav(13,ip)=partsav(13,ip) + ems*
1331 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1332 partsav(14,ip)=partsav(14,ip) + half * ems *
1333 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1334 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1335 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1341 IF(sh3tree(3,n) == 0 .OR. sh3tree(3,n) == -1 )
THEN
1353 partsav(1,ip)=partsav(1,ip) + ems
1354 partsav(2,ip)=partsav(2,ip) + ems*
1355 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)+p3*xreftg(3,1,i))
1356 partsav(3,ip)=partsav(3,ip) + ems*
1357 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1358 partsav(4,ip)=partsav(4,ip) + ems*
1359 . (p1*xreftg(1,3,i)+p2*xreftg(2,3,i)+p3*xreftg(3,3,i))
1360 xx = p1*xreftg(1,1,i)**2
1361 . + p2*xreftg(2,1,i)**2
1362 . + p3*xreftg(3,1,i)**2
1363 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1364 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1365 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1366 yy = p1*xreftg(1,2,i)**2
1367 . + p2*xreftg(2,2,i)**2
1368 . + p3*xreftg(3,2,i)**2
1369 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1370 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1371 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1372 zz = p1*xreftg(1,3,i)**2
1373 . + p2*xreftg(2,3,i)**2
1374 . + p3*xreftg(3,3,i)**2
1375 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1376 . + p2*xreftg(2,3,i)*xreftg
1377 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1379 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1380 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1381 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1382 partsav(8,ip) =partsav(8,ip) - ems * xy
1383 partsav(9,ip) =partsav(9,ip) - ems * yz
1384 partsav(10,ip)=partsav(10,ip) - ems * zx
1386 partsav(11,ip)=partsav(11,ip) + ems*
1387 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1388 partsav(12,ip)=partsav(12,ip) + ems*
1389 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1390 partsav(13,ip)=partsav(13,ip) + ems*
1391 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1392 partsav(14,ip)=partsav(14,ip) + half * ems *
1393 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1394 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1395 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1402 IF(nintemp > 0 )
THEN
1404 IF(temp(ix1(i))== zero) temp(ix1(i)) = pm(79,imat)
1405 IF(temp(ix2(i))== zero) temp(ix2(i)) = pm(79,imat)
1406 IF(temp(ix3(i))== zero) temp(ix3(i)) = pm(79,imat)
1410 temp(ix1(i))= pm(79,imat)
1411 temp(ix2(i))= pm(79,imat)
1412 temp(ix3(i))= pm(79,imat)
1419 IF (igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))
THEN
1422 et=stack%PM(2,isubstack)*thk(i)
1424 nshnod(ix1(i))=nshnod(ix1(i))+1
1425 nshnod(ix2(i))=nshnod(ix2(i))+1
1426 nshnod(ix3(i))=nshnod(ix3(i))+1
1429 ELSEIF(igtyp == 11 .AND. igmat > 0)
THEN
1433 et=geo(ipgmat +2 ,iprop)*thk(i)
1435 nshnod(ix1(i))=nshnod(ix1(i))+1
1436 nshnod(ix2(i))=nshnod(ix2(i))+1
1437 nshnod(ix3(i))=nshnod(ix3(i))+1
1443 et=pm(20,imat)*thk(i)
1445 nshnod(ix1(i))=nshnod(ix1(i))+1
1446 nshnod(ix2(i))=nshnod(ix2(i))+1
1447 nshnod(ix3(i))=nshnod(ix3(i))+1
1451 DEALLOCATE(thk_it,pos_it, laypos,ratio_thkly )