46
47
48
49 USE elbufdef_mod
54 use element_mod , only : nixtg
55
56
57
58
59
60#include "implicit_f.inc"
61
62
63
64#include "mvsiz_p.inc"
65
66
67
68#include "com04_c.inc"
69#include "param_c.inc"
70#include "remesh_c.inc"
71#include "scr03_c.inc"
72#include "scr12_c.inc"
73#include "scr17_c.inc"
74#include "vect01_c.inc"
75#include "drape_c.inc"
76
77
78
79 INTEGER (NIXTG,*),IPART(*),IMAT,IPROP,
80 . IGEO(NPROPGI,*),IGTYP, NSHNOD(*), SH3TREE(KSH3TREE,*),
81 . SH3TRIM(*),ISUBSTACK,NLAY,IDRAPE,
82 . PERTURB(NPERTURB)
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
94
95
96
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,
100 . NLAY_MAX
101 INTEGER ID
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
111
112 TYPE (DRAPE_PLY_) , POINTER :: DRAPE_PLY
113
115 . a_gauss(9,9),w_gauss(9,9)
116
117 DATA a_gauss /
118 1 0. ,0. ,0. ,
119 1 0. ,0. ,0. ,
120 1 0. ,0. ,0. ,
121 2 -.577350269189626,0.577350269189626,0. ,
122 2 0. ,0. ,0. ,
123 2 0. ,0. ,0. ,
124 3 -.774596669241483,0. ,0.774596669241483,
125 3 0. ,0. ,0. ,
126 3 0. ,0. ,0. ,
127 4 -.861136311594053,-.339981043584856,0.339981043584856,
128 4 0.861136311594053,0. ,0. ,
129 4 0. ,0. ,0. ,
130 5 -.906179845938664,-.538469310105683,0. ,
131 5 0.538469310105683,0.906179845938664,0. ,
132 5 0. ,0. ,0. ,
133 6 -.932469514203152,-.661209386466265,-.238619186083197,
134 6 0.238619186083197,0.661209386466265,0.932469514203152,
135 6 0. ,0. ,0. ,
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/
145 DATA w_gauss /
146 1 2. ,0. ,0. ,
147 1 0. ,0. ,0. ,
148 1 0. ,0. ,0. ,
149 2 1. ,1. ,0. ,
150 2 0. ,0. ,0. ,
151 2 0. ,0. ,0. ,
152 3 0.555555555555556,0.888888888888889,0.555555555555556,
153 3 0. ,0. ,0. ,
154 3 0. ,0. ,0. ,
155 4 0.347854845137454,0.652145154862546,0.652145154862546,
156 4 0.347854845137454,0. ,0. ,
157 4 0. ,0. ,0. ,
158 5 0.236926885056189,0.478628670499366,0.568888888888889,
159 5 0.478628670499366,0.236926885056189,0. ,
160 5 0. ,0. ,0. ,
161 6 0.171324492379170,0.360761573048139,0.467913934572691,
162 6 0.467913934572691,0.360761573048139,0.171324492379170,
163 6 0. ,0. ,0. ,
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/
173
174 igtyp = nint(geo(12,iprop))
175 igmat = igeo(98,iprop)
176 max_nptt = 1
177 IF(igtyp == 51 .OR. igtyp == 52) THEN
178 DO ilay=1,nlay
179 max_nptt =
max(max_nptt ,elbuf_str%BUFLY(ilay)%NPTT)
180 ENDDO
181 ALLOCATE(thk_it(max_nptt),pos_it(max_nptt))
182 ELSE
183 ALLOCATE(thk_it(0),pos_it(0))
184 ENDIF
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
188 DO i=lft,llt
189 IF (ndrape == 0) THEN
190 IF (thke(i)== zero) THEN
191 thk(i) = stack%GEO(1,isubstack)
192 thke(i)= thk(i)
193 ELSE
194 thk(i)=thke(i)
195 ENDIF
196 thki(i) = stack%GEO(1,isubstack)
197 ELSEIF (ndrape > 0) THEN
198 thk(i) = thke(i)
199 thki(i)= thk(i)
200 ENDIF
201
202 IF (nperturb /= 0) THEN
203 DO j=1,nperturb
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)
209 ENDIF
210 ENDDO
211 ENDIF
212
213 rho(i)= pm(89,imat)
214 rhocp(i) = pm(69,imat)
215 IF(thk(i)==zero.AND.igtyp/=0)THEN
217 . msgtype=msgerror,
218 . anmode=aninfo_blind_1,
219 . i1=ixtg(6,i))
220 ENDIF
221 ENDDO
222 ELSEIF (igtyp == 52 .OR.
223 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
224 DO i=lft,llt
225 IF (ndrape == 0) THEN
226 IF (thke(i)== zero) THEN
227 thk(i) = stack%GEO(1,isubstack)
228 thke(i)= thk(i)
229 ELSE
230 thk(i)=thke(i)
231 ENDIF
232 thki(i) = stack%GEO(1,isubstack)
233 ELSEIF (ndrape > 0) THEN
234 thk(i) = thke(i)
235 thki(i)= thk(i)
236 ENDIF
237
238 IF (nperturb /= 0) THEN
239 DO j=1,nperturb
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)
245 ENDIF
246 ENDDO
247 ENDIF
248
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
253 . msgtype=msgerror,
254 . anmode=aninfo_blind_1,
255 . i1=ixtg(6,i))
256 ENDIF
257 ENDDO
258 ELSE
259 DO i=lft,llt
260 IF(thke(i)==zero)THEN
261 thk(i) =geo(1,iprop)
262 thke(i)=thk(i)
263 ELSE
264 thk(i)=thke(i)
265 ENDIF
266 thki(i) = thk(i)
267
268 IF (nperturb /= 0) THEN
269 DO j=1,nperturb
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)
275 ENDIF
276 ENDDO
277 ENDIF
278
279 rho(i)= pm(89,imat)
280 rhocp(i) = pm(69,imat)
281 IF(thk(i)==zero.AND.igtyp/=0)THEN
283 . msgtype=msgerror,
284 . anmode=aninfo_blind_1,
285 . i1=ixtg(6,i))
286 ENDIF
287 ENDDO
288 ENDIF
289
290
291
292 IF (igtyp==11 .OR. igtyp==16) THEN
293 ipthk = 300
294 ippos = 400
295 ipmat = 100
296 DO i=lft,llt
297 em(i)=zero
298 DO n=1,npt
299 i1=ipthk+n
300 i2=ipmat+n
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)
305 ENDDO
306
307 mst = rho(i)*thk(i)*
area(i)
308 dms = abs(em(i)-mst)/mst
309
310
311
312 IF(igmat <= 0) THEN
313 IF (dms > em02) THEN
314 idprop = nint(geo(40,iprop))
316 . igeo(npropgi-ltitr+1,iprop),ltitr)
318 . msgtype=msgwarning,
319 . anmode=aninfo_blind_2,
320 . i1=igeo(1,iprop),
321 . c1=titr,i2=ixtg(nixtg,nft+i),
322 . r1=em(i),r2=mst)
323 ENDIF
324 ENDIF
325 ENDDO
326 ELSEIF (igtyp == 17) THEN
327 ippid = 2
328 ipmat = ippid + npt
329 ipang = 1
330 ipthk = ipang + npt
331 ippos = ipthk + npt
332 nthk = ippos + npt
333 IF(idrape == 0) THEN
334 DO i=lft,llt
335 em(i)=zero
336 DO n=1,npt
337 i1=ipthk + n
338 i2=ippid + n
339 i3=ipmat + n
340 thickt = stack%GEO(1,isubstack)
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)
344 ENDDO
345
346 mst = rho(i)*thk(i)*
area(i)
347 dms = abs(em(i)-mst)/mst
348 IF(igmat < 0) THEN
349 IF (dms > em02) THEN
350 idprop = nint(geo(40,iprop))
352 . igeo(npropgi-ltitr+1,iprop),ltitr)
354 . msgtype=msgwarning,
355 . anmode=aninfo_blind_2,
356 . i1=igeo(1,iprop),
357 . c1=titr,i2=ixtg(nixtg,nft+i),
358 . r1=em(i),r2=mst)
359 ENDIF
360 ENDIF
361 ENDDO
362 ELSE
363 DO i=lft,llt
364 em(i)=zero
365 ie = indx(nft + i)
366 IF(ie == 0) THEN
367 DO n=1,npt
368 i1=ipthk + n
369 i2=ippid + n
370 i3=ipmat + n
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)
375 ENDDO
376 ELSE
377 thick_drape = drape(ie)%THICK
378 scale = thk(i)/thick_drape
379 DO n=1,npt
380 i1=ipthk + n
381 i2=ippid + n
382 i3=ipmat + n
383 ip = drape(ie)%INDX_PLY(n)
384 IF(ip > 0) THEN
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,isubstack)
391 thkly = thkly*thinning*scale
392 matly = stack%IGEO(i3,isubstack)
393 em(i) = em(i) + pm(89,matly)*thkly*
area(i)
394 ELSE
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)
400 ENDIF
401 ENDDO
402 ENDIF
403
404 mst = rho(i)*thk(i)*
area(i)
405 dms = abs(em(i)-mst)/mst
406 IF(igmat < 0) THEN
407 IF (dms > em02) THEN
408 idprop = nint(geo(40,iprop))
410 . igeo(npropgi-ltitr+1,iprop),ltitr)
412 . msgtype=msgwarning,
413 . anmode=aninfo_blind_2,
414 . i1=igeo(1,iprop),
415 . c1=titr,i2=ixtg(nixtg,nft+i),
416 . r1=em(i),r2=mst)
417 ENDIF
418 ENDIF
419 ENDDO
420 ENDIF
421 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
422 ipang = 1
423 ippid = 2
424 ipmat = ippid + nlay
425 ipthk = ipang + nlay
426 ippos = ipthk + nlay
427 nthk = ippos + nlay
428
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)
433 IF(idrape == 0) THEN
434 DO i=lft,llt
435 scale = thk(i)/thickt
436 em(i)=zero
437 ipt_all = 0
438 DO ilay=1,nlay
439 nptt = elbuf_str%BUFLY(ilay)%NPTT
440 i1=ipthk + ilay
441 i2=ippid + ilay
442 i3=ipmat + ilay
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)
448 iint = igeo(47,iprop)
449 IF (iint == 1) THEN
450 DO it=1,nptt
451 thk_it(it) = thkly/nptt
452 ENDDO
453 ELSEIF (iint == 2) THEN
454 DO it=1,nptt
455 thk_it(it) = half*thkly*w_gauss(it,nptt)
456 ENDDO
457 ENDIF
458
459 DO it=1,nptt
460 thk_it(it) = thk_it(it)*scale
461 ems_nptt = pm(89,matly)*thk_it(it)*
area(i)
462 em(i) = em(i) + ems_nptt
463 ENDDO
464
465 ENDDO
466
467 mst = rho(i)*thk(i)*
area(i)
468 dms = abs(em(i)-mst)/mst
469 IF(igtyp == 51 .AND.igmat < 0) THEN
470 IF (dms > em02) THEN
471 idprop = nint(geo(40,iprop))
473 . igeo(npropgi-ltitr+1,iprop),ltitr)
475 . msgtype=msgwarning,
476 . anmode=aninfo_blind_2,
477 . i1=igeo(1,iprop),
478 . c1=titr,i2=ixtg(nixtg,nft+i),
479 . r1=em(i),r2=mst)
480 ENDIF
481 ENDIF
482 ENDDO
483 ELSE
484 thickt = stack%GEO(1,isubstack)
485 DO i=lft,llt
486 em(i)=zero
487 ie = indx(nft + i)
488 ipt_all = 0
489 scale = thk(i)/thickt
490 IF(ie == 0) THEN
491 DO ilay=1,nlay
492 nptt = elbuf_str%BUFLY(ilay)%NPTT
493 i1=ipthk + ilay
494 i2=ippid + ilay
495 i3=ipmat + ilay
496 ipid = stack%IGEO(i2,isubstack)
497 iint = igeo(47,iprop)
498 thickt = stack%GEO(1,isubstack)
499 IF(iint == 1) THEN
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)
504 DO it=1,nptt
505 thk_it(it) = thkly/nptt
506 ENDDO
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)
512 DO it=1,nptt
513 thk_it(it) = half*thkly*w_gauss(it,nptt)
514 ENDDO
515 ENDIF
516
517 DO it=1,nptt
518 thk_it(it) = scale*thk_it(it)
519 ems_nptt = pm(89,matly)*thk_it(it)*
area(i)
520 em(i) = em(i) + ems_nptt
521 ENDDO
522 ENDDO
523 ELSE
524 ipt_all = 0
525 thick_drape = drape(ie)%THICK
526 scale = thk(i)/thick_drape
527 DO ilay=1,nlay
528 nptt = elbuf_str%BUFLY(ilay)%NPTT
529 i1=ipthk + ilay
530 i2=ippid + ilay
531 i3=ipmat + ilay
532 ipid = stack%IGEO(i2,isubstack) ! or ipid = ix(nixc-1,nft+i)
533 iint = igeo(47,iprop)
534 thickt = stack%GEO(1,isubstack)
535 ip = drape(ie)%INDX_PLY(ilay)
536 IF(ip == 0) THEN
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)
541 IF(iint == 1) THEN
542 DO it=1,nptt
543 thk_it(it) = thkly/nptt
544 ENDDO
545 ELSEIF(iint == 2)THEN
546 DO it=1,nptt
547 thk_it(it) = half*thkly*w_gauss(it,nptt)
548 ENDDO
549 ENDIF
550 ELSE
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)
557 IF(iint == 1) THEN
558 DO it=1,nptt
559 thinning = drape_ply%RDRAPE(it,1)
560 thk_it(it) = thkly*thinning/nptt
561 ENDDO
562 ELSEIF(iint == 2)THEN
563 DO it=1,nptt
564 thinning = drape_ply%RDRAPE(it,1)
565 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
566 ENDDO
567 ENDIF
568 ENDIF
569 DO it=1,nptt
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
573 ENDDO
574 ENDDO
575 ENDIF
576
577 mst = rho(i)*thk(i)*
area(i)
578 dms = abs(em(i)-mst)/mst
579 IF(igtyp == 51 .AND.igmat < 0) THEN
580 IF (dms > em02) THEN
581 idprop = nint(geo(40,iprop))
583 . igeo(npropgi-ltitr+1,iprop),ltitr)
585 . msgtype=msgwarning,
586 . anmode=aninfo_blind_2,
587 . i1=igeo(1,iprop),
588 . c1=titr,i2=ixtg(nixtg,nft+i),
589 . r1=em(i),r2=mst)
590 ENDIF
591 ENDIF
592 ENDDO
593 ENDIF
594 ELSE
595 DO i=lft,llt
596 em(i) = rho(i)*thk(i)*
area(i)
597 ENDDO
598 DO i=lft,llt
599 emscp(i)=rhocp(i)*thk(i)*
area(i)
600 ENDDO
601 ENDIF
602
603
604
605 IF (igtyp==11 .OR. igtyp==16) THEN
606 DO i=lft,llt
607 xi(i)=zero
608 DO n=1,npt
609 i1=ipthk+n
610 i2=ipmat+n
611 i3=ippos+n
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)
618 xi(i) = xi(i) + xil
619 ENDDO
620
621 IF(igmat <= 0) THEN
622 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
623 dms = abs(xi(i)-ins)/ins
624 IF(dms > em02) THEN
625 idprop = nint(geo(40,iprop))
627 . igeo(npropgi-ltitr+1,iprop),ltitr)
629 . msgtype=msgwarning,
630 . anmode=aninfo_blind_2,
631 . i1=igeo(1,iprop),
632 . c1=titr,i2=ixtg(nixtg,nft+i))
633 ENDIF
634 ENDIF
635 ENDDO
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
642 DO i=lft,llt
643 xi(i)=zero
644 DO n=1,npt
645 i1 =ippid + n
646 i2 =ipmat + n
647 i3 =ipthk + n
648 i4 =ippos + n
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)
655 xi(i) = xi(i) + xil
656 ENDDO
657
658 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
659 dms = abs(xi(i)-ins)/ins
660
661 IF(igmat <= 0) THEN
662 IF(dms > em02) THEN
663 idprop = nint(geo(40,iprop))
665 . igeo(npropgi-ltitr+1,iprop),ltitr)
667 . msgtype=msgwarning,
668 . anmode=aninfo_blind_2,
669 . i1=igeo(1,iprop),
670 . c1=titr,i2=ixtg(nixtg,nft+i))
671 ENDIF
672 ENDIF
673 ENDDO
674 ELSE
675 DO i=lft,llt
676 xi(i)=zero
677 ie = indx(nft + i)
678 IF(ie == 0) THEN
679 DO n=1,npt
680 i1 =ippid + n
681 i2 =ipmat + n
682 i3 =ipthk + n
683 i4 =ippos + n
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)
690 xi(i) = xi(i) + xil
691 ENDDO
692 ELSE
693 thick_drape = drape(ie)%THICK
694 scale = thk(i)/thick_drape
695 DO n=1,npt
696 ip = drape(ie)%INDX_PLY(n)
697 i1 =ippid + n
698 i2 =ipmat + n
699 i3 =ipthk + n
700 i4 =ippos + n
701 IF(ip == 0 ) THEN
702 thkly = stack%GEO(i3,isubstack)*thickt
703 matly = stack%IGEO(i2,isubstack)
704 ratio_thkly(n,i) = thkly/ thick_drape
705 ELSE
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
712 ENDIF
713 IF (n == 1) THEN
714 laypos(n,i) = zshift + half*ratio_thkly(n,i)
715 ELSE
716 laypos(n,i) = laypos(n-1,i) + half*(ratio_thkly(n,i)+ratio_thkly(n-1,i))
717 ENDIF
718 posly = laypos(n,i)*thk(i)
719 thkly = thkly*scale
720 msl = pm(89,matly)*thkly*
area(i)
721 xil = msl*(
area(i)/nine_over_2+thkly*thkly
722 xi(i) = xi(i) + xil
723 ENDDO
724 ENDIF
725
726 ins = em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
727 dms = abs(xi(i)-ins)/ins
728
729 IF(igmat <= 0) THEN
730 IF(dms > em02) THEN
731 idprop = nint(geo(40,iprop))
733 . igeo(npropgi-ltitr+1,iprop),ltitr)
735 . msgtype=msgwarning,
736 . anmode=aninfo_blind_2,
737 . i1=igeo(1,iprop),
738 . c1=titr,i2=ixtg(nixtg,nft+i))
739 ENDIF
740 ENDIF
741 ENDDO
742 ENDIF
743 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
744 ipang = 1
745 ippid = 2
746 ipmat = ippid + nlay
747 ipthk = ipang + nlay
748 ippos = ipthk + nlay
749 nthk = ippos + nlay
750
751
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)
756 IF(idrape == 0) THEN
757 thickt = stack%GEO(1,isubstack)
758 DO i=lft,llt
759 xi(i)=zero
760 ipt_all = 0
761 scale = thk(i)/thickt
762 DO ilay=1,nlay
763 nptt = elbuf_str%BUFLY(ilay)%NPTT
764 i1 =ippid + ilay
765 i2 =ipmat + ilay
766 i3 =ipthk + ilay
767 i4 =ippos + ilay
768 ipid_ly = stack%IGEO(i1,isubstack)
769 matly = stack%IGEO(i2,isubstack)
770 iint = igeo(47,iprop)
771 IF(iint == 1) THEN
772 thickt = stack%GEO(1,isubstack)
773 thkly = stack%GEO(i3,isubstack)*thickt
774 DO it=1,nptt
775 ipt = ipt_all + it
776 thk_it(it) = thkly/nptt
777 ratio_thkly(ipt,i) = thk_it(it)/thickt
778 IF (ipt == 1) THEN
779 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
780 ELSE
781 laypos(ipt,i) = laypos(ipt-1,i) +
782 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
783 ENDIF
784 pos_it(it) = laypos(ipt,i
785 ENDDO
786 ELSEIF(iint == 2) THEN
787 thickt = stack%GEO(1,isubstack)
788 thkly = stack%GEO(i3,isubstack)*thickt
789 DO it=1,nptt
790 ipt = ipt_all + it
791 thk_it(it) = half*thkly*w_gauss(it,nptt)
792 ratio_thkly(ipt,i) = thk_it(it)/thickt
793 IF (ipt == 1) THEN
794 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
795 ELSE
796 laypos(ipt,i) = laypos(ipt-1,i) +
797 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt
798 ENDIF
799 pos_it(it) = laypos(ipt,i)*thk(i)
800 ENDDO
801 ENDIF
802
803 DO it=1,nptt
804 ms_nptt = pm(89,matly)*thk_it(it)*
area(i)
805 thk_it(it) = thk_it(it)*scale
806 xil = ms_nptt*(
area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
807 . pos_it(it)*pos_it(it))
808 xi(i) = xi(i) + xil
809 ENDDO
810 ipt_all = ipt_all + nptt
811 ENDDO
812
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
816 IF (dms > em02) THEN
817 idprop = nint(geo(40,iprop))
819 . igeo(npropgi-ltitr+1,iprop),ltitr)
821 . msgtype=msgwarning,
822 . anmode=aninfo_blind_2,
823 . i1=igeo(1,iprop),
824 . c1=titr,i2=ixtg(nixtg,nft+i))
825 ENDIF
826 ENDIF
827 ENDDO
828 ELSE
829 thickt = stack%GEO(1,isubstack)
830 DO i=lft,llt
831 xi(i)=zero
832 ipt_all = 0
833 ie = indx(nft + i)
834 scale = thk(i)/thickt
835 IF(ie == 0) THEN
836 DO ilay=1,nlay
837 nptt = elbuf_str%BUFLY(ilay)%NPTT
838 i1 =ippid + ilay
839 i2 =ipmat + ilay
840 i3 =ipthk + ilay
841 i4 =ippos + ilay
842 ipid_ly = stack%IGEO(i1,isubstack)
843 matly = stack%IGEO(i2,isubstack)
844 iint = igeo(47,iprop)
845 IF(iint == 1) THEN
846 thickt = stack%GEO(1,isubstack)
847 thkly = stack%GEO(i3,isubstack)*thickt
848 DO it=1,nptt
849 ipt = ipt_all + it
850 thk_it(it) = thkly/nptt
851 ratio_thkly(ipt,i) = thk_it(it)/thickt
852 IF (ipt == 1) THEN
853 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
854 ELSE
855 laypos(ipt,i) = laypos(ipt-1,i) +
856 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
857 ENDIF
858 pos_it(it) = laypos(ipt,i)*thk(i)
859 ENDDO
860 ELSEIF(iint == 2) THEN
861 thickt = stack%GEO(1,isubstack)
862 thkly = stack%GEO(i3,isubstack)*thickt
863 DO it=1,nptt
864 ipt = ipt_all + it
865 thk_it(it) = half*thkly*w_gauss(it,nptt)
866 ratio_thkly(ipt,i) = thk_it(it)/thickt
867 IF (ipt == 1) THEN
868 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
869 ELSE
870 laypos(ipt,i) = laypos(ipt-1,i) +
871 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
872 ENDIF
873 pos_it(it) = laypos(ipt,i)*thk(i)
874 ENDDO
875 ENDIF
876 DO it=1,nptt
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))
881 xi(i) = xi(i) + xil
882 ENDDO
883 ipt_all = ipt_all + nptt
884 ENDDO
885 ELSE
886 ipt_all = 0
887 thick_drape = drape(ie)%THICK
888 scale = thk(i)/thick_drape
889 DO ilay=1,nlay
890 ip = drape(ie)%INDX_PLY(ilay)
891 nptt = elbuf_str%BUFLY(ilay)%NPTT
892 i1 =ippid + ilay
893 i2 =ipmat + ilay
894 i3 =ipthk + ilay
895 i4 =ippos + ilay
896 ipid_ly = stack%IGEO(i1,isubstack)
897 matly = stack%IGEO(i2,isubstack)
898 iint = igeo(47,iprop)
899 IF(ip > 0) THEN
900 drape_ply => drape(ie)%DRAPE_PLY(ip)
901 nslice = drape_ply%NSLICE
902 IF(iint == 1) THEN
903 thickt = stack%GEO(1,isubstack)
904 thkly = stack%GEO(i3,isubstack)*thickt
905 DO it=1,nptt
906 ipt = ipt_all + it
907 thinning = drape_ply%RDRAPE(it,1)
908 thk_it(it) = thkly*thinning/nptt
909
910 IF (ipt == 1) THEN
911 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
912 ELSE
913 laypos(ipt,i) = laypos(ipt-1,i) +
914 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
915 ENDIF
916 pos_it(it) = laypos(ipt,i)*thk(i)
917 ENDDO
918 ELSEIF(iint == 2) THEN
919 thickt = stack%GEO(1,isubstack)
920 thkly = stack%GEO(i3,isubstack)*thickt
921 DO it=1,nptt
922 ipt = ipt_all + it
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
926 IF (ipt == 1) THEN
927 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
928 ELSE
929 laypos(ipt,i) = laypos(ipt-1,i) +
930 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
931 ENDIF
932 pos_it(it) = laypos(ipt,i)*thk(i)
933 ENDDO
934 ENDIF
935 ELSE
936 IF(iint == 1) THEN
937 thickt = stack%GEO(1,isubstack)
938 thkly = stack%GEO(i3,isubstack)*thickt
939 DO it=1,nptt
940 ipt = ipt_all + it
941 thk_it(it) = thkly/nptt
942 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
943 IF (ipt == 1) THEN
944 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
945 ELSE
946 laypos(ipt,i) = laypos(ipt-1,i) +
947 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
948 ENDIF
949 pos_it(it) = laypos(ipt,i)*thk(i)
950 ENDDO
951 ELSEIF(iint == 2) THEN
952 thickt = stack%GEO(1,isubstack)
953 thkly = stack%GEO(i3,isubstack)*thickt
954 DO it=1,nptt
955 ipt = ipt_all + it
956 thk_it(it) = half*thkly*w_gauss(it,nptt)
957 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
958 IF (ipt == 1) THEN
959 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
960 ELSE
961 laypos(ipt,i) = laypos(ipt-1,i) +
962 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
963 ENDIF ! IF (n == 1)
964 pos_it(it) = laypos(ipt,i)*thk(i)
965 ENDDO
966 ENDIF
967 ENDIF
968 DO it=1,nptt
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))
973 xi(i) = xi(i) + xil
974 ENDDO
975 ipt_all = ipt_all + nptt
976 ENDDO
977 ENDIF
978
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
982 IF (dms > em02) THEN
983 idprop = nint(geo(40,iprop))
985 . igeo(npropgi-ltitr+1,iprop),ltitr)
987 . msgtype=msgwarning,
988 . anmode=aninfo_blind_2,
989 . i1=igeo(1,iprop),
990 . c1=titr,i2=ixtg(nixtg,nft+i))
991 ENDIF
992 ENDIF
993 ENDDO
994 ENDIF
995 ELSE
996 DO i=lft,llt
997
998
999 xi(i)=em(i)*(
area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
1000 ENDDO
1001 ENDIF
1002
1003
1004
1005
1006 IF(jthe == 0 ) THEN
1007 DO i=lft,llt
1008 mstg(i)=em(i)
1009 intg(i)=xi(i)
1010 END DO
1011 ELSE
1012 DO i=lft,llt
1013 mstg(i)=em(i)
1014 intg(i)=xi(i)
1015 mcptg(i) = emscp(i)
1016 END DO
1017 ENDIF
1018
1019 DO i=lft,llt
1020 aa = x2(i)
1021 a2 = aa**2
1022 b2 = (x2(i)-x3(i))**2 + y3(i)**2
1023 bb = sqrt(b2)
1024 c2 = x3(i)**2 + y3(i)**2
1025 cc = sqrt(c2)
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
1029 ptg(1,i)=p1
1030 ptg(2,i)=p2
1031 ptg(3,i)=p3
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,
1040 . msgtype=msgerror,
1041 . anmode=aninfo,
1042 . i1=ixtg(6,nft+i))
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
1051 . msgtype=msgerror,
1052 . anmode=aninfo,
1053 . i1=ixtg(6,nft+i))
1054 ENDIF
1055 END DO
1056
1057
1058 IF(nadmesh==0)THEN
1059 IF (nxref == 0) THEN
1060 DO i=lft,llt
1061 p1 = ptg(1,i)
1062 p2 = ptg(2,i)
1063 p3 = ptg(3,i)
1064 ems = em(i)
1065 xin = xi(i)
1066
1067 i1 = ix1(i)
1068 i2 = ix2(i)
1069 i3 = ix3(i)
1070
1071 ip=ipart(i)
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+x(3,i2)*x(1,i2)*p2
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
1097
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) + half * ems *
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)
1108 ENDDO
1109 ELSE
1110 DO i=lft,llt
1111 p1 = ptg(1,i)
1112 p2 = ptg(2,i)
1113 p3 = ptg(3,i)
1114 ems = em(i)
1115 xin = xi(i)
1116
1117 i1 = ix1(i)
1118 i2 = ix2(i)
1119 i3 = ix3(i)
1120 ip=ipart(i)
1121
1122 partsav(1,ip)=partsav(1,ip) + ems
1123 partsav(2,ip)=partsav(2,ip) + ems*
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)
1147
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
1154
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(3,i3)*v(3,i3))*p3)
1165 ENDDO
1166 ENDIF
1167 ELSEIF (istatcnd==0)THEN
1168 IF (nxref == 0) THEN
1169 DO i=lft,llt
1170 n=nft+i
1171 IF(sh3tree(3,n) >= 0)THEN
1172 p1 = ptg(1,i)
1173 p2 = ptg(2,i)
1174 p3 = ptg(3,i)
1175 ems = em(i)
1176 xin = xi(i)
1177
1178 i1 = ix1(i)
1179 i2 = ix2(i)
1180 i3 = ix3(i)
1181 ip = ipart(i)
1182
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)
1202
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
1209
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)
1220 END IF
1221 END DO
1222 ELSE
1223 DO i=lft,llt
1224 n=nft+i
1225 IF (sh3tree(3,n) >= 0)THEN
1226 p1 = ptg(1,i)
1227 p2 = ptg(2,i)
1228 p3 = ptg(3,i)
1229 ems = em(i)
1230 xin = xi(i)
1231
1232 i1 = ix1(i)
1233 i2 = ix2(i)
1234 i3 = ix3(i)
1235 ip = ipart(i)
1236
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)
1262
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
1269
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)
1280 END IF
1281 END DO
1282 ENDIF
1283 ELSE
1284 IF (nxref == 0) THEN
1285 DO i=lft,llt
1286 n=nft+i
1287 IF(sh3tree(3,n) == 0 .OR. sh3tree(3,n) == -1 )THEN
1288 p1 = ptg(1,i)
1289 p2 = ptg(2,i)
1290 p3 = ptg(3,i)
1291 ems = em(i)
1292 xin = xi(i)
1293
1294 i1 = ix1(i)
1295 i2 = ix2(i)
1296 i3 = ix3(i)
1297 ip = ipart(i)
1298
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)
1318
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
1325
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)
1336 END IF
1337 END DO
1338 ELSE
1339 DO i=lft,llt
1340 n=nft+i
1341 IF(sh3tree(3,n) == 0 .OR. sh3tree(3,n) == -1 )THEN
1342 p1
1343 p2 = ptg(2,i)
1344 p3 = ptg(3,i)
1345 ems = em(i)
1346 xin = xi(i)
1347
1348 i1 = ix1(i)
1349 i2 = ix2(i)
1350 i3 = ix3(i)
1351 ip = ipart(i)
1352
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(2,1,i)
1377 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1378
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
1385
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)
1396 END IF
1397 END DO
1398 ENDIF
1399 END IF
1400
1401 IF(jthe > 0 ) THEN
1402 IF(nintemp > 0 ) THEN
1403 DO i= lft,llt
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)
1407 ENDDO
1408 ELSE
1409 DO i=lft,llt
1410 temp(ix1(i))= pm(79,imat)
1411 temp(ix2(i))= pm(79,imat)
1412 temp(ix3(i))= pm(79,imat)
1413 ENDDO
1414 ENDIF
1415 ENDIF
1416
1417
1418
1419 IF (igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))THEN
1420 IF(i7stifs/=0)THEN
1421 DO i=lft,llt
1422 et=stack%PM(2,isubstack)*thk(i)
1423 sttg(i)=et
1424 nshnod(ix1(i))=nshnod(ix1(i))+1
1425 nshnod(ix2(i))=nshnod(ix2(i))+1
1426 nshnod(ix3(i))=nshnod(ix3(i))+1
1427 ENDDO
1428 ENDIF
1429 ELSEIF(igtyp == 11 .AND. igmat > 0) THEN
1430 ipgmat = 700
1431 IF(i7stifs/=0)THEN
1432 DO i=lft,llt
1433 et=geo(ipgmat +2 ,iprop)*thk(i)
1434 sttg(i)=et
1435 nshnod(ix1(i))=nshnod(ix1(i))+1
1436 nshnod(ix2(i))=nshnod(ix2(i))+1
1437 nshnod(ix3(i))=nshnod(ix3(i))+1
1438 ENDDO
1439 ENDIF
1440 ElSE
1441 IF(i7stifs/=0)THEN
1442 DO i=lft,llt
1443 et=pm(20,imat)*thk(i)
1444 sttg(i)=et
1445 nshnod(ix1(i))=nshnod(ix1(i))+1
1446 nshnod(ix2(i))=nshnod(ix2(i))+1
1447 nshnod(ix3(i))=nshnod(ix3(i))+1
1448 ENDDO
1449 ENDIF
1450 ENDIF
1451 DEALLOCATE(thk_it,pos_it, laypos,ratio_thkly )
1452
1453 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
integer, parameter nchartitle
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)