OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cinmas.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| double_flot_ieee ../starter/source/elements/shell/coque/cinmas.F
25!||====================================================================
26 SUBROUTINE double_flot_ieee(JFT ,JLT ,I8 ,R8, I8F)
27C-----------------------------------------------
28C I m p l i c i t T y p e s
29C-----------------------------------------------
30#include "implicit_f.inc"
31C-----------------------------------------------
32C G l o b a l P a r a m e t e r s
33C-----------------------------------------------
34#include "mvsiz_p.inc"
35C-----------------------------------------------
36C D u m m y A r g u m e n t s
37C-----------------------------------------------
38 INTEGER JFT, JLT
39 integer*8 I8(*),I8F(3,*)
40 my_real
41 . r8(mvsiz)
42C-----------------------------------------------
43C L o c a l V a r i a b l e s
44C-----------------------------------------------
45c___________________________________________________
46 double precision
47 . r8_local,r8_deuxp43,aa
48 INTEGER*8 I8_DEUXP43
49 DATA i8_deuxp43 /'80000000000'x/
50 DATA r8_deuxp43 /'42A0000000000000'x/
51 INTEGER I
52c___________________________________________________
53C-----------------------------------------------
54C
55 DO i=jft,jlt
56c___________________________________________________
57 i8f(1,i) = r8(i)
58 aa = i8f(1,i)
59 r8_local = (r8(i) - aa) * r8_deuxp43
60 i8f(2,i) = r8_local
61 aa = i8f(2,i)
62 r8_local = (r8_local - aa) * r8_deuxp43
63 i8f(3,i) = r8_local + half
64 ENDDO
65c___________________________________________________
66 RETURN
67 END
68C
69!||====================================================================
70!|| cinmas ../starter/source/elements/shell/coque/cinmas.F
71!||--- called by ------------------------------------------------------
72!|| cbainit3 ../starter/source/elements/shell/coqueba/cbainit3.F
73!|| cinit3 ../starter/source/elements/shell/coque/cinit3.F
74!|| inirig_mat ../starter/source/elements/initia/inirig_mat.F
75!|| inivoid ../starter/source/elements/initia/inivoid.F
76!||--- calls -----------------------------------------------------
77!|| ancmsg ../starter/source/output/message/message.F
78!|| fretitl2 ../starter/source/starter/freform.F
79!||--- uses -----------------------------------------------------
80!|| drape_mod ../starter/share/modules1/drape_mod.F
81!|| message_mod ../starter/share/message_module/message_mod.F
82!|| stack_mod ../starter/share/modules1/stack_mod.F
83!||====================================================================
84 SUBROUTINE cinmas(
85 . X ,XREFC ,IX ,GEO ,PM ,
86 . MS ,TINER ,THKE ,IHBE ,PARTSAV ,
87 . V ,IPART ,MSC ,INC ,AREA ,
88 . I8MI ,IGEO ,ETNOD ,IMID ,IPROP ,
89 . NSHNOD ,STC ,SH4TREE ,MCP ,MCPS ,
90 . TEMP ,MS_LAYER,ZI_LAYER,MS_LAYERC,ZI_LAYERC,
91 . MSZ2C ,ZPLY ,ISUBSTACK,NLAY ,ELBUF_STR,
92 . STACK ,THKI ,RNOISE ,DRAPE ,NINTEMP ,
93 . PERTURB ,IX1 ,IX2 ,IX3 ,IX4 ,
94 . IDRAPE ,INDX)
95C-----------------------------------------------
96C M o d u l e s
97C-----------------------------------------------
98 USE elbufdef_mod
99 USE message_mod
100 USE stack_mod
101 USE drape_mod
102 USE names_and_titles_mod , ONLY : nchartitle
103 use element_mod , only : nixc
104C----------------------------------------------
105C initialization of nodal masses and inertias
106C----------------------------------------------
107C I m p l i c i t T y p e s
108C-----------------------------------------------
109#include "implicit_f.inc"
110C-----------------------------------------------
111C G l o b a l P a r a m e t e r s
112C-----------------------------------------------
113#include "mvsiz_p.inc"
114C-----------------------------------------------
115C C o m m o n B l o c k s
116C-----------------------------------------------
117#include "com01_c.inc"
118#include "com04_c.inc"
119#include "remesh_c.inc"
120#include "param_c.inc"
121#include "scr03_c.inc"
122#include "scr12_c.inc"
123#include "scr17_c.inc"
124#include "scr22_c.inc"
125#include "vect01_c.inc"
126#include "drape_c.inc"
127C-----------------------------------------------
128C D u m m y A r g u m e n t s
129C-----------------------------------------------
130 INTEGER IX(NIXC,*),IPART(*),IHBE,IMID,IPROP,
131 . IGEO(NPROPGI,*),NSHNOD(*), SH4TREE(KSH4TREE,*),
132 . ISUBSTACK,NLAY,PERTURB(NPERTURB),IDRAPE
133 INTEGER, DIMENSION(MVSIZ), INTENT(IN) :: IX1,IX2,IX3,IX4
134 INTEGER ,INTENT(IN) :: NINTEMP
135 INTEGER*8 I8MI(6,*)
136 INTEGER, DIMENSION(SCDRAPE) :: INDX
137 my_real, DIMENSION(MVSIZ), INTENT(IN) :: AREA
138 my_real
139 . X(3,*), GEO(NPROPG,*), PM(NPROPM,*), MS(*), TINER(*),THKE(*),
140 . v(3,*),partsav(20,*),msc(*),inc(*),
141 . etnod(*), stc(*),mcp(*),mcps(*),temp(*),
142 . ms_layer(numnod,*), zi_layer(numnod,*),xrefc(4,3,*),
143 . ms_layerc(numelc,*),zi_layerc(numelc,*),msz2c(*),zply(*),thki(*),
144 . rnoise(nperturb,*)
145 my_real
146 . zshift
147 TYPE(elbuf_struct_), TARGET :: ELBUF_STR
148 TYPE (STACK_PLY) :: STACK
149 TYPE (DRAPE_) , DIMENSION(NUMELC_DRAPE+NUMELTG_DRAPE), TARGET :: DRAPE
150C-----------------------------------------------
151C L o c a l V a r i a b l e s
152C-----------------------------------------------
153 INTEGER I,II,J,N,IP,I1,I2,I3,I4,MATLY,IPTHK,IPMAT,IPPOS,
154 . ERRM,IDPROP,IPID,IPPID,ISTACK(MVSIZ,NPT),ISLV,IPID_LY,
155 . ipang,nthk,igtyp,igmat,ilay,nptt,it,iint,ipgmat,max_nptt,
156 . nslice,ipos,ipt,ipt_all,ie,nlay_max
157 CHARACTER(LEN=NCHARTITLE) :: TITR
158 my_real
159 . fac,mst,dms,thickt,thkly,laypos,ins,mzi,ems_nptt,et,zj,
160 . xx,yy,zz,xy,yz,zx,msl,xil,ms_nptt,thinning,f_offset,scale,
161 . thick_drape
162 my_real
163 . ems_layer(mvsiz,nlay)
164 my_real, DIMENSION(:), ALLOCATABLE :: posly,ratio_thkly, thk_it,pos_it
165 my_real, DIMENSION(MVSIZ) :: rho,rhocp,ems,emscp,xi,thk
166C
167
168 TYPE (DRAPE_PLY_), POINTER :: DRAPE_PLY
169 my_real A_GAUSS(9,9),W_GAUSS(9,9)
170C-----------------------------------------------
171 DATA A_GAUSS /
172 1 0. ,0. ,0. ,
173 1 0. ,0. ,0. ,
174 1 0. ,0. ,0. ,
175 2 -.577350269189626,0.577350269189626,0. ,
176 2 0. ,0. ,0. ,
177 2 0. ,0. ,0. ,
178 3 -.774596669241483,0. ,0.774596669241483,
179 3 0. ,0. ,0. ,
180 3 0. ,0. ,0. ,
181 4 -.861136311594053,-.339981043584856,0.339981043584856,
182 4 0.861136311594053,0. ,0. ,
183 4 0. ,0. ,0. ,
184 5 -.906179845938664,-.538469310105683,0. ,
185 5 0.538469310105683,0.906179845938664,0. ,
186 5 0. ,0. ,0. ,
187 6 -.932469514203152,-.661209386466265,-.238619186083197,
188 6 0.238619186083197,0.661209386466265,0.932469514203152,
189 6 0. ,0. ,0. ,
190 7 -.949107912342759,-.741531185599394,-.405845151377397,
191 7 0. ,0.405845151377397,0.741531185599394,
192 7 0.949107912342759,0. ,0. ,
193 8 -.960289856497536,-.796666477413627,-.525532409916329,
194 8 -.183434642495650,0.183434642495650,0.525532409916329,
195 8 0.796666477413627,0.960289856497536,0. ,
196 9 -.968160239507626,-.836031107326636,-.613371432700590,
197 9 -.324253423403809,0. ,0.324253423403809,
198 9 0.613371432700590,0.836031107326636,0.968160239507626/
199 DATA w_gauss /
200 1 2. ,0. ,0. ,
201 1 0. ,0. ,0. ,
202 1 0. ,0. ,0. ,
203 2 1. ,1. ,0. ,
204 2 0. ,0. ,0. ,
205 2 0. ,0. ,0. ,
206 3 0.555555555555556,0.888888888888889,0.555555555555556,
207 3 0. ,0. ,0. ,
208 3 0. ,0. ,0. ,
209 4 0.347854845137454,0.652145154862546,0.652145154862546,
210 4 0.347854845137454,0. ,0. ,
211 4 0. ,0. ,0. ,
212 5 0.236926885056189,0.478628670499366,0.568888888888889,
213 5 0.478628670499366,0.236926885056189,0. ,
214 5 0. ,0. ,0. ,
215 6 0.171324492379170,0.360761573048139,0.467913934572691,
216 6 0.467913934572691,0.360761573048139,0.171324492379170,
217 6 0. ,0. ,0. ,
218 7 0.129484966168870,0.279705391489277,0.381830050505119,
219 7 0.417959183673469,0.381830050505119,0.279705391489277,
220 7 0.129484966168870,0. ,0. ,
221 8 0.101228536290376,0.222381034453374,0.313706645877887,
222 8 0.362683783378362,0.362683783378362,0.313706645877887,
223 8 0.222381034453374,0.101228536290376,0. ,
224 9 0.081274388361574,0.180648160694857,0.260610696402935,
225 9 0.312347077040003,0.330239355001260,0.312347077040003,
226 9 0.260610696402935,0.180648160694857,0.081274388361574/
227C=======================================================================
228 igtyp=nint(geo(12,iprop))
229 igmat = igeo(98,iprop)
230 ipang = 1 ! NPT
231 ipthk = ipang + nlay ! NPT
232 ippos = ipthk + nlay ! NPT
233
234C
235 max_nptt = 1
236 IF(igtyp == 51 .OR. igtyp == 52) THEN
237 DO ilay=1,nlay ! NPT
238 max_nptt = max(max_nptt ,elbuf_str%BUFLY(ilay)%NPTT)
239 ENDDO
240 ALLOCATE(thk_it(max_nptt),pos_it(max_nptt))
241 ELSE
242 ALLOCATE(thk_it(0),pos_it(0))
243 ENDIF
244 nlay_max = max(nlay, npt)
245 ALLOCATE(posly(nlay_max*max_nptt),ratio_thkly(nlay_max*max_nptt))
246c-----------------------------
247 IF ((igtyp == 17 .OR. igtyp == 51) .AND. igmat < 0 ) THEN
248c-----------------------------
249 DO i=lft,llt
250 IF (ndrape == 0) THEN
251 IF (thke(i)== zero) THEN
252 thk(i) = stack%GEO(1,isubstack)
253 thke(i)= thk(i)
254 ELSE
255 thk(i)=thke(i)
256 ENDIF
257 thki(i) = stack%GEO(1,isubstack)
258 ELSEIF (ndrape > 0) THEN
259 thk(i) = thke(i)
260 thki(i)= thk(i)
261 ENDIF ! IF (NDRAPE == 0)
262c
263 IF (nperturb /= 0) THEN
264 DO j=1,nperturb
265 IF (perturb(j) == 1 .AND. rnoise(j,i+nft) /= zero) THEN
266 thk(i) = thk(i) * rnoise(j,i+nft)
267 thke(i) = thke(i) * rnoise(j,i+nft)
268 thki(i) = thki(i) * rnoise(j,i+nft)
269 ENDIF
270 ENDDO
271 ENDIF
272c
273 rho(i) = pm(89,imid)
274 rhocp(i) = pm(69,imid)
275 IF(thk(i) == zero.AND.igtyp/=0)THEN
276 CALL ancmsg(msgid=495,
277 . msgtype=msgerror,
278 . anmode=aninfo_blind_1,
279 . i1=ix(7,nft+i))
280 ENDIF
281 ENDDO
282c-----------------------------
283 ELSEIF (igtyp == 52 .OR.
284 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))THEN
285c-----------------------------
286 DO i=lft,llt
287 IF (ndrape == 0) THEN
288 IF (thke(i)== zero) THEN
289 thk(i) = stack%GEO(1,isubstack)
290 thke(i)= thk(i)
291 ELSE
292 thk(i)=thke(i)
293 ENDIF
294 thki(i) = stack%GEO(1,isubstack)
295 ELSEIF (ndrape > 0) THEN
296 thk(i) = thke(i)
297 thki(i)= thk(i)
298 ENDIF ! IF (NDRAPE == 0)
299 rho(i) = stack%PM(1,isubstack)
300 rhocp(i) = stack%PM(15,isubstack)* rho(i)/stack%PM(14,isubstack)
301c
302 IF (nperturb /= 0) THEN
303 DO j=1,nperturb
304 IF (perturb(j) == 1 .AND. rnoise(j,i+nft) /= zero) THEN
305 thk(i) = thk(i) * rnoise(j,i+nft)
306 thke(i) = thke(i) * rnoise(j,i+nft)
307 thki(i) = thki(i) * rnoise(j,i+nft)
308 ENDIF
309 ENDDO
310 ENDIF
311c
312 IF(thk(i) == zero.AND.igtyp/=0)THEN
313 CALL ancmsg(msgid=495,
314 . msgtype=msgerror,
315 . anmode=aninfo_blind_1,
316 . i1=ix(7,nft+i))
317 ENDIF
318 ENDDO
319c-----------------------------
320 ELSE ! IGTYP : std shell property
321c-----------------------------
322 DO i=lft,llt
323
324 IF (thke(i) == zero) THEN
325 thk(i) = geo(1,iprop)
326 thke(i)= thk(i)
327 ELSE
328 thk(i)=thke(i)
329 ENDIF
330c
331 thki(i) = thk(i)
332c
333 IF (nperturb /= 0) THEN
334 DO j=1,nperturb
335 IF (perturb(j) == 1 .AND. rnoise(j,i+nft) /= zero) THEN
336 thk(i) = thk(i) * rnoise(j,i+nft)
337 thke(i) = thke(i) * rnoise(j,i+nft)
338 thki(i) = thki(i) * rnoise(j,i+nft)
339 ENDIF
340 ENDDO
341 ENDIF
342c
343 rho(i) = pm(89,imid)
344 rhocp(i) = pm(69,imid)
345 IF(thk(i) == zero.AND.igtyp/=0)THEN
346 CALL ancmsg(msgid=495,
347 . msgtype=msgerror,
348 . anmode=aninfo_blind_1,
349 . i1=ix(7,nft+i))
350 ENDIF
351 ENDDO
352 ENDIF
353C----------------------------------------------
354C MASSES ELEMENTS / 4
355C----------------------------------------------
356 IF (igtyp == 11 .OR. igtyp == 16) THEN
357c--------------------------------------------
358 ipthk = 300
359 ipmat = 100
360 DO i=lft,llt
361 ems(i)=zero
362 DO n=1,npt
363 i1=ipthk+n
364 i2=ipmat+n
365 thkly = geo(i1,iprop)*thk(i)
366 matly = igeo(i2,iprop)
367 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth
368 ENDDO
369 mst = rho(i)*thk(i)*area(i)*fourth
370 dms = abs(ems(i)-mst)/mst
371C! for igty=11 ---> IGMAT=-1 (old G-mat), 1 for new g-mat
372C! for igtyp=16 --->IGMAT= 0
373C! for igtyp=17 --->IGMAT= 0
374 IF(igmat <= 0) THEN
375 IF (dms > em02) errm = 1
376 IF (dms > em02) THEN
377 idprop = igeo(1,iprop)
378 CALL fretitl2(titr,
379 . igeo(npropgi-ltitr+1,iprop),ltitr)
380 CALL ancmsg(msgid=519,
381 . msgtype=msgwarning,
382 . anmode=aninfo_blind_2,
383 . i1=idprop,
384 . c1=titr,i2=ix(nixc,nft+i),
385 . r1=ems(i),r2=mst)
386 ENDIF
387 ENDIF
388 ENDDO
389c--------------------------------------------
390 ELSEIF(igtyp == 17 ) THEN
391c--------------------------------------------
392 ippid = 2
393 ipmat = ippid + npt
394 ipang = 1
395 ipthk = ipang + npt
396 ippos = ipthk + npt
397 nthk = ippos + npt
398C
399 IF(idrape == 0) THEN
400 DO i=lft,llt
401 ems(i)=zero
402 mzi = zero
403 DO n=1,npt
404 i1 = ippid + n
405 i2 = ipmat + n
406 i3 = ipthk + n
407 i4 = ippos + n
408 thickt = stack%GEO(1 ,isubstack)
409 thkly = stack%GEO(i3 ,isubstack)*thk(i)
410 matly = stack%IGEO(i2,isubstack)
411 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth
412 ems_layer(i,n) = pm(89,matly)*thkly*area(i)*fourth
413 mzi = mzi + ems_layer(i,n)*stack%GEO(i4 ,isubstack)
414 ipid = stack%IGEO(i1, isubstack)
415 istack(i,n) = igeo(102,ipid)
416 ENDDO
417 mst = rho(i)*thk(i)*area(i)*fourth
418 dms = abs(ems(i)-mst)/mst
419 IF(dms > em02) errm = 1
420 IF(igmat <= 0) THEN
421 IF(dms > em02) THEN
422 idprop = igeo(1,iprop)
423 CALL fretitl2(titr,
424 . igeo(npropgi-ltitr+1,iprop),ltitr)
425 CALL ancmsg(msgid=519,
426 . msgtype=msgwarning,
427 . anmode=aninfo_blind_2,
428 . i1=idprop,
429 . c1=titr,i2=ix(nixc,nft+i),
430 . r1=ems(i),r2=mst)
431 ENDIF
432 ENDIF
433 !!
434 DO n=1,npt
435 i4 = ippos + n
436 IF(abs(stack%GEO(i4 ,isubstack)) < em15)
437 . stack%GEO(i4 ,isubstack)=zero
438 ENDDO
439 IF(abs(mzi) > em02) THEN
440 DO n=1,npt
441 i4 = ippos + n
442 dms = mzi/ems(i)
443 stack%GEO(i4 ,isubstack) = stack%GEO(i4 ,isubstack) - dms
444 ENDDO
445 ENDIF
446 ENDDO
447 ELSE ! idrape > 0
448 DO i=lft,llt
449 ems(i)=zero
450 mzi = zero
451 ie = indx(nft+i)
452 IF(ie == 0) THEN
453 DO n=1,npt
454 i1 = ippid + n
455 i2 = ipmat + n
456 i3 = ipthk + n
457 i4 = ippos + n
458 thickt = stack%GEO(1 ,isubstack)
459 thkly = stack%GEO(i3 ,isubstack)*thk(i)
460 matly = stack%IGEO(i2,isubstack)
461 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth
462 ems_layer(i,n) = pm(89,matly)*thkly*area(i)*fourth
463 mzi = mzi + ems_layer(i,n)*stack%GEO(i4 ,isubstack)
464 ipid = stack%IGEO(i1, isubstack)
465 istack(i,n) = igeo(102,ipid)
466 ENDDO
467 ELSE
468 thick_drape = drape(ie)%THICK
469 scale = thk(i)/thick_drape
470 DO n=1,npt
471 i1 = ippid + n
472 i2 = ipmat + n
473 i3 = ipthk + n
474 i4 = ippos + n
475 ip = drape(ie)%INDX_PLY(n)
476 IF(ip > 0 ) THEN
477 drape_ply => drape(ie)%DRAPE_PLY(ip)
478 thinning = drape_ply%RDRAPE(1,1) ! one slice by layer
479 thickt = stack%GEO(1 ,isubstack)
480 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
481 thkly = thkly*thinning
482 !! SCALE = THK(I)/THICK_DRAPE
483 matly = stack%IGEO(i2,isubstack)
484 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth*scale
485 ems_layer(i,n) = pm(89,matly)*thkly*area(i)*fourth *scale
486 ratio_thkly(n) = thkly/thick_drape
487 IF (n == 1) THEN
488 posly(n) = -half + half*ratio_thkly(n)
489 ELSE
490 posly(n) = posly(n-1)
491 . + half*(ratio_thkly(n)+ratio_thkly(n-1))
492 ENDIF ! IF (N == 1)
493
494 mzi = mzi + ems_layer(i,n)*posly(n)
495 ipid = stack%IGEO(i1, isubstack)
496 istack(i,n) = igeo(102,ipid)
497 ELSE ! IP=0
498 thickt = stack%GEO(1 ,isubstack)
499 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
500 matly = stack%IGEO(i2,isubstack)
501 !! SCALE = THK(I)/THICK_DRAPE
502 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth*scale
503 ems_layer(i,n) = pm(89,matly)*thkly*area(i)*fourth*scale
504 ratio_thkly(n) = thkly/thick_drape
505 IF (n == 1) THEN
506 posly(n) = -half + half*ratio_thkly(n)
507 ELSE
508 posly(n) = posly(n-1)
509 . + half*(ratio_thkly(n)+ratio_thkly(n-1))
510 ENDIF ! IF (N == 1)
511 mzi = mzi + ems_layer(i,n)*posly(n)
512 ipid = stack%IGEO(i1, isubstack)
513 istack(i,n) = igeo(102,ipid)
514 ENDIF
515 ENDDO
516 ENDIF ! ie
517 mst = rho(i)*thk(i)*area(i)*fourth
518 dms = abs(ems(i)-mst)/mst
519 IF(dms > em02) errm = 1
520 IF(igmat <= 0) THEN
521 IF(dms > em02) THEN
522 idprop = igeo(1,iprop)
523 CALL fretitl2(titr,
524 . igeo(npropgi-ltitr+1,iprop),ltitr)
525 CALL ancmsg(msgid=519,
526 . msgtype=msgwarning,
527 . anmode=aninfo_blind_2,
528 . i1=idprop,
529 . c1=titr,i2=ix(nixc,nft+i),
530 . r1=ems(i),r2=mst)
531 ENDIF
532 ENDIF
533C correction of the neutral surface position at the centroid
534 DO n=1,npt
535 i4 = ippos + n
536 IF(abs(stack%GEO(i4 ,isubstack)) < em15)
537 . stack%GEO(i4 ,isubstack)=zero
538 ENDDO
539 IF(abs(mzi) > em02) THEN
540 DO n=1,npt
541 i4 = ippos + n
542 dms = mzi/ems(i)
543 stack%GEO(i4 ,isubstack) = stack%GEO(i4 ,isubstack) - dms
544 ENDDO
545 ENDIF
546 ENDDO
547 ENDIF ! IDRAPE
548c--------------------------------------------
549 ELSEIF (igtyp == 51 .OR. igtyp == 52 ) THEN
550c--------------------------------------------
551
552 ipang = 1
553 ippid = 2
554 ipmat = ippid + nlay ! NPT
555 ipthk = ipang + nlay ! NPT
556 ippos = ipthk + nlay ! NPT
557 nthk = ippos + nlay ! NPT
558C
559 ipos = igeo(99,iprop)
560 thickt = stack%GEO(1,isubstack)
561 zshift = geo(199,iprop)
562 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
563 IF(idrape == 0 ) THEN
564 DO i=lft,llt
565 ems(i)=zero
566 mzi = zero
567 ipt_all = 0
568 scale = thk(i)/thickt
569 DO ilay=1,nlay ! NPT
570 nptt = elbuf_str%BUFLY(ilay)%NPTT
571 i1 = ippid + ilay
572 i2 = ipmat + ilay
573 i3 = ipthk + ilay
574 i4 = ippos + ilay
575 ipid = stack%IGEO(i1,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
576 iint = igeo(47,iprop)
577 thickt = stack%GEO(1,isubstack)
578 IF(iint == 1) THEN
579 thkly = stack%GEO(i3,isubstack)*thickt
580 laypos = stack%GEO(i4,isubstack)
581 ipid_ly = stack%IGEO(i1,isubstack)
582 matly = stack%IGEO(i2,isubstack)
583 DO it=1,nptt
584 ipt = ipt_all + it
585 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
586 ratio_thkly(ipt) = thk_it(it)/thickt
587 IF (ipt == 1) THEN
588 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
589 ELSE
590 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
591 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
592 ENDIF ! IF (IPT == 1)
593 pos_it(it) = posly(ipt)*thk(i)
594 ENDDO
595 ELSEIF(iint == 2) THEN
596 thkly = stack%GEO(i3,isubstack)*thickt
597 laypos = stack%GEO(i4,isubstack)
598 ipid_ly = stack%IGEO(i1,isubstack)
599 matly = stack%IGEO(i2,isubstack)
600 DO it=1,nptt
601 ipt = ipt_all + it
602 thk_it(it) = half*thkly*w_gauss(it,nptt) ! Gauss distribution of NPTT through layer
603 ratio_thkly(ipt) = thk_it(it)/thickt
604 IF (ipt == 1 ) THEN
605 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
606 ELSE
607 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
608 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
609 ENDIF ! IF (IPT == 1)
610 pos_it(it) = posly(ipt)*thk(i)
611 ENDDO
612 ENDIF ! iint
613 !!SCALE = THK(I)/THICKT
614 DO it=1,nptt
615 thk_it(it) = thk_it(it)*scale
616 ems_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
617 ems(i) = ems(i) + ems_nptt
618 mzi = mzi + ems_nptt * pos_it(it)/thk(i)
619 ENDDO
620 ipt_all = ipt_all + nptt
621 ENDDO ! NLAY
622 !!
623 mst = rho(i)*thk(i)*area(i)*fourth
624 dms = abs(ems(i)-mst)/mst
625 IF (dms > em02) errm = 1
626 IF(igtyp == 51 .AND. igmat < 0) THEN
627 IF (dms > em02) THEN
628 idprop = igeo(1,iprop)
629 CALL fretitl2(titr,
630 . igeo(npropgi-ltitr+1,iprop),ltitr)
631 CALL ancmsg(msgid=519,
632 . msgtype=msgwarning,
633 . anmode=aninfo_blind_2,
634 . i1=idprop,
635 . c1=titr,i2=ix(nixc,nft+i),
636 . r1=ems(i),r2=mst)
637 ENDIF
638 ENDIF
639 ENDDO ! LFT,LLT
640C------
641 IF( jthe > 0 ) THEN
642 DO i=lft,llt
643 emscp(i)=rhocp(i)*thk(i)*area(i)*fourth
644 ENDDO
645 ENDIF
646 ELSE ! idrape > 0
647 DO i=lft,llt
648 ems(i)=zero
649 mzi = zero
650 ipt_all = 0
651 ie = indx(nft+i)
652 IF(ie == 0) THEN
653 scale = thk(i)/thickt
654 DO ilay=1,nlay ! NPT
655 nptt = elbuf_str%BUFLY(ilay)%NPTT
656 i1 = ippid + ilay
657 i2 = ipmat + ilay
658 i3 = ipthk + ilay
659 i4 = ippos + ilay
660 ipid = stack%IGEO(i1,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
661 iint = igeo(47,iprop)
662 thickt = stack%GEO(1,isubstack)
663 IF(iint == 1) THEN
664 thkly = stack%GEO(i3,isubstack)*thickt
665 laypos = stack%GEO(i4,isubstack)
666 ipid_ly = stack%IGEO(i1,isubstack)
667 matly = stack%IGEO(i2,isubstack)
668 DO it=1,nptt
669 ipt = ipt_all + it
670 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
671 ratio_thkly(ipt) = thk_it(it)/thickt
672 IF (ipt == 1) THEN
673 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
674 ELSE
675 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
676 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
677 ENDIF ! IF (IPT == 1)
678 pos_it(it) = posly(ipt)*thk(i)
679 ENDDO
680 ELSEIF(iint == 2) THEN
681 thkly = stack%GEO(i3,isubstack)*thickt
682 laypos = stack%GEO(i4,isubstack)
683 ipid_ly = stack%IGEO(i1,isubstack)
684 matly = stack%IGEO(i2,isubstack)
685 DO it=1,nptt
686 ipt = ipt_all + it
687 thk_it(it) = half*thkly*w_gauss(it,nptt) ! Gauss distribution of NPTT through layer
688 ratio_thkly(ipt) = thk_it(it)/thickt
689 IF (ipt == 1 ) THEN
690 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
691 ELSE
692 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
693 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
694 ENDIF ! IF (IPT == 1)
695 pos_it(it) = posly(ipt)*thk(i)
696 ENDDO
697 ENDIF ! iint
698 !! SCALE = THK(I)/THICK
699 DO it=1,nptt
700 thk_it(it) = scale*thk_it(it)
701 ems_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
702 ems(i) = ems(i) + ems_nptt
703 mzi = mzi + ems_nptt * pos_it(it)/thk(i)
704 ENDDO
705 ipt_all = ipt_all + nptt
706 ENDDO ! NLAY
707 ELSE ! IE
708 thick_drape = drape(ie)%THICK
709 scale = thk(i)/thick_drape
710 DO ilay=1,nlay ! NPT
711 ip = drape(ie)%INDX_PLY(ilay)
712 nptt = elbuf_str%BUFLY(ilay)%NPTT ! = NSLICE
713 IF(ip > 0 ) THEN
714 drape_ply => drape(ie)%DRAPE_PLY(ip)
715 nslice = drape_ply%NSLICE
716 i1 = ippid + ilay
717 i2 = ipmat + ilay
718 i3 = ipthk + ilay
719 i4 = ippos + ilay
720 ipid = stack%IGEO(i1,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
721 iint = igeo(47,iprop)
722 thickt = stack%GEO(1,isubstack)
723 IF(iint == 1) THEN
724 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
725 ipid_ly = stack%IGEO(i1,isubstack)
726 matly = stack%IGEO(i2,isubstack)
727 DO it=1,nptt
728 ipt = ipt_all + it
729 thinning = drape_ply%RDRAPE(it,1)
730 thk_it(it) = thkly*thinning/nptt
731 ratio_thkly(ipt) = thk_it(it)/thick_drape
732 IF (ipt == 1) THEN
733 posly(ipt) = zshift + half*ratio_thkly(ipt)
734 ELSE
735 posly(ipt) = posly(ipt-1)
736 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
737 ENDIF ! IF (IPT == 1)
738 pos_it(it) = posly(ipt)*thk(i)
739 ENDDO ! nptt
740 ELSEIF(iint == 2) THEN
741 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
742 ipid_ly = stack%IGEO(i1,isubstack)
743 matly = stack%IGEO(i2,isubstack)
744 DO it=1,nptt
745 ipt = ipt_all + it
746 thinning = drape_ply%RDRAPE(it,1)
747 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
748 ratio_thkly(ipt) = thk_it(it)/thick_drape
749 IF (ipt == 1) THEN
750 posly(ipt) = zshift + half*ratio_thkly(ipt)
751 ELSE
752 posly(ipt) = posly(ipt-1)
753 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
754 ENDIF ! IF (IPT == 1)
755 pos_it(it) = posly(ipt)*thk(i)
756 ENDDO ! nptt
757 ENDIF ! iint
758 !! SCALE = THK(I)/THICK_DRAPE
759 DO it=1,nptt
760 thk_it(it)= scale*thk_it(it)
761 ems_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
762 ems(i) = ems(i) + ems_nptt
763 mzi = mzi + ems_nptt * pos_it(it)/thk(i)
764 ENDDO
765 ELSE ! Ip not draped
766 nptt = elbuf_str%BUFLY(ilay)%NPTT
767 i1 = ippid + ilay
768 i2 = ipmat + ilay
769 i3 = ipthk + ilay
770 i4 = ippos + ilay
771 ipid = stack%IGEO(i1,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
772 iint = igeo(47,iprop)
773 thickt = stack%GEO(1,isubstack)
774 IF(iint == 1) THEN
775 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
776 ipid_ly = stack%IGEO(i1,isubstack)
777 matly = stack%IGEO(i2,isubstack)
778 DO it=1,nptt
779 ipt = ipt_all + it
780 thk_it(it) = thkly/nptt
781 ratio_thkly(ipt) = thk_it(it)/thick_drape
782 IF (ipt == 1) THEN
783 posly(ipt) = zshift + half*ratio_thkly(ipt)
784 ELSE
785 posly(ipt) = posly(ipt-1)
786 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
787 ENDIF ! IF (IPT == 1)
788 pos_it(it) = posly(ipt)*thk(i)
789 ENDDO ! nptt
790 ELSEIF(iint == 2) THEN
791 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
792 ipid_ly = stack%IGEO(i1,isubstack)
793 matly = stack%IGEO(i2,isubstack)
794 DO it=1,nptt
795 ipt = ipt_all + it
796 thk_it(it) = half*thkly*w_gauss(it,nptt)
797 ratio_thkly(ipt) = thk_it(it)/thk(i)/thick_drape
798 IF (ipt == 1) THEN
799 posly(ipt) = zshift + half*ratio_thkly(ipt)
800 ELSE
801 posly(ipt) = posly(ipt-1)
802 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
803 ENDIF ! IF (IPT == 1)
804 pos_it(it) = posly(ipt)*thk(i)
805 ENDDO ! nptt
806 ENDIF ! iint
807 !! SCALE = THK(I)/THICK_DRAPE
808 DO it=1,nptt
809 thk_it(it ) = scale*thk_it(it)
810 ems_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
811 ems(i) = ems(i) + ems_nptt
812 mzi = mzi + ems_nptt * pos_it(it)/thk(i)
813 ENDDO
814 ENDIF ! IP
815 ipt_all = ipt_all + nptt
816 ENDDO ! NLAY
817 ENDIF ! ie
818
819 mst = rho(i)*thk(i)*area(i)*fourth
820 dms = abs(ems(i)-mst)/mst
821 IF (dms > em02) errm = 1
822 IF(igtyp == 51 .AND. igmat < 0) THEN
823 IF (dms > em02) THEN
824 idprop = igeo(1,iprop)
825 CALL fretitl2(titr,
826 . igeo(npropgi-ltitr+1,iprop),ltitr)
827 CALL ancmsg(msgid=519,
828 . msgtype=msgwarning,
829 . anmode=aninfo_blind_2,
830 . i1=idprop,
831 . c1=titr,i2=ix(nixc,nft+i),
832 . r1=ems(i),r2=mst)
833 ENDIF
834 ENDIF
835C
836 ENDDO
837C------
838 IF( jthe > 0 ) THEN
839 DO i=lft,llt
840 emscp(i)=rhocp(i)*thk(i)*area(i)*fourth
841 ENDDO
842 ENDIF
843 ENDIF ! idrape
844c--------------------------------------------
845 ELSE ! IGTYP
846c--------------------------------------------
847 DO i=lft,llt
848 ems(i)=rho(i)*thk(i)*area(i)*fourth
849 ENDDO
850C
851 IF( jthe > 0 ) THEN
852 DO i=lft,llt
853 emscp(i) = rhocp(i)*thk(i)*area(i)*fourth
854 ENDDO
855 ENDIF
856C
857 ENDIF ! IGTYP
858C----------------------------------------------
859C initialization of nodal masses
860c elt mass (parith/on spmd + adaptive meshing)
861C----------------------------------------------
862 IF(jthe == 0 ) THEN
863 DO i=lft,llt
864 msc(i)=ems(i)
865 ENDDO
866 ELSE
867 DO i=lft,llt
868 msc(i)=ems(i)
869 mcps(i) = emscp(i)
870 ENDDO
871 ENDIF
872c
873 IF(ishxfem_ply > 0 )THEN
874 ippid = 2
875 ipmat = ippid + npt
876 ipang = 1
877 ipthk = ipang + npt
878 ippos = ipthk + npt
879 nthk = ippos + npt
880 DO i=lft,llt
881 ii = nft + i
882 DO j=1,npt
883 ip = istack(i,j)
884 ms_layerc(ii,ip)= ems_layer(i,j)
885C
886 i4 = ippos + j
887 thickt = stack%GEO(1 ,isubstack)
888 islv = igeo(32,iprop)
889 IF(islv == 0)THEN
890 zi_layerc(ii,ip) = thickt*stack%GEO(i4 ,isubstack)
891 zply(ip) = zi_layerc(ii,ip)
892 ENDIF
893cc ZJ = ZI_LAYERC(II,IP)
894 zj = thickt*stack%GEO(i4 ,isubstack)
895 msz2c(ii) = msz2c(ii) + ems_layer(i,j)*zj*zj
896 ENDDO
897 ENDDO
898 ENDIF
899C
900 IF (jthe > 0 .or. nintemp > 0) THEN
901 IF (nintemp > 0 ) THEN
902 DO i= lft,llt
903 IF(temp(ix1(i))== zero) temp(ix1(i)) = pm(79,imid)
904 IF(temp(ix2(i))== zero) temp(ix2(i)) = pm(79,imid)
905 IF(temp(ix3(i))== zero) temp(ix3(i)) = pm(79,imid)
906 IF(temp(ix4(i))== zero) temp(ix4(i)) = pm(79,imid)
907 ENDDO
908 ELSE
909 DO i=lft,llt
910 temp(ix1(i))= pm(79,imid)
911 temp(ix2(i))= pm(79,imid)
912 temp(ix3(i))= pm(79,imid)
913 temp(ix4(i))= pm(79,imid)
914 ENDDO
915 ENDIF
916 ENDIF
917C----------------------------------------------
918C INERTIES ELEMENTS /4
919C----------------------------------------------
920 IF(iner_9_12 /= zero)THEN
921 fac=iner_9_12
922 ELSEIF(ihbe>=11)THEN
923 fac=twelve
924 ELSE
925 fac=nine
926 ENDIF
927 IF (igtyp == 11 .OR. igtyp == 16) THEN
928 ipthk = 300
929 ippos = 400
930 DO i=lft,llt
931 xi(i)=zero
932 DO n=1,npt
933 i1=ipthk+n
934 i3=ippos+n
935 thickt= geo(200,iprop)
936 thkly = geo(i1,iprop)*thk(i)
937 laypos = geo(i3,iprop)*thk(i)
938 i2=ipmat+n
939 matly = igeo(i2,iprop)
940 msl = pm(89,matly)*thkly*area(i)*fourth
941 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
942 xi(i) = xi(i) + xil
943 ENDDO
944 ENDDO
945 ELSEIF(igtyp == 17) THEN
946 ippid = 2
947 ipmat = ippid + npt
948 ipang = 1
949 ipthk = ipang + npt
950 ippos = ipthk + npt
951 nthk = ippos + npt
952 ipos = igeo(99,iprop)
953 thickt = stack%GEO(1,isubstack)
954 zshift = geo(199,iprop)
955 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
956C
957 IF(idrape == 0) THEN
958 DO i=lft,llt
959 xi(i)=zero
960 DO n=1,npt
961 i1 = ippid + n
962 i2 = ipmat + n
963 i3 = ipthk + n
964 i4 = ippos + n
965 thickt= stack%GEO(1,isubstack)
966 thkly = stack%GEO(i3 ,isubstack )*thk(i)
967 laypos = stack%GEO(i4,isubstack)*thk(i)
968 matly = stack%IGEO(i2 ,isubstack)
969 msl = pm(89,matly)*thkly*area(i)*fourth
970 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
971 xi(i) = xi(i) + xil
972 ENDDO
973 ENDDO
974 ELSE ! IDRAPE > 0
975 DO i=lft,llt
976 xi(i)=zero
977 ie = indx(nft + i)
978 IF(ie == 0) THEN
979 DO n=1,npt
980 i1 = ippid + n
981 i2 = ipmat + n
982 i3 = ipthk + n
983 i4 = ippos + n
984 thickt= stack%GEO(1,isubstack)
985 thkly = stack%GEO(i3 ,isubstack )*thk(i)
986 laypos = stack%GEO(i4,isubstack)*thk(i)
987 matly = stack%IGEO(i2 ,isubstack)
988 msl = pm(89,matly)*thkly*area(i)*fourth
989 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
990 xi(i) = xi(i) + xil
991 ENDDO
992 ELSE
993 thick_drape = drape(ie)%THICK
994 scale = thk(i)/thick_drape
995 DO n=1,npt
996 ip = drape(ie)%INDX_PLY(n)
997 IF(ip > 0) THEN
998 drape_ply => drape(ie)%DRAPE_PLY(ip)
999 i1 = ippid + n
1000 i2 = ipmat + n
1001 i3 = ipthk + n
1002 i4 = ippos + n
1003 thinning = drape_ply%RDRAPE(1,1)
1004 thickt= stack%GEO(1,isubstack)
1005 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1006 thkly = thkly*thinning
1007 matly = stack%IGEO(i2,isubstack)
1008 ratio_thkly(n) = thkly/thick_drape
1009 IF (n == 1) THEN
1010 posly(n) = zshift + half*ratio_thkly(n)
1011 ELSE
1012 posly(n) = posly(n-1)
1013 . + half*(ratio_thkly(n)+ratio_thkly(n-1))
1014 ENDIF ! IF (N == 1)
1015 laypos = posly(n)*thk(i)
1016 thkly = thkly*scale
1017 msl = pm(89,matly)*thkly*area(i)*fourth
1018 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
1019 xi(i) = xi(i) + xil
1020 ELSE !IP ==0
1021 i1 = ippid + n
1022 i2 = ipmat + n
1023 i3 = ipthk + n
1024 i4 = ippos + n
1025 thickt= stack%GEO(1,isubstack)
1026 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1027 matly = stack%IGEO(i2,isubstack)
1028 ratio_thkly(n) = thkly/thick_drape
1029 IF (n == 1) THEN
1030 posly(n) = zshift + half*ratio_thkly(n)
1031 ELSE
1032 posly(n) = posly(n-1)
1033 . + half*(ratio_thkly(n)+ratio_thkly(n-1))
1034 ENDIF ! IF (N == 1)
1035 laypos = posly(n)*thk(i)
1036 thkly = thkly*scale
1037 msl = pm(89,matly)*thkly*area(i)*fourth
1038 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
1039 xi(i) = xi(i) + xil
1040 ENDIF
1041 ENDDO
1042 ENDIF
1043 ENDDO
1044 ENDIF
1045 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
1046 ipang = 1
1047 ippid = 2
1048 ipmat = ippid + nlay ! NPT
1049 ipthk = ipang + nlay ! NPT
1050 ippos = ipthk + nlay ! NPT
1051 nthk = ippos + nlay ! NPT
1052 ipos = igeo(99,iprop)
1053 thickt = stack%GEO(1,isubstack)
1054 zshift = geo(199,iprop)
1055 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
1056C
1057 IF(idrape == 0) THEN
1058 DO i=lft,llt
1059 scale = thk(i)/thickt
1060 xi(i)=zero
1061 ipt_all = 0
1062 DO ilay=1, nlay ! NPT
1063 nptt = elbuf_str%BUFLY(ilay)%NPTT
1064 i1 = ippid + ilay
1065 i2 = ipmat + ilay
1066 i3 = ipthk + ilay
1067 i4 = ippos + ilay
1068 ipid_ly = stack%IGEO(i1,isubstack)
1069 matly = stack%IGEO(i2,isubstack)
1070 iint = igeo(47,iprop)
1071 IF(iint == 1) THEN
1072 thickt = stack%GEO(1,isubstack)
1073 thkly = stack%GEO(i3,isubstack)*thickt
1074 laypos = stack%GEO(i4 ,isubstack)*thickt
1075 DO it=1,nptt
1076 ipt = ipt_all + it
1077 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
1078 ratio_thkly(ipt) = thk_it(it)/thickt
1079 IF (ipt == 1 ) THEN
1080 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
1081 ELSE
1082 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
1083 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
1084 ENDIF ! IF (IPT == 1)
1085 pos_it(it) = posly(ipt)*thk(i)
1086 ENDDO
1087 ELSEIF(iint == 2) THEN
1088 thickt = stack%GEO(1,isubstack)
1089 thkly = stack%GEO(i3,isubstack)*thickt
1090 laypos = stack%GEO(i4 ,isubstack)*thickt
1091 DO it=1,nptt
1092 ipt = ipt_all + it
1093 thk_it(it) = half*thkly*w_gauss(it,nptt)
1094 ratio_thkly(ipt) = thk_it(it)/thickt
1095 IF (ipt == 1) THEN
1096 posly(ipt) = zshift + half*ratio_thkly(ipt)
1097 ELSE
1098 posly(ipt) = posly(ipt-1)
1099 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1100 ENDIF ! IF (IPT 1)
1101 pos_it(it) = posly(ipt)*thk(i)
1102 ENDDO ! nptt
1103 ENDIF ! iint
1104C
1105 !!SCALE = THK(I)/THICKT
1106 DO it=1,nptt
1107 thk_it(it) = scale*thk_it(it)
1108 ms_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
1109 xil = ms_nptt*(area(i)/fac+thk_it(it)*thk_it(it)*one_over_12+
1110 . pos_it(it)*pos_it(it))
1111 xi(i) = xi(i) + xil
1112 ENDDO
1113 ipt_all = ipt_all + nptt
1114 ENDDO ! NLAY
1115 ENDDO ! LFT,LLT
1116C
1117 ELSE ! IDRAPE > 0
1118 DO i=lft,llt
1119 xi(i)=zero
1120 ipt_all = 0
1121 ie = indx(nft + i)
1122 IF(ie == 0) THEN
1123 scale = thk(i)/thickt
1124 DO ilay=1, nlay ! NPT
1125 nptt = elbuf_str%BUFLY(ilay)%NPTT
1126 i1 = ippid + ilay
1127 i2 = ipmat + ilay
1128 i3 = ipthk + ilay
1129 i4 = ippos + ilay
1130 ipid_ly = stack%IGEO(i1,isubstack)
1131 matly = stack%IGEO(i2,isubstack)
1132 iint = igeo(47,iprop)
1133 IF(iint == 1) THEN
1134 thickt = stack%GEO(1,isubstack)
1135 thkly = stack%GEO(i3,isubstack)*thickt
1136 laypos = stack%GEO(i4 ,isubstack)*thickt
1137 DO it=1,nptt
1138 ipt = ipt_all + it
1139 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
1140 ratio_thkly(ipt) = thk_it(it)/thickt
1141 IF (ipt == 1 ) THEN
1142 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
1143 ELSE
1144 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
1145 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
1146 ENDIF ! IF (IPT == 1)
1147 pos_it(it) = posly(ipt)*thk(i)
1148 ENDDO
1149 ELSEIF(iint == 2) THEN
1150 thickt = stack%GEO(1,isubstack)
1151 thkly = stack%GEO(i3,isubstack)*thickt
1152 laypos = stack%GEO(i4 ,isubstack)*thickt
1153 DO it=1,nptt
1154 ipt = ipt_all + it
1155 thk_it(it) = half*thkly*w_gauss(it,nptt)
1156 ratio_thkly(ipt) = thk_it(it)/thickt
1157 IF (ipt == 1) THEN
1158 posly(ipt) = zshift + half*ratio_thkly(ipt)
1159 ELSE
1160 posly(ipt) = posly(ipt-1)
1161 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1162 ENDIF ! IF (IPT 1)
1163 pos_it(it) = posly(ipt)*thk(i)
1164 ENDDO ! nptt
1165 ENDIF ! iint
1166C
1167 !!SCALE = THK(I)/THICKT
1168 DO it=1,nptt
1169 thk_it(it) = scale*thk_it(it)
1170 ms_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
1171 xil = ms_nptt*(area(i)/fac+thk_it(it)*thk_it(it)*one_over_12+
1172 . pos_it(it)*pos_it(it))
1173 xi(i) = xi(i) + xil
1174 ENDDO
1175 ipt_all = ipt_all + nptt
1176 ENDDO ! NLAY
1177 ELSE ! IE
1178 thick_drape = drape(ie)%THICK
1179 scale = thk(i)/thick_drape
1180 DO ilay=1, nlay ! NPT
1181 nptt = elbuf_str%BUFLY(ilay)%NPTT
1182 ip = drape(ie)%INDX_PLY(ilay)
1183 IF(ip > 0) THEN
1184 drape_ply => drape(ie)%DRAPE_PLY(ip)
1185 nslice = drape_ply%NSLICE ! NPTT
1186 i1 = ippid + ilay
1187 i2 = ipmat + ilay
1188 i3 = ipthk + ilay
1189 i4 = ippos + ilay
1190 ipid_ly = stack%IGEO(i1,isubstack)
1191 matly = stack%IGEO(i2,isubstack)
1192 iint = igeo(47,iprop)
1193 IF(iint == 1) THEN
1194 thickt = stack%GEO(1,isubstack)
1195 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1196 ipid_ly = stack%IGEO(i1,isubstack)
1197 matly = stack%IGEO(i2,isubstack)
1198 DO it=1,nptt
1199 ipt = ipt_all + it
1200 thinning = drape_ply%RDRAPE(it,1)
1201 thk_it(it) = thkly*thinning/nptt
1202 ratio_thkly(ipt) = thk_it(it)/thick_drape
1203 IF (ipt == 1) THEN
1204 posly(ipt) = zshift + half*ratio_thkly(ipt)
1205 ELSE
1206 posly(ipt) = posly(ipt-1)
1207 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1208 ENDIF ! IF (IPT == 1)
1209 pos_it(it) = posly(ipt)*thk(i)
1210 ENDDO ! nptt
1211 ELSEIF(iint == 2) THEN
1212 thickt = stack%GEO(1,isubstack)
1213 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1214 DO it=1,nptt
1215 ipt = ipt_all + it
1216 thinning = drape_ply%RDRAPE(it,1)
1217 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
1218 ratio_thkly(ipt) = thk_it(it)/thick_drape
1219 IF (ipt == 1) THEN
1220 posly(ipt) = zshift + half*ratio_thkly(ipt)
1221 ELSE
1222 posly(ipt) = posly(ipt-1)
1223 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1224 ENDIF ! IF (IPT 1)
1225 pos_it(it) = posly(ipt)*thk(i)
1226 ENDDO ! nptt
1227 ENDIF ! iint
1228 !!SCALE = THK(I)/THICK_DRAPE
1229 DO it=1,nptt
1230 thk_it(it) = scale * thk_it(it)
1231 ms_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
1232 xil = ms_nptt*(area(i)/fac+thk_it(it)*thk_it(it)*one_over_12+
1233 . pos_it(it)*pos_it(it))
1234 xi(i) = xi(i) + xil
1235 ENDDO
1236 ELSE ! IP == 0
1237 nptt = elbuf_str%BUFLY(ilay)%NPTT
1238 i1 = ippid + ilay
1239 i2 = ipmat + ilay
1240 i3 = ipthk + ilay
1241 i4 = ippos + ilay
1242 ipid_ly = stack%IGEO(i1,isubstack)
1243 matly = stack%IGEO(i2,isubstack)
1244 iint = igeo(47,iprop)
1245 IF(iint == 1) THEN
1246 thickt = stack%GEO(1,isubstack)
1247 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1248 ipid_ly = stack%IGEO(i1,isubstack)
1249 matly = stack%IGEO(i2,isubstack)
1250 DO it=1,nptt
1251 ipt = ipt_all + it
1252 thk_it(it) = thkly/nptt
1253 ratio_thkly(ipt) = thk_it(it)/thick_drape
1254 IF (ipt == 1) THEN
1255 posly(ipt) = zshift + half*ratio_thkly(ipt)
1256 ELSE
1257 posly(ipt) = posly(ipt-1)
1258 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1259 ENDIF ! IF (IPT == 1)
1260 pos_it(it) = posly(ipt)*thk(i)
1261 ENDDO ! nptt
1262 ELSEIF(iint == 2) THEN
1263 thickt = stack%GEO(1,isubstack)
1264 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1265 DO it=1,nptt
1266 ipt = ipt_all + it
1267 thk_it(it) = half*thkly*w_gauss(it,nptt)
1268 ratio_thkly(ipt) = thk_it(it)/thick_drape
1269 IF (ipt == 1) THEN
1270 posly(ipt) = zshift + half*ratio_thkly(ipt)
1271 ELSE
1272 posly(ipt) = posly(ipt-1)
1273 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1274 ENDIF ! IF (IPT 1)
1275 pos_it(it) = posly(ipt)*thk(i)
1276 ENDDO ! nptt
1277 ENDIF ! iint
1278 !!SCALE = THK(I)/THICK_DRAPE
1279 DO it=1,nptt
1280 thk_it(it) = scale*thk_it(it)
1281 ms_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
1282 xil = ms_nptt*(area(i)/fac+thk_it(it)*thk_it(it)*one_over_12+
1283 . pos_it(it)*pos_it(it))
1284 xi(i) = xi(i) + xil
1285 ENDDO
1286 ENDIF ! IP
1287 ipt_all = ipt_all + nptt
1288 ENDDO ! NLAY
1289 ENDIF ! IE
1290 ENDDO ! LFT,LLT
1291 ENDIF ! idrape
1292 ELSE
1293 zshift = geo(199,iprop)
1294 f_offset= zshift*zshift
1295 DO i=lft,llt
1296 xi(i)=ems(i)*(area(i)/fac+thk(i)*thk(i)*(one_over_12+f_offset))
1297 ENDDO
1298 ENDIF
1299C--- check stability condition
1300 IF (igtyp == 11 .OR. igtyp == 16 .OR. igtyp == 17
1301 . .OR. igtyp == 51)THEN
1302 IF(igmat <= 0) THEN !
1303 DO i=lft,llt
1304 ins = ems(i)*(area(i)/fac+thk(i)*thk(i)*one_over_12)
1305 dms = abs(xi(i)-ins)/ins
1306 IF(dms > em02) THEN
1307 idprop = igeo(1,iprop)
1308 CALL fretitl2(titr,
1309 . igeo(npropgi-ltitr+1,iprop),ltitr)
1310 CALL ancmsg(msgid=520,
1311 . msgtype=msgwarning,
1312 . anmode=aninfo_blind_2,
1313 . i1=idprop,
1314 . c1=titr,i2=ix(nixc,nft+i))
1315 ENDIF
1316 ENDDO
1317 ENDIF
1318 ENDIF
1319C----------------------------------------------
1320C-----------------------
1321 DO i=lft,llt
1322 inc(i)=xi(i)
1323 ENDDO
1324C-----------------------
1325 IF(nadmesh==0)THEN
1326 IF (nxref == 0) THEN
1327 DO i=lft,llt
1328 i1 = ix1(i)
1329 i2 = ix2(i)
1330 i3 = ix3(i)
1331 i4 = ix4(i)
1332 ip = ipart(i)
1333C
1334
1335 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1336 partsav(2,ip)=partsav(2,ip) + ems(i)*
1337 . (x(1,i1)+x(1,i2)+x(1,i3)+x(1,i4))
1338 partsav(3,ip)=partsav(3,ip) + ems(i)*
1339 . (x(2,i1)+x(2,i2)+x(2,i3)+x(2,i4))
1340 partsav(4,ip)=partsav(4,ip) + ems(i)*
1341 . (x(3,i1)+x(3,i2)+x(3,i3)+x(3,i4))
1342 xx = (x(1,i1)*x(1,i1)+x(1,i2)*x(1,i2)
1343 . +x(1,i3)*x(1,i3)+x(1,i4)*x(1,i4))
1344 xy = (x(1,i1)*x(2,i1)+x(1,i2)*x(2,i2)
1345 . +x(1,i3)*x(2,i3)+x(1,i4)*x(2,i4))
1346 yy = (x(2,i1)*x(2,i1)+x(2,i2)*x(2,i2)
1347 . +x(2,i3)*x(2,i3)+x(2,i4)*x(2,i4))
1348 yz = (x(2,i1)*x(3,i1)+x(2,i2)*x(3,i2)
1349 . +x(2,i3)*x(3,i3)+x(2,i4)*x(3,i4))
1350 zz = (x(3,i1)*x(3,i1)+x(3,i2)*x(3,i2)
1351 . +x(3,i3)*x(3,i3)+x(3,i4)*x(3,i4))
1352 zx = (x(3,i1)*x(1,i1)+x(3,i2)*x(1,i2)
1353 . +x(3,i3)*x(1,i3)+x(3,i4)*x(1,i4))
1354C
1355 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i) * (yy+zz)
1356 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i) * (zz+xx)
1357 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i) * (xx+yy)
1358 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1359 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1360 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1361C
1362 partsav(11,ip)=partsav(11,ip) + ems(i)*
1363 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1364 partsav(12,ip)=partsav(12,ip) + ems(i)*
1365 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1366 partsav(13,ip)=partsav(13,ip) + ems(i)*
1367 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1368C
1369 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1370 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1371 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1372 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1373 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1374 END DO
1375 ELSE
1376 DO i=lft,llt
1377 i1 = ix1(i)
1378 i2 = ix2(i)
1379 i3 = ix3(i)
1380 i4 = ix4(i)
1381 ip = ipart(i)
1382C
1383 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1384 partsav(2,ip)=partsav(2,ip) + ems(i)*
1385 . (xrefc(1,1,i)+xrefc(2,1,i)+xrefc(3,1,i)+xrefc(4,1,i))
1386 partsav(3,ip)=partsav(3,ip) + ems(i)*
1387 . (xrefc(1,2,i)+xrefc(2,2,i)+xrefc(3,2,i)+xrefc(4,2,i))
1388 partsav(4,ip)=partsav(4,ip) + ems(i)*
1389 . (xrefc(1,3,i)+xrefc(2,3,i)+xrefc(3,3,i)+xrefc(4,3,i))
1390 xx = xrefc(1,1,i)**2 + xrefc(2,1,i)**2
1391 . + xrefc(3,1,i)**2 + xrefc(4,1,i)
1392 xy = xrefc(1,1,i)*xrefc(1,2,i) + xrefc(2,1,i)*xrefc(2,2,i)
1393 . + xrefc(3,1,i)*xrefc(3,2,i) + xrefc(4,1,i)*xrefc(4,2,i)
1394 yy = xrefc(1,2,i)**2 + xrefc(2,2,i)**2
1395 . + xrefc(3,2,i)**2 + xrefc(4,2,i)
1396 yz = xrefc(1,2,i)*xrefc(1,3,i) + xrefc(2,2,i)*xrefc(2,3,i)
1397 . + xrefc(3,2,i)*xrefc(3,3,i) + xrefc(4,2,i)*xrefc(4,3,i)
1398 zz = xrefc(1,3,i)**2 + xrefc(2,3,i)**2
1399 . + xrefc(3,3,i)**2 + xrefc(4,3,i)
1400 zx = xrefc(1,3,i)*xrefc(1,1,i) + xrefc(2,3,i)*xrefc(2,1,i)
1401 . + xrefc(3,3,i)*xrefc(3,1,i) + xrefc(4,3,i)*xrefc(4,1,i)
1402C
1403 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1404 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1405 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1406 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1407 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1408 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1409C
1410 partsav(11,ip)=partsav(11,ip) + ems(i)*
1411 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1412 partsav(12,ip)=partsav(12,ip) + ems(i)*
1413 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1414 partsav(13,ip)=partsav(13,ip) + ems(i)*
1415 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1416C
1417 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1418 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1419 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1420 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1421 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1422 END DO
1423 ENDIF
1424 ELSEIF(istatcnd==0) THEN
1425 IF (nxref == 0) THEN
1426 DO i=lft,llt
1427 n=nft+i
1428 IF(sh4tree(3,n) >= 0)THEN
1429 i1 = ix1(i)
1430 i2 = ix2(i)
1431 i3 = ix3(i)
1432 i4 = ix4(i)
1433 ip = ipart(i)
1434C
1435 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1436 partsav(2,ip)=partsav(2,ip) + ems(i)*
1437 . (x(1,i1)+x(1,i2)+x(1,i3)+x(1,i4))
1438 partsav(3,ip)=partsav(3,ip) + ems(i)*
1439 . (x(2,i1)+x(2,i2)+x(2,i3)+x(2,i4))
1440 partsav(4,ip)=partsav(4,ip) + ems(i)*
1441 . (x(3,i1)+x(3,i2)+x(3,i3)+x(3,i4))
1442 xx = (x(1,i1)*x(1,i1)+x(1,i2)*x(1,i2)
1443 . +x(1,i3)*x(1,i3)+x(1,i4)*x(1,i4))
1444 xy = (x(1,i1)*x(2,i1)+x(1,i2)*x(2,i2)
1445 . +x(1,i3)*x(2,i3)+x(1,i4)*x(2,i4))
1446 yy = (x(2,i1)*x(2,i1)+x(2,i2)*x(2,i2)
1447 . +x(2,i3)*x(2,i3)+x(2,i4)*x(2,i4))
1448 yz = (x(2,i1)*x(3,i1)+x(2,i2)*x(3,i2)
1449 . +x(2,i3)*x(3,i3)+x(2,i4)*x(3,i4))
1450 zz = (x(3,i1)*x(3,i1)+x(3,i2)*x(3,i2)
1451 . +x(3,i3)*x(3,i3)+x(3,i4)*x(3,i4))
1452 zx = (x(3,i1)*x(1,i1)+x(3,i2)*x(1,i2)
1453 . +x(3,i3)*x(1,i3)+x(3,i4)*x(1,i4))
1454C
1455 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1456 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1457 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1458 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1459 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1460 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1461C
1462 partsav(11,ip)=partsav(11,ip) + ems(i)*
1463 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1464 partsav(12,ip)=partsav(12,ip) + ems(i)*
1465 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1466 partsav(13,ip)=partsav(13,ip) + ems(i)*
1467 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1468C
1469 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1470 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1471 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1472 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1473 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1474 END IF
1475 END DO
1476 ELSE
1477 DO i=lft,llt
1478 n=nft+i
1479 IF(sh4tree(3,n) >= 0)THEN
1480 i1 = ix1(i)
1481 i2 = ix2(i)
1482 i3 = ix3(i)
1483 i4 = ix4(i)
1484 ip = ipart(i)
1485C
1486 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1487 partsav(2,ip)=partsav(2,ip) + ems(i)*
1488 . (xrefc(1,1,i)+xrefc(2,1,i)+xrefc(3,1,i)+xrefc(4,1,i))
1489 partsav(3,ip)=partsav(3,ip) + ems(i)*
1490 . (xrefc(1,2,i)+xrefc(2,2,i)+xrefc(3,2,i)+xrefc(4,2,i))
1491 partsav(4,ip)=partsav(4,ip) + ems(i)*
1492 . (xrefc(1,3,i)+xrefc(2,3,i)+xrefc(3,3,i)+xrefc(4,3,i))
1493 xx = xrefc(1,1,i)**2 + xrefc(2,1,i)**2
1494 . + xrefc(3,1,i)**2 + xrefc(4,1,i)
1495 xy = xrefc(1,1,i)*xrefc(1,2,i) + xrefc(2,1,i)*xrefc(2,2,i)
1496 . + xrefc(3,1,i)*xrefc(3,2,i) + xrefc(4,1,i)*xrefc(4,2,i)
1497 yy = xrefc(1,2,i)**2 + xrefc(2,2,i)**2
1498 . + xrefc(3,2,i)**2 + xrefc(4,2,i)
1499 yz = xrefc(1,2,i)*xrefc(1,3,i) + xrefc(2,2,i)*xrefc(2,3,i)
1500 . + xrefc(3,2,i)*xrefc(3,3,i) + xrefc(4,2,i)*xrefc(4,3,i)
1501 zz = xrefc(1,3,i)**2 + xrefc(2,3,i)**2
1502 . + xrefc(3,3,i)**2 + xrefc(4,3,i)
1503 zx = xrefc(1,3,i)*xrefc(1,1,i) + xrefc(2,3,i)*xrefc(2,1,i)
1504 . + xrefc(3,3,i)*xrefc(3,1,i) + xrefc(4,3,i)*xrefc(4,1,i)
1505C
1506 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1507 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1508 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1509 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1510 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1511 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1512C
1513 partsav(11,ip)=partsav(11,ip) + ems(i)*
1514 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1515 partsav(12,ip)=partsav(12,ip) + ems(i)*
1516 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1517 partsav(13,ip)=partsav(13,ip) + ems(i)*
1518 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1519C
1520 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1521 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1522 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1523 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1524 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1525 END IF
1526 END DO
1527 ENDIF
1528 ELSE
1529 IF (nxref == 0) THEN
1530 DO i=lft,llt
1531 n=nft+i
1532 IF(sh4tree(3,n) == 0 .OR. sh4tree(3,n) == -1)THEN
1533 i1 = ix1(i)
1534 i2 = ix2(i)
1535 i3 = ix3(i)
1536 i4 = ix4(i)
1537 ip = ipart(i)
1538C
1539 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1540 partsav(2,ip)=partsav(2,ip) + ems(i)*
1541 . (x(1,i1)+x(1,i2)+x(1,i3)+x(1,i4))
1542 partsav(3,ip)=partsav(3,ip) + ems(i)*
1543 . (x(2,i1)+x(2,i2)+x(2,i3)+x(2,i4))
1544 partsav(4,ip)=partsav(4,ip) + ems(i)*
1545 . (x(3,i1)+x(3,i2)+x(3,i3)+x(3,i4))
1546 xx = (x(1,i1)*x(1,i1)+x(1,i2)*x(1,i2)
1547 . +x(1,i3)*x(1,i3)+x(1,i4)*x(1,i4))
1548 xy = (x(1,i1)*x(2,i1)+x(1,i2)*x(2,i2)
1549 . +x(1,i3)*x(2,i3)+x(1,i4)*x(2,i4))
1550 yy = (x(2,i1)*x(2,i1)+x(2,i2)*x(2,i2)
1551 . +x(2,i3)*x(2,i3)+x(2,i4)*x(2,i4))
1552 yz = (x(2,i1)*x(3,i1)+x(2,i2)*x(3,i2)
1553 . +x(2,i3)*x(3,i3)+x(2,i4)*x(3,i4))
1554 zz = (x(3,i1)*x(3,i1)+x(3,i2)*x(3,i2)
1555 . +x(3,i3)*x(3,i3)+x(3,i4)*x(3,i4))
1556 zx = (x(3,i1)*x(1,i1)+x(3,i2)*x(1,i2)
1557 . +x(3,i3)*x(1,i3)+x(3,i4)*x(1,i4))
1558C
1559 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1560 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1561 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1562 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1563 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1564 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1565C
1566 partsav(11,ip)=partsav(11,ip) + ems(i)*
1567 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1568 partsav(12,ip)=partsav(12,ip) + ems(i)*
1569 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1570 partsav(13,ip)=partsav(13,ip) + ems(i)*
1571 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1572C
1573 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1574 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1575 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1576 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1577 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1578 END IF
1579 END DO
1580 ELSE
1581 DO i=lft,llt
1582 n=nft+i
1583 IF(sh4tree(3,n) == 0 .OR. sh4tree(3,n) == -1)THEN
1584 i1 = ix1(i)
1585 i2 = ix2(i)
1586 i3 = ix3(i)
1587 i4 = ix4(i)
1588 ip = ipart(i)
1589C
1590 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1591 partsav(2,ip)=partsav(2,ip) + ems(i)*
1592 . (xrefc(1,1,i)+xrefc(2,1,i)+xrefc(3,1,i)+xrefc(4,1,i))
1593 partsav(3,ip)=partsav(3,ip) + ems(i)*
1594 . (xrefc(1,2,i)+xrefc(2,2,i)+xrefc(3,2,i)+xrefc(4,2,i))
1595 partsav(4,ip)=partsav(4,ip) + ems(i)*
1596 . (xrefc(1,3,i)+xrefc(2,3,i)+xrefc(3,3,i)+xrefc(4,3,i))
1597 xx = xrefc(1,1,i)**2 + xrefc(2,1,i)**2
1598 . + xrefc(3,1,i)**2 + xrefc(4,1,i)
1599 xy = xrefc(1,1,i)*xrefc(1,2,i) + xrefc(2,1,i)*xrefc(2,2,i)
1600 . + xrefc(3,1,i)*xrefc(3,2,i) + xrefc(4,1,i)*xrefc(4,2,i)
1601 yy = xrefc(1,2,i)**2 + xrefc(2,2,i)**2
1602 . + xrefc(3,2,i)**2 + xrefc(4,2,i)
1603 yz = xrefc(1,2,i)*xrefc(1,3,i) + xrefc(2,2,i)*xrefc(2,3,i)
1604 . + xrefc(3,2,i)*xrefc(3,3,i) + xrefc(4,2,i)*xrefc(4,3,i)
1605 zz = xrefc(1,3,i)**2 + xrefc(2,3,i)**2
1606 . + xrefc(3,3,i)**2 + xrefc(4,3,i)
1607 zx = xrefc(1,3,i)*xrefc(1,1,i) + xrefc(2,3,i)*xrefc(2,1,i)
1608 . + xrefc(3,3,i)*xrefc(3,1,i) + xrefc(4,3,i)*xrefc(4,1,i)
1609C
1610 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1611 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1612 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1613 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1614 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1615 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1616C
1617 partsav(11,ip)=partsav(11,ip) + ems(i)*
1618 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1619 partsav(12,ip)=partsav(12,ip) + ems(i)*
1620 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1621 partsav(13,ip)=partsav(13,ip) + ems(i)*
1622 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1623C
1624 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1625 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1626 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1627 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1628 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1629 END IF
1630 END DO
1631 ENDIF
1632 END IF
1633
1634C----------------------------------------------
1635C initialization of nodal rigidities for interfaces
1636C----------------------------------------------
1637 IF (igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
1638 IF(i7stifs/=0)THEN
1639 DO i=lft,llt
1640 et= stack%PM(2,isubstack)*thk(i)
1641 stc(i)=et
1642 nshnod(ix1(i))=nshnod(ix1(i))+1
1643 nshnod(ix2(i))=nshnod(ix2(i))+1
1644 nshnod(ix3(i))=nshnod(ix3(i))+1
1645 nshnod(ix4(i))=nshnod(ix4(i))+1
1646 ENDDO
1647 ENDIF
1648 ELSEIF(igtyp == 11 .AND. igmat > 0) THEN
1649 ipgmat = 700
1650 IF(i7stifs/=0)THEN
1651 DO i=lft,llt
1652 et=geo(ipgmat +2 ,iprop)*thk(i)
1653 stc(i)=et
1654 nshnod(ix1(i))=nshnod(ix1(i))+1
1655 nshnod(ix2(i))=nshnod(ix2(i))+1
1656 nshnod(ix3(i))=nshnod(ix3(i))+1
1657 nshnod(ix4(i))=nshnod(ix4(i))+1
1658 ENDDO
1659 ENDIF
1660 ELSE ! igmat = 0
1661 IF(i7stifs/=0)THEN
1662 DO i=lft,llt
1663 et=pm(20,imid)*thk(i)
1664 stc(i)=et
1665 nshnod(ix1(i))=nshnod(ix1(i))+1
1666 nshnod(ix2(i))=nshnod(ix2(i))+1
1667 nshnod(ix3(i))=nshnod(ix3(i))+1
1668 nshnod(ix4(i))=nshnod(ix4(i))+1
1669 ENDDO
1670 ENDIF
1671 ENDIF ! igmat
1672 DEALLOCATE(posly,ratio_thkly,thk_it,pos_it )
1673C-----------
1674 RETURN
1675 END
cinmas
subroutine cinmas(x, xrefc, ix, geo, pm, ms, tiner, thke, ihbe, partsav, v, ipart, msc, inc, area, i8mi, igeo, etnod, imid, iprop, nshnod, stc, sh4tree, mcp, mcps, temp, ms_layer, zi_layer, ms_layerc, zi_layerc, msz2c, zply, isubstack, nlay, elbuf_str, stack, thki, rnoise, drape, nintemp, perturb, ix1, ix2, ix3, ix4, idrape, indx)
Definition cinmas.F:95
double_flot_ieee
subroutine double_flot_ieee(jft, jlt, i8, r8, i8f)
Definition cinmas.F:27
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21
drape_mod
Definition drape_mod.F:54
message_mod
Definition message_mod.F:1257
names_and_titles_mod
Definition names_and_titles_mod.F:1002
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1008
stack_mod
Definition stack_mod.F:75
ancmsg
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)
Definition message.F:895
fretitl2
subroutine fretitl2(titr, iasc, l)
Definition freform.F:799