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)
59#include "implicit_f.inc"
69#include "remesh_c.inc"
73#include "vect01_c.inc"
78 INTEGER IXTG(NIXTG,*),IPART(*),IMAT,IPROP,
79 . IGEO(NPROPGI,*),IGTYP, NSHNOD(*), SH3TREE(KSH3TREE
82INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ)
83 INTEGER,
DIMENSION (STDRAPE) :: INDX
84 INTEGER ,
INTENT(IN) :: NINTEMP
86 . x(3,*), geo(npropg,*), pm(npropm,*), ms(*), tiner(*),thke(*),
87 . v(3,*),partsav(20,*),mstg(*),intg(*),ptg(3,*),
88 . etnod(*), sttg(*),mcp(*),mcptg(*),temp(*),xreftg(3,3,*),thki(*),
89 . rnoise(nperturb,*),
area(mvsiz)
90 TYPE(elbuf_struct_),
TARGET :: ELBUF_STR
91 TYPE (STACK_PLY) :: STACK
92 TYPE (DRAPE_),
DIMENSION(NUMELC_DRAPE + NUMELTG_DRAPE),
TARGET :: DRAPE
96 INTEGER I,J,K,N,IP,I1,I2,I3,I4,IPTHK,IPMAT,IPPOS,MATLY,IDPROP,
97 . ipid,ippid,iigeo,iadi,iadr,ipang,nthk,igmat,ilay,nptt,
98 . it,max_nptt,iint,ipid_ly,ipgmat,ipos,ipt,ipt_all,nslice,ie,
101 CHARACTER(LEN=NCHARTITLE) :: TITR
102 my_real XX,YY,ZZ,XY,YZ,ZX,ZSHIFT
103 my_real x2(mvsiz), x3(mvsiz), y3(mvsiz),
104 . rho(mvsiz), thk(mvsiz), em(mvsiz), xi(mvsiz),
105 . p1, p2, p3, aa, bb, cc, a2, b2, c2, xin, ems, msl, xil,
106 . mst,dms,ins,posly,thkly,thickt, et,emscp(mvsiz),rhocp(mvsiz),
107 . ems_nptt,ms_nptt,pos_0,thinning,scale,thick_drape
108 my_real,
DIMENSION(:) ,
ALLOCATABLE :: thk_it, pos_it
109 my_real,
DIMENSION(:,:),
ALLOCATABLE :: laypos, ratio_thkly
111 TYPE (DRAPE_PLY_) ,
POINTER :: DRAPE_PLY
114 . A_GAUSS(9,9),W_GAUSS(9,9)
120 2 -.577350269189626,0.577350269189626,0.
123 3 -.774596669241483,0. ,0.774596669241483,
126 4 -.861136311594053,-.339981043584856,0.339981043584856,
127 4 0.861136311594053,0. ,0. ,
129 5 -.906179845938664,-.538469310105683,0. ,
130 5 0.538469310105683,0.906179845938664,0. ,
132 6 -.932469514203152,-.661209386466265,-.238619186083197,
133 6 0.238619186083197,0.661209386466265,0.932469514203152,
135 7 -.949107912342759,-.741531185599394,-.405845151377397,
136 7 0. ,0.405845151377397,0.741531185599394,
137 7 0.949107912342759,0. ,0. ,
138 8 -.960289856497536,-.796666477413627,-.525532409916329,
139 8 -.183434642495650,0.183434642495650,0.525532409916329,
140 8 0.796666477413627,0.960289856497536,0. ,
141 9 -.968160239507626,-.836031107326636,-.613371432700590,
142 9 -.324253423403809,0. ,0.324253423403809,
143 9 0.613371432700590,0.836031107326636,0.968160239507626/
151 3 0.555555555555556,0.888888888888889,0.555555555555556,
154 4 0.347854845137454,0.652145154862546,0.652145154862546,
155 4 0.347854845137454,0. ,0. ,
157 5 0.236926885056189,0.478628670499366,0.568888888888889,
158 5 0.478628670499366,0.236926885056189,0. ,
160 6 0.171324492379170,0.360761573048139,0.467913934572691,
161 6 0.467913934572691,0.360761573048139,0.171324492379170,
163 7 0.129484966168870,0.279705391489277,0.381830050505119,
164 7 0.417959183673469,0.381830050505119,0.279705391489277,
165 7 0.129484966168870,0. ,0. ,
166 8 0.101228536290376,0.222381034453374,0.313706645877887,
167 8 0.362683783378362,0.362683783
168 8 0.222381034453374,0.101228536290376,0. ,
169 9 0.081274388361574,0.180648160694857,0.260610696402935,
170 9 0.312347077040003,0.330239355001260,0.312347077040003,
171 9 0.260610696402935,0.180648160694857,0.081274388361574/
173 igtyp = nint(geo(12,iprop))
174 igmat = igeo(98,iprop)
176 IF(igtyp == 51 .OR. igtyp == 52)
THEN
178 max_nptt =
max(max_nptt ,elbuf_str%BUFLY(ilay)%NPTT)
180 ALLOCATE(thk_it(max_nptt),pos_it(max_nptt))
182 ALLOCATE(thk_it(0),pos_it(0))
184 nlay_max =
max(nlay, npt)
185 ALLOCATE( laypos(nlay_max*max_nptt,mvsiz),ratio_thkly(nlay_max*max_nptt,mvsiz) )
186 IF ((igtyp == 17 .OR. igtyp == 51) .AND. igmat < 0)
THEN
188 IF (ndrape == 0)
THEN
189 IF (thke(i)== zero)
THEN
190 thk(i) = stack%GEO(1,isubstack)
195 thki(i) = stack%GEO(1,isubstack)
196 ELSEIF (ndrape > 0)
THEN
201 IF (nperturb /= 0)
THEN
203 IF (perturb(j) == 1 .AND.
204 . rnoise(j,numelc+i+nft) /= zero)
THEN
205 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
206 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
207 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
213 rhocp(i) = pm(69,imat)
214 IF(thk(i)==zero.AND.igtyp/=0)
THEN
217 . anmode=aninfo_blind_1,
221 ELSEIF (igtyp == 52 .OR.
222 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))
THEN
224 IF (ndrape == 0)
THEN
225 IF (thke(i)== zero)
THEN
226 thk(i) = stack%GEO(1,isubstack)
231 thki(i) = stack%GEO(1,isubstack)
232 ELSEIF (ndrape > 0)
THEN
237 IF (nperturb /= 0)
THEN
239 IF (perturb(j) == 1 .AND.
240 . rnoise(j,numelc+i+nft) /= zero)
THEN
241 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
242 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
243 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
248 rho(i) = stack%PM(1,isubstack)
249 rhocp(i) = stack%PM(15,isubstack)* rho(i)/stack%PM(14,isubstack)
250 IF (thk(i) == zero .AND. igtyp /= 0)
THEN
253 . anmode=aninfo_blind_1,
259 IF(thke(i)==zero)
THEN
265 thki(i) =geo(1,iprop)
267 IF (nperturb /= 0)
THEN
269 IF (perturb(j) == 1 .AND.
270 . rnoise(j,numelc+i+nft) /= zero)
THEN
271 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
272 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
273 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
279 rhocp(i) = pm(69,imat)
280 IF(thk(i)==zero.AND.igtyp/=0)
THEN
283 . anmode=aninfo_blind_1,
291 IF (igtyp==11 .OR. igtyp==16)
THEN
300 thickt = geo(200,iprop)
301 thkly = geo(i1,iprop)*thk(i)
302 matly = igeo(i2,iprop)
303 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
306 mst = rho(i)*thk(i)*
area(i)
307 dms = abs(em(i)-mst)/mst
315 . igeo(npropgi-ltitr+1,iprop),ltitr)
317 . msgtype=msgwarning,
318 . anmode=aninfo_blind_2,
320 . c1=titr,i2=ixtg(nixtg,nft+i),
325 ELSEIF (igtyp == 17)
THEN
339 thickt = stack%GEO(1,isubstack)
340 thkly = stack%GEO(i1,isubstack)*thk(i)
341 matly = stack%IGEO(i3,isubstack)
342 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
345 mst = rho(i)*thk(i)*
area(i)
346 dms = abs(em(i)-mst)/mst
349 idprop = nint(geo(40,iprop))
351 . igeo(npropgi-ltitr+1,iprop),ltitr)
353 . msgtype=msgwarning,
354 . anmode=aninfo_blind_2,
356 . c1=titr,i2=ixtg(nixtg,nft+i),
370 thickt = stack%GEO(1,isubstack)
371 thkly = stack%GEO(i1,isubstack)*thk(i)
372 matly = stack%IGEO(i3,isubstack)
373 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
376 thick_drape = drape(ie)%THICK
377 scale = thk(i)/thick_drape
382 ip = drape(ie)%INDX_PLY(n)
384 drape_ply => drape(ie)%DRAPE_PLY
385 nslice = drape_ply%NSLICE
386 thinning = drape_ply%RDRAPE(1,1)
387 thickt = stack%GEO(1,isubstack)
388 thkly = stack%GEO(i1,isubstack)*thickt
389 matly = stack%IGEO(i3,isubstack)
390 thkly = thkly*thinning*scale
391 matly = stack%IGEO(i3,isubstack)
392 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
394 thickt = stack%GEO(1,isubstack)
395 thkly = stack%GEO(i1,isubstack)*thickt*scale
396 matly = stack%IGEO(i3,isubstack)
397 matly = stack%IGEO(i3,isubstack)
398 em(i) = em(i) + pm(89,matly)*thkly*
area(i
403 mst = rho(i)*thk(i)*
area(i)
404 dms = abs(em(i)-mst)/mst
407 idprop = nint(geo(40,iprop))
409 . igeo(npropgi-ltitr+1,iprop),ltitr)
411 . msgtype=msgwarning,
412 . anmode=aninfo_blind_2,
414 . c1=titr,i2=ixtg(nixtg,nft+i),
420 ELSEIF (igtyp == 51 .OR. igtyp == 52)
THEN
428 ipos = igeo(99,iprop)
429 thickt = stack%GEO(1,isubstack)
430 zshift = geo(199,iprop)
431 IF(ipos == 2 ) zshift = - zshift /
max(thickt,em20)
434 scale = thk(i)/thickt
438 nptt = elbuf_str%BUFLY(ilay)%NPTT
442 thickt = stack%GEO(1,isubstack)
443 thkly = stack%GEO(i1,isubstack)*thickt
444 ipid_ly= stack%IGEO(i2,isubstack)
445 matly = stack%IGEO(i3,isubstack)
446 ipid = stack%IGEO(ippid,isubstack) ! or ipid = ix(nixc-1,nft+i)
447 iint = igeo(47,iprop)
450 thk_it(it) = thkly/nptt
452 ELSEIF (iint == 2)
THEN
454 thk_it(it) = half*thkly*w_gauss(it,nptt)
459 thk_it(it) = thk_it(it)*scale
460 ems_nptt = pm(89,matly)*thk_it(it)*
area(i)
461 em(i) = em(i) + ems_nptt
466 mst = rho(i)*thk(i)*
area(i)
467 dms = abs(em(i)-mst)/mst
468 IF(igtyp == 51 .AND.igmat < 0)
THEN
470 idprop = nint(geo(40,iprop))
472 . igeo(npropgi-ltitr+1,iprop
474 . msgtype=msgwarning,
475 . anmode=aninfo_blind_2,
477 . c1=titr,i2=ixtg(nixtg,nft+i),
483 thickt = stack%GEO(1,isubstack)
488 scale = thk(i)/thickt
491 nptt = elbuf_str%BUFLY(ilay)%NPTT
495 ipid = stack%IGEO(i2,isubstack)
496 iint = igeo(47,iprop)
497 thickt = stack%GEO(1,isubstack)
499 thickt = stack%GEO(1,isubstack)
500 thkly = stack%GEO(i1,isubstack)*thickt
501 ipid_ly= stack%IGEO(i2,isubstack)
502 matly = stack%IGEO(i3,isubstack)
504 thk_it(it) = thkly/nptt
506 ELSEIF(iint == 2)
THEN
507 thickt = stack%GEO(1,isubstack)
508 thkly = stack%GEO(i1,isubstack)*thickt
509 ipid_ly= stack%IGEO(i2,isubstack)
510 matly = stack%IGEO(i3,isubstack)
512 thk_it(it) = half*thkly*w_gauss(it,nptt)
517 thk_it(it) = scale*thk_it(it)
518 ems_nptt = pm(89,matly)*thk_it(it)*
area(i)
519 em(i) = em(i) + ems_nptt
524 thick_drape = drape(ie)%THICK
525 scale = thk(i)/thick_drape
527 nptt = elbuf_str%BUFLY(ilay)%NPTT
532 iint = igeo(47,iprop)
533 thickt = stack%GEO(1,isubstack)
534 ip = drape(ie)%INDX_PLY(ilay)
536 thickt = stack%GEO(1,isubstack)
537 thkly = stack%GEO(i1,isubstack)*thickt
538 ipid_ly= stack%IGEO(i2,isubstack)
539 matly = stack%IGEO(i3,isubstack)
542 thk_it(it) = thkly/nptt
544 ELSEIF(iint == 2)
THEN
546 thk_it(it) = half*thkly*w_gauss(it,nptt)
550 drape_ply => drape(ie)%DRAPE_PLY(ip)
551 nslice = drape_ply%NSLICE
552 thickt = stack%GEO(1,isubstack)
553 thkly = stack%GEO(i1,isubstack)*thickt
554 ipid_ly= stack%IGEO(i2
555 matly = stack%IGEO(i3,isubstack)
558 thinning = drape_ply%RDRAPE(it,1)
559 thk_it(it) = thkly*thinning/nptt
561 ELSEIF(iint == 2)
THEN
563 thinning = drape_ply%RDRAPE(it,1)
564 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
570 ems_nptt = pm(89,matly)*thk_it(it)*
area(i)
571 em(i) = em(i) + ems_nptt
576 mst = rho(i)*thk(i)*
area(i)
577 dms = abs(em(i)-mst)/mst
578 IF(igtyp == 51 .AND.igmat < 0)
THEN
580 idprop = nint(geo(40,iprop))
582 . igeo(npropgi-ltitr+1,iprop),ltitr)
584 . msgtype=msgwarning,
585 . anmode=aninfo_blind_2,
587 . c1=titr,i2=ixtg(nixtg,nft+i),
595 em(i) = rho(i)*thk(i)*
area(i)
598 emscp(i)=rhocp(i)*thk(i)*
area(i)
604 IF (igtyp==11 .OR. igtyp==16)
THEN
611 thickt= geo(200,iprop)
612 thkly = geo(i1,iprop)*thk(i)
613 posly = geo(i3,iprop)*thk(i)
614 matly = igeo(i2,iprop)
615 msl = pm(89,matly)*thkly*
area(i)
616 xil = msl*(
area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
621 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
622 dms = abs(xi(i)-ins)/ins
624 idprop = nint(geo(40,iprop))
626 . igeo(npropgi-ltitr+1,iprop),ltitr)
628 . msgtype=msgwarning,
629 . anmode=aninfo_blind_2,
631 . c1=titr,i2=ixtg(nixtg,nft+i))
635 ELSEIF (igtyp == 17)
THEN
636 ipos = igeo(99,iprop)
637 thickt = stack%GEO(1,isubstack)
638 zshift = geo(199,iprop)
639 IF(ipos == 2 ) zshift = - zshift /
max(thickt,em20)
640 IF(idrape == 0 )
THEN
648 thickt= stack%GEO(1,isubstack)
649 thkly = stack%GEO(i3,isubstack)*thk(i)
650 posly = stack%GEO(i4,isubstack)*thk(i)
651 matly = stack%IGEO(i2,isubstack)
652 msl = pm(89,matly)*thkly*
area(i)
653 xil = msl*(
area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
657 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
658 dms = abs(xi(i)-ins)/ins
662 idprop = nint(geo(40,iprop))
664 . igeo(npropgi-ltitr+1,iprop),ltitr)
666 . msgtype=msgwarning,
667 . anmode=aninfo_blind_2,
669 . c1=titr,i2=ixtg(nixtg,nft+i))
683 thickt= stack%GEO(1,isubstack
684 thkly = stack%GEO(i3,isubstack)*thk(i)
685 posly = stack%GEO(i4,isubstack)*thk(i)
686 matly = stack%IGEO(i2,isubstack)
687 msl = pm(89,matly)*thkly*
area(i)
688 xil = msl*(
area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
692 thick_drape = drape(ie)%THICK
693 scale = thk(i)/thick_drape
695 ip = drape(ie)%INDX_PLY(n)
701 thkly = stack%GEO(i3,isubstack)*thickt
702 matly = stack%IGEO(i2,isubstack)
703 ratio_thkly(n,i) = thkly/ thick_drape
705 drape_ply => drape(ie)%DRAPE_PLY(ip)
706 thinning = drape_ply%RDRAPE(1,1)
707 thkly = stack%GEO(i3,isubstack)*thickt
708 thkly = thkly*thinning
709 matly = stack%IGEO(i2,isubstack)
710 ratio_thkly(n,i) = thkly/thick_drape
713 laypos(n,i) = zshift + half*ratio_thkly(n,i)
715 laypos(n,i) = laypos(n-1,i) + half*(ratio_thkly(n,i)+ratio_thkly(n-1,i))
717 posly = laypos(n,i)*thk(i)
719 msl = pm(89,matly)*thkly*
area(i)
720 xil = msl*(
area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
725 ins = em(i)*(
area(i)/nine_over_2
726 dms = abs(xi(i)-ins)/ins
730 idprop = nint(geo(40,iprop))
732 . igeo(npropgi-ltitr+1,iprop),ltitr)
734 . msgtype=msgwarning,
735 . anmode=aninfo_blind_2,
737 . c1=titr,i2=ixtg(nixtg,nft+i))
742 ELSEIF (igtyp == 51 .OR. igtyp == 52)
THEN
751 ipos = igeo(99,iprop)
752 thickt = stack%GEO(1,isubstack)
753 zshift = geo(199,iprop)
754 IF(ipos == 2 ) zshift = - zshift /
max(thickt,em20)
756 thickt = stack%GEO(1,isubstack)
760 scale = thk(i)/thickt
762 nptt = elbuf_str%BUFLY(ilay)%NPTT
767 ipid_ly = stack%IGEO(i1,isubstack)
768 matly = stack%IGEO(i2,isubstack)
769 iint = igeo(47,iprop)
771 thickt = stack%GEO(1,isubstack)
772 thkly = stack%GEO(i3,isubstack)*thickt
775 thk_it(it) = thkly/nptt
776 ratio_thkly(ipt,i) = thk_it(it)/thickt
778 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
780 laypos(ipt,i) = laypos(ipt-1,i) +
781 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
783 pos_it(it) = laypos(ipt,i)*thk(i)
785 ELSEIF(iint == 2)
THEN
786 thickt = stack%GEO(1,isubstack)
787 thkly = stack%GEO(i3,isubstack)*thickt
790 thk_it(it) = half*thkly*w_gauss(it,nptt)
791 ratio_thkly(ipt,i) = thk_it(it)/thickt
793 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
795 laypos(ipt,i) = laypos(ipt-1,i) +
796 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
798 pos_it(it) = laypos(ipt,i)*thk(i)
804 thk_it(it) = thk_it(it)*scale
805 xil = ms_nptt*(
area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
806 . pos_it(it)*pos_it(it))
809 ipt_all = ipt_all + nptt
812 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
813 dms = abs(xi(i)-ins)/ins
814 IF(igtyp == 51 .AND. igmat < 0)
THEN
816 idprop = nint(geo(40,iprop))
818 . igeo(npropgi-ltitr+1,iprop),ltitr)
820 . msgtype=msgwarning,
821 . anmode=aninfo_blind_2,
823 . c1=titr,i2=ixtg(nixtg,nft+i))
828 thickt = stack%GEO(1,isubstack)
833 scale = thk(i)/thickt
836 nptt = elbuf_str%BUFLY(ilay)%NPTT
841 ipid_ly = stack%IGEO(i1,isubstack)
842 matly = stack%IGEO(i2,isubstack)
843 iint = igeo(47,iprop)
845 thickt = stack%GEO(1,isubstack)
846 thkly = stack%GEO(i3,isubstack)*thickt
849 thk_it(it) = thkly/nptt
850 ratio_thkly(ipt,i) = thk_it(it)/thickt
852 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
854 laypos(ipt,i) = laypos(ipt-1,i) +
855 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
857 pos_it(it) = laypos(ipt,i)*thk(i)
859 ELSEIF(iint == 2)
THEN
860 thickt = stack%GEO(1,isubstack)
861 thkly = stack%GEO(i3,isubstack)*thickt
864 thk_it(it) = half*thkly*w_gauss(it,nptt)
865 ratio_thkly(ipt,i) = thk_it(it)/thickt
869 laypos(ipt,i) = laypos(ipt-1,i) +
870 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
872 pos_it(it) = laypos(ipt,i)*thk(i)
876 thk_it(it) = scale*thk_it(it)
877 ms_nptt = pm(89,matly)*thk_it(it)*
area(i)
878 xil = ms_nptt*(
area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
882 ipt_all = ipt_all + nptt
886 thick_drape = drape(ie)%THICK
887 scale = thk(i)/thick_drape
889 ip = drape(ie)%INDX_PLY(ilay)
890 nptt = elbuf_str%BUFLY(ilay)%NPTT
895 ipid_ly = stack%IGEO(i1,isubstack)
896 matly = stack%IGEO(i2,isubstack)
897 iint = igeo(47,iprop)
899 drape_ply => drape(ie)%DRAPE_PLY(ip)
900 nslice = drape_ply%NSLICE
902 thickt = stack%GEO(1,isubstack)
903 thkly = stack%GEO(i3,isubstack)*thickt
906 thinning = drape_ply%RDRAPE(it,1)
907 thk_it(it) = thkly*thinning/nptt
908 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
912 laypos(ipt,i) = laypos(ipt-1,i) +
913 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
915 pos_it(it) = laypos(ipt,i)*thk(i)
917 ELSEIF(iint == 2)
THEN
918 thickt = stack%GEO(1,isubstack)
919 thkly = stack%GEO(i3,isubstack)*thickt
922 thinning = drape_ply%RDRAPE(it,1)
923 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
924 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
926 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
928 laypos(ipt,i) = laypos(ipt-1,i) +
929 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
931 pos_it(it) = laypos(ipt,i)*thk(i)
936 thickt = stack%GEO(1,isubstack)
937 thkly = stack%GEO(i3,isubstack)*thickt
940 thk_it(it) = thkly/nptt
941 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
943 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
945 laypos(ipt,i) = laypos(ipt-1,i) +
946 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
948 pos_it(it) = laypos(ipt,i)*thk(i)
950 ELSEIF(iint == 2)
THEN
951 thickt = stack%GEO(1,isubstack)
952 thkly = stack%GEO(i3,isubstack)*thickt
955 thk_it(it) = half*thkly*w_gauss(it,nptt)
956 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
958 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
960 laypos(ipt,i) = laypos(ipt-1,i
961 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
963 pos_it(it) = laypos(ipt,i)*thk(i)
968 thk_it(it) = scale*thk_it(it)
969 ms_nptt = pm(89,matly)*thk_it(it)*
area(i)
970 xil = ms_nptt*(
area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
971 . pos_it(it)*pos_it(it))
974 ipt_all = ipt_all + nptt
978 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
979 dms = abs(xi(i)-ins)/ins
980 IF(igtyp == 51 .AND. igmat < 0)
THEN
982 idprop = nint(geo(40,iprop))
984 . igeo(npropgi-ltitr+1,iprop),ltitr)
986 . msgtype=msgwarning,
987 . anmode=aninfo_blind_2,
989 . c1=titr,i2=ixtg(nixtg,nft+i))
998 xi(i)=em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
1021 b2 = (x2(i)-x3(i))**2 + y3(i)**2
1023 c2 = x3(i)**2 + y3(i)**2
1025 p1 = acos((a2 + c2 - b2)/(two * aa * cc)) / pi
1026 p2 = acos((a2 + b2 - c2)/(two * aa * bb)) / pi
1027 p3 = acos((b2 + c2 - a2)/(two * bb * cc)) / pi
1031 IF ( ( ((a2 + c2 - b2)/(two * aa * cc) <= -one) .OR.
1032 . ( (a2 + c2 - b2)/(two * aa * cc) >= one) .OR.
1033 . ( (a2 + b2 - c2)/(two * aa * bb) <= -one) .OR.
1034 . ( (a2 + b2 - c2)/(two * aa * bb) >= one) .OR.
1035 . ( (b2 + c2 - a2)/(two * bb * cc) <= -one) .OR.
1036 . ( (b2 + c2 - a2)/(two * bb * cc) >= one)) .AND.
1037 . lsh3trim /= 0 )
THEN
1038 IF (sh3trim(nft+i) /= -1)
CALL ancmsg(msgid=750,
1042 ELSEIF ( ( ((a2 + c2 - b2)/(two * aa * cc) <= -one) .OR.
1043 . ( (a2 + c2 - b2)/(two * aa * cc) >= one) .OR.
1044 . ( (a2 + b2 - c2)/(two * aa * bb) <= -one) .OR.
1045 . ( (a2 + b2 - c2)/(two * aa * bb) >= one) .OR.
1046 . ( (b2 + c2 - a2)/(two * bb * cc) <= -one) .OR.
1047 . ( (b2 + c2 - a2)/(two * bb * cc) >= one)) .AND.
1048 . lsh3trim == 0 )
THEN
1058 IF (nxref == 0)
THEN
1071 partsav(1,ip)=partsav(1,ip) + ems
1072 partsav(2,ip)=partsav(2,ip) + ems*
1073 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1074 partsav(3,ip)=partsav(3,ip) + ems*
1075 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1076 partsav(4,ip)=partsav(4,ip)
1077 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1078 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1079 . +x(1,i3)*x(1,i3)*p3)
1080 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1081 . +x(1,i3)*x(2,i3)*p3)
1082 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1083 . +x(2,i3)*x(2,i3)*p3)
1084 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1085 . +x(2,i3)*x(3,i3)*p3)
1086 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1087 . +x(3,i3)*x(3,i3)*p3)
1088 zx = (x(3,i1)*x(1,i1)*p1+x(3,i2)*x(1,i2)*p2
1089 . +x(3,i3)*x(1,i3)*p3)
1090 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1091 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1092 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1093 partsav(8,ip) =partsav(8,ip) - ems * xy
1094 partsav(9,ip) =partsav(9,ip) - ems * yz
1095 partsav(10,ip)=partsav(10,ip) - ems * zx
1097 partsav(11,ip)=partsav(11,ip) + ems*
1098 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1099 partsav(12,ip)=partsav(12,ip) + ems*
1100 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1101 partsav(13,ip)=partsav(13,ip) + ems*
1102 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1103 partsav(14,ip)=partsav(14,ip) + half * ems *
1104 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1105 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1106 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1121 partsav(1,ip)=partsav(1,ip) + ems
1122 partsav(2,ip)=partsav(2,ip) + ems*
1123 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)
1125 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1126 partsav(4,ip)=partsav(4,ip) + ems*
1128 xx = p1*xreftg(1,1,i)**2
1129 . + p2*xreftg(2,1,i)**2
1130 . + p3*xreftg(3,1,i)**2
1131 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1132 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1133 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1134 yy = p1*xreftg(1,2,i)**2
1135 . + p2*xreftg(2,2,i)**2
1136 . + p3*xreftg(3,2,i)**2
1137 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1138 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1139 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1140 zz = p1*xreftg(1,3,i)**2
1141 . + p2*xreftg(2,3,i)**2
1142 . + p3*xreftg(3,3,i)**2
1143 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1144 . + p2*xreftg(2,3,i)*xreftg(2,1,i)
1145 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1147 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1148 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1149 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1150 partsav(8,ip) =partsav(8,ip) - ems * xy
1151 partsav(9,ip) =partsav(9,ip) - ems * yz
1152 partsav(10,ip)=partsav(10,ip) - ems * zx
1154 partsav(11,ip)=partsav(11,ip) + ems*
1155 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1156 partsav(12,ip)=partsav(12,ip) + ems*
1157 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1158 partsav(13,ip)=partsav(13,ip) + ems*
1159 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1160 partsav(14,ip)=partsav(14,ip) + half * ems *
1161 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1162 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1163 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1166 ELSEIF (istatcnd==0)
THEN
1167 IF (nxref == 0)
THEN
1170 IF(sh3tree(3,n) >= 0)
THEN
1182 partsav(1,ip)=partsav(1,ip) + ems
1183 partsav(2,ip)=partsav(2,ip) + ems*
1184 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1185 partsav(3,ip)=partsav(3,ip) + ems*
1186 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1187 partsav(4,ip)=partsav(4,ip) + ems*
1188 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1189 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1190 . +x(1,i3)*x(1,i3)*p3)
1191 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1192 . +x(1,i3)*x(2,i3)*p3)
1193 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1194 . +x(2,i3)*x(2,i3)*p3)
1195 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1196 . +x(2,i3)*x(3,i3)*p3)
1197 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1198 . +x(3,i3)*x(3,i3)*p3)
1199 zx = (x(3,i1)*x(1,i1)*p1+x(3,i2)*x(1,i2)*p2
1200 . +x(3,i3)*x(1,i3)*p3)
1202 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1203 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1204 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1205 partsav(8,ip) =partsav(8,ip) - ems * xy
1206 partsav(9,ip) =partsav(9,ip) - ems * yz
1207 partsav(10,ip)=partsav(10,ip) - ems * zx
1209 partsav(11,ip)=partsav(11,ip) + ems*
1210 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1211 partsav(12,ip)=partsav(12,ip) + ems*
1212 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1213 partsav(13,ip)=partsav(13,ip) + ems*
1214 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1215 partsav(14,ip)=partsav(14,ip) + half * ems *
1216 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1217 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1218 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1224 IF (sh3tree(3,n) >= 0)
THEN
1236 partsav(1,ip)=partsav(1,ip) + ems
1237 partsav(2,ip)=partsav(2,ip) + ems*
1238 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)+p3*xreftg(3,1,i))
1239 partsav(3,ip)=partsav(3,ip) + ems*
1240 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1241 partsav(4,ip)=partsav(4,ip) + ems*
1242 . (p1*xreftg(1,3,i)+p2*xreftg(2,3,i)+p3*xreftg(3,3,i))
1243 xx = p1*xreftg(1,1,i)**2
1244 . + p2*xreftg(2,1,i)**2
1245 . + p3*xreftg(3,1,i)**2
1246 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1247 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1248 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1249 yy = p1*xreftg(1,2,i)**2
1250 . + p2*xreftg(2,2,i)**2
1251 . + p3*xreftg(3,2,i)**2
1252 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1253 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1254 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1255 zz = p1*xreftg(1,3,i)**2
1256 . + p2*xreftg(2,3,i)**2
1257 . + p3*xreftg(3,3,i)**2
1258 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1259 . + p2*xreftg(2,3,i)*xreftg(2,1,i)
1260 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1262 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1263 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1264 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1265 partsav(8,ip) =partsav(8,ip) - ems * xy
1266 partsav(9,ip) =partsav(9,ip) - ems * yz
1267 partsav(10,ip)=partsav(10,ip) - ems * zx
1269 partsav(11,ip)=partsav(11,ip) + ems*
1270 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1271 partsav(12,ip)=partsav(12,ip) + ems*
1272 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1273 partsav(13,ip)=partsav(13,ip) + ems
1274 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1275 partsav(14,ip)=partsav(14,ip) + half * ems *
1276 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1277 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1278 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1283 IF (nxref == 0)
THEN
1286 IF(sh3tree(3,n) == 0 .OR. sh3tree(3,n) == -1 )
THEN
1298 partsav(1,ip)=partsav(1,ip) + ems
1299 partsav(2,ip)=partsav(2,ip) + ems*
1300 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1301 partsav(3,ip)=partsav(3,ip) + ems*
1302 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1303 partsav(4,ip)=partsav(4,ip) + ems*
1304 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1305 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1306 . +x(1,i3)*x(1,i3)*p3)
1307 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1308 . +x(1,i3)*x(2,i3)*p3)
1309 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1310 . +x(2,i3)*x(2,i3)*p3)
1311 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1312 . +x(2,i3)*x(3,i3)*p3)
1313 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1314 . +x(3,i3)*x(3,i3)*p3)
1315 zx = (x(3,i1)*x(1,i1)*p1+x(3,i2)*x(1,i2)*p2
1316 . +x(3,i3)*x(1,i3)*p3)
1318 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1320 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1321 partsav(8,ip) =partsav(8,ip) - ems * xy
1322 partsav(9,ip) =partsav(9,ip) - ems * yz
1323 partsav(10,ip)=partsav(10,ip) - ems * zx
1325 partsav(11,ip)=partsav(11,ip) + ems*
1326 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1327 partsav(12,ip)=partsav(12,ip) + ems*
1328 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1329 partsav(13,ip)=partsav(13,ip) + ems*
1330 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1331 partsav(14,ip)=partsav(14,ip) + half * ems *
1332 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1333 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1334 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1340 IF(sh3tree(3,n) == 0 .OR. sh3tree(3,n) == -1 )
THEN
1352 partsav(1,ip)=partsav(1,ip) + ems
1353 partsav(2,ip)=partsav(2,ip) + ems*
1354 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)+p3*xreftg(3,1,i))
1355 partsav(3,ip)=partsav(3,ip) + ems*
1356 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1357 partsav(4,ip)=partsav(4,ip) + ems*
1358 . (p1*xreftg(1,3,i)+p2*xreftg(2,3,i)+p3*xreftg(3,3,i))
1359 xx = p1*xreftg(1,1,i)**2
1360 . + p2*xreftg(2,1,i)**2
1361 . + p3*xreftg(3,1,i)**2
1362 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1363 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1364 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1365 yy = p1*xreftg(1,2,i)**2
1366 . + p2*xreftg(2,2,i)**2
1367 . + p3*xreftg(3,2,i)**2
1368 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1369 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1370 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1371 zz = p1*xreftg(1,3,i)**2
1372 . + p2*xreftg(2,3,i)**2
1373 . + p3*xreftg(3,3,i)**2
1374 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1375 . + p2*xreftg(2,3,i)*xreftg(2,1,i)
1376 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1378 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1379 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1380 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1381 partsav(8,ip) =partsav(8,ip) - ems * xy
1382 partsav(9,ip) =partsav(9,ip) - ems * yz
1383 partsav(10,ip)=partsav(10,ip) - ems * zx
1385 partsav(11,ip)=partsav(11,ip) + ems*
1386 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1387 partsav(12,ip)=partsav(12,ip) + ems*
1388 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1389 partsav(13,ip)=partsav(13,ip) + ems*
1390 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1391 partsav(14,ip)=partsav(14,ip) + half * ems *
1392 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1393 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1394 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1401 IF(nintemp > 0 )
THEN
1403 IF(temp(ix1(i))== zero) temp(ix1(i)) = pm(79,imat)
1404 IF(temp(ix2(i))== zero) temp(ix2(i)) = pm(79,imat)
1405 IF(temp(ix3(i))== zero) temp(ix3(i)) = pm(79,imat)
1409 temp(ix1(i))= pm(79,imat)
1410 temp(ix2(i))= pm(79,imat)
1411 temp(ix3(i))= pm(79,imat)
1418 IF (igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))
THEN
1421 et=stack%PM(2,isubstack)*thk(i)
1423 nshnod(ix1(i))=nshnod(ix1(i))+1
1424 nshnod(ix2(i))=nshnod(ix2(i))+1
1425 nshnod(ix3(i))=nshnod(ix3(i))+1
1428 ELSEIF(igtyp == 11 .AND. igmat > 0)
THEN
1432 et=geo(ipgmat +2 ,iprop)*thk(i)
1434 nshnod(ix1(i))=nshnod(ix1(i))+1
1435 nshnod(ix2(i))=nshnod(ix2(i))+1
1436 nshnod(ix3(i))=nshnod(ix3(i))+1
1442 et=pm(20,imat)*thk(i)
1444 nshnod(ix1(i))=nshnod(ix1(i))+1
1445 nshnod(ix2(i))=nshnod(ix2(i))+1
1446 nshnod(ix3(i))=nshnod(ix3(i))+1
1450 DEALLOCATE(thk_it,pos_it, laypos,ratio_thkly )