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
103C----------------------------------------------
104C INITIALISATION DES MASSES ET INERTIES NODALES
105C----------------------------------------------
106C I m p l i c i t T y p e s
107C-----------------------------------------------
108#include "implicit_f.inc"
109C-----------------------------------------------
110C G l o b a l P a r a m e t e r s
111C-----------------------------------------------
112#include "mvsiz_p.inc"
113C-----------------------------------------------
114C C o m m o n B l o c k s
115C-----------------------------------------------
116#include "com01_c.inc"
117#include "com04_c.inc"
118#include "remesh_c.inc"
119#include "param_c.inc"
120#include "scr03_c.inc"
121#include "scr12_c.inc"
122#include "scr17_c.inc"
123#include "scr22_c.inc"
124#include "vect01_c.inc"
125#include "drape_c.inc"
126C-----------------------------------------------
127C D u m m y A r g u m e n t s
128C-----------------------------------------------
129 INTEGER IX(NIXC,*),IPART(*),IHBE,IMID,IPROP,
130 . IGEO(NPROPGI,*),NSHNOD(*), SH4TREE(KSH4TREE,*),
131 . ISUBSTACK,NLAY,PERTURB(NPERTURB),IDRAPE
132 INTEGER, DIMENSION(MVSIZ), INTENT(IN) :: IX1,IX2,IX3,IX4
133 INTEGER ,INTENT(IN) :: NINTEMP
134 INTEGER*8 I8MI(6,*)
135 INTEGER, DIMENSION(SCDRAPE) :: INDX
136 my_real, DIMENSION(MVSIZ), INTENT(IN) :: AREA
137 my_real
138 . X(3,*), GEO(NPROPG,*), PM(NPROPM,*), MS(*), TINER(*),THKE(*),
139 . v(3,*),partsav(20,*),msc(*),inc(*),
140 . etnod(*), stc(*),mcp(*),mcps(*),temp(*),
141 . ms_layer(numnod,*), zi_layer(numnod,*),xrefc(4,3,*),
142 . ms_layerc(numelc,*),zi_layerc(numelc,*),msz2c(*),zply(*),thki(*),
143 . rnoise(nperturb,*)
144 my_real
145 . zshift
146 TYPE(elbuf_struct_), TARGET :: ELBUF_STR
147 TYPE (STACK_PLY) :: STACK
148 TYPE (DRAPE_) , DIMENSION(NUMELC_DRAPE+NUMELTG_DRAPE), TARGET :: DRAPE
149C-----------------------------------------------
150C L o c a l V a r i a b l e s
151C-----------------------------------------------
152 INTEGER I,II,J,K,N,IP,I1,I2,I3,I4,MATLY,IPTHK,IPMAT,IPPOS,ID,
153 . ERRM,IDPROP,IPID,IPPID,ISTACK(MVSIZ,NPT),ISLV,IPID_LY,
154 . ipang,nthk,igtyp,igmat,ilay,nptt,it,iint,ipgmat,max_nptt,
155 . nslice,ipos,ipt,ipt_all,ie,nlay_max
156 CHARACTER(LEN=NCHARTITLE) :: TITR
157 my_real
158 . fac,mst,dms,thickt,thkly,laypos,mse,ins,mzi,ems_nptt,et,zj,
159 . xx,yy,zz,xy,yz,zx,msl,xil,ms_nptt,pos_0,thinning,f_offset,scale,
160 . thick_drape
161 my_real
162 . ems_layer(mvsiz,nlay)
163 my_real, DIMENSION(:), ALLOCATABLE :: posly,ratio_thkly, thk_it,pos_it
164 my_real, DIMENSION(MVSIZ) :: rho,rhocp,ems,emscp,xi,thk
165C
166
167 TYPE (DRAPE_PLY_), POINTER :: DRAPE_PLY
168 my_real A_GAUSS(9,9),W_GAUSS(9,9)
169C-----------------------------------------------
170 DATA A_GAUSS /
171 1 0. ,0. ,0. ,
172 1 0. ,0. ,0. ,
173 1 0. ,0. ,0. ,
174 2 -.577350269189626,0.577350269189626,0. ,
175 2 0. ,0. ,0. ,
176 2 0. ,0. ,0. ,
177 3 -.774596669241483,0. ,0.774596669241483,
178 3 0. ,0. ,0. ,
179 3 0. ,0. ,0. ,
180 4 -.861136311594053,-.339981043584856,0.339981043584856,
181 4 0.861136311594053,0. ,0. ,
182 4 0. ,0. ,0. ,
183 5 -.906179845938664,-.538469310105683,0. ,
184 5 0.538469310105683,0.906179845938664,0. ,
185 5 0. ,0. ,0. ,
186 6 -.932469514203152,-.661209386466265,-.238619186083197,
187 6 0.238619186083197,0.661209386466265,0.932469514203152,
188 6 0. ,0. ,0. ,
189 7 -.949107912342759,-.741531185599394,-.405845151377397,
190 7 0. ,0.405845151377397,0.741531185599394,
191 7 0.949107912342759,0. ,0. ,
192 8 -.960289856497536,-.796666477413627,-.525532409916329,
193 8 -.183434642495650,0.183434642495650,0.525532409916329,
194 8 0.796666477413627,0.960289856497536,0. ,
195 9 -.968160239507626,-.836031107326636,-.613371432700590,
196 9 -.324253423403809,0. ,0.324253423403809,
197 9 0.613371432700590,0.836031107326636,0.968160239507626/
198 DATA w_gauss /
199 1 2. ,0. ,0. ,
200 1 0. ,0. ,0. ,
201 1 0. ,0. ,0. ,
202 2 1. ,1. ,0. ,
203 2 0. ,0. ,0. ,
204 2 0. ,0. ,0. ,
205 3 0.555555555555556,0.888888888888889,0.555555555555556,
206 3 0. ,0. ,0. ,
207 3 0. ,0. ,0. ,
208 4 0.347854845137454,0.652145154862546,0.652145154862546,
209 4 0.347854845137454,0. ,0. ,
210 4 0. ,0. ,0. ,
211 5 0.236926885056189,0.478628670499366,0.568888888888889,
212 5 0.478628670499366,0.236926885056189,0. ,
213 5 0. ,0. ,0. ,
214 6 0.171324492379170,0.360761573048139,0.467913934572691,
215 6 0.467913934572691,0.360761573048139,0.171324492379170,
216 6 0. ,0. ,0. ,
217 7 0.129484966168870,0.279705391489277,0.381830050505119,
218 7 0.417959183673469,0.381830050505119,0.279705391489277,
219 7 0.129484966168870,0. ,0. ,
220 8 0.101228536290376,0.222381034453374,0.313706645877887,
221 8 0.362683783378362,0.362683783378362,0.313706645877887,
222 8 0.222381034453374,0.101228536290376,0. ,
223 9 0.081274388361574,0.180648160694857,0.260610696402935,
224 9 0.312347077040003,0.330239355001260,0.312347077040003,
225 9 0.260610696402935,0.180648160694857,0.081274388361574/
226C=======================================================================
227 igtyp=nint(geo(12,iprop))
228 igmat = igeo(98,iprop)
229 ipang = 1 ! NPT
230 ipthk = ipang + nlay ! NPT
231 ippos = ipthk + nlay ! NPT
232
233C
234 max_nptt = 1
235 IF(igtyp == 51 .OR. igtyp == 52) THEN
236 DO ilay=1,nlay ! NPT
237 max_nptt = max(max_nptt ,elbuf_str%BUFLY(ilay)%NPTT)
238 ENDDO
239 ALLOCATE(thk_it(max_nptt),pos_it(max_nptt))
240 ELSE
241 ALLOCATE(thk_it(0),pos_it(0))
242 ENDIF
243 nlay_max = max(nlay, npt)
244 ALLOCATE(posly(nlay_max*max_nptt),ratio_thkly(nlay_max*max_nptt))
245c-----------------------------
246 IF ((igtyp == 17 .OR. igtyp == 51) .AND. igmat < 0 ) THEN
247c-----------------------------
248 DO i=lft,llt
249 IF (ndrape == 0) THEN
250 IF (thke(i)== zero) THEN
251 thk(i) = stack%GEO(1,isubstack)
252 thke(i)= thk(i)
253 ELSE
254 thk(i)=thke(i)
255 ENDIF
256 thki(i) = stack%GEO(1,isubstack)
257 ELSEIF (ndrape > 0) THEN
258 thk(i) = thke(i)
259 thki(i)= thk(i)
260 ENDIF ! IF (NDRAPE == 0)
261c
262 IF (nperturb /= 0) THEN
263 DO j=1,nperturb
264 IF (perturb(j) == 1 .AND. rnoise(j,i+nft) /= zero) THEN
265 thk(i) = thk(i) * rnoise(j,i+nft)
266 thke(i) = thke(i) * rnoise(j,i+nft)
267 thki(i) = thki(i) * rnoise(j,i+nft)
268 ENDIF
269 ENDDO
270 ENDIF
271c
272 rho(i) = pm(89,imid)
273 rhocp(i) = pm(69,imid)
274 IF(thk(i) == zero.AND.igtyp/=0)THEN
275 CALL ancmsg(msgid=495,
276 . msgtype=msgerror,
277 . anmode=aninfo_blind_1,
278 . i1=ix(7,nft+i))
279 ENDIF
280 ENDDO
281c-----------------------------
282 ELSEIF (igtyp == 52 .OR.
283 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))THEN
284c-----------------------------
285 DO i=lft,llt
286 IF (ndrape == 0) THEN
287 IF (thke(i)== zero) THEN
288 thk(i) = stack%GEO(1,isubstack)
289 thke(i)= thk(i)
290 ELSE
291 thk(i)=thke(i)
292 ENDIF
293 thki(i) = stack%GEO(1,isubstack)
294 ELSEIF (ndrape > 0) THEN
295 thk(i) = thke(i)
296 thki(i)= thk(i)
297 ENDIF ! IF (NDRAPE == 0)
298 rho(i) = stack%PM(1,isubstack)
299 rhocp(i) = stack%PM(15,isubstack)* rho(i)/stack%PM(14,isubstack)
300c
301 IF (nperturb /= 0) THEN
302 DO j=1,nperturb
303 IF (perturb(j) == 1 .AND. rnoise(j,i+nft) /= zero) THEN
304 thk(i) = thk(i) * rnoise(j,i+nft)
305 thke(i) = thke(i) * rnoise(j,i+nft)
306 thki(i) = thki(i) * rnoise(j,i+nft)
307 ENDIF
308 ENDDO
309 ENDIF
310c
311 IF(thk(i) == zero.AND.igtyp/=0)THEN
312 CALL ancmsg(msgid=495,
313 . msgtype=msgerror,
314 . anmode=aninfo_blind_1,
315 . i1=ix(7,nft+i))
316 ENDIF
317 ENDDO
318c-----------------------------
319 ELSE ! IGTYP : std shell property
320c-----------------------------
321 DO i=lft,llt
322
323 IF (thke(i) == zero) THEN
324 thk(i) = geo(1,iprop)
325 thke(i)= thk(i)
326 ELSE
327 thk(i)=thke(i)
328 ENDIF
329c
330 thki(i) = geo(1,iprop)
331c
332 IF (nperturb /= 0) THEN
333 DO j=1,nperturb
334 IF (perturb(j) == 1 .AND. rnoise(j,i+nft) /= zero) THEN
335 thk(i) = thk(i) * rnoise(j,i+nft)
336 thke(i) = thke(i) * rnoise(j,i+nft)
337 thki(i) = thki(i) * rnoise(j,i+nft)
338 ENDIF
339 ENDDO
340 ENDIF
341c
342 rho(i) = pm(89,imid)
343 rhocp(i) = pm(69,imid)
344 IF(thk(i) == zero.AND.igtyp/=0)THEN
345 CALL ancmsg(msgid=495,
346 . msgtype=msgerror,
347 . anmode=aninfo_blind_1,
348 . i1=ix(7,nft+i))
349 ENDIF
350 ENDDO
351 ENDIF
352C----------------------------------------------
353C MASSES ELEMENTS / 4
354C----------------------------------------------
355 IF (igtyp == 11 .OR. igtyp == 16) THEN
356c--------------------------------------------
357 ipthk = 300
358 ipmat = 100
359 DO i=lft,llt
360 ems(i)=zero
361 DO n=1,npt
362 i1=ipthk+n
363 i2=ipmat+n
364 thkly = geo(i1,iprop)*thk(i)
365 matly = igeo(i2,iprop)
366 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth
367 ENDDO
368 mst = rho(i)*thk(i)*area(i)*fourth
369 dms = abs(ems(i)-mst)/mst
370C! for igty=11 ---> IGMAT=-1 (old G-mat), 1 for new g-mat
371C! for igtyp=16 --->IGMAT= 0
372C! for igtyp=17 --->IGMAT= 0
373 IF(igmat <= 0) THEN
374 IF (dms > em02) errm = 1
375 IF (dms > em02) THEN
376 idprop = igeo(1,iprop)
377 CALL fretitl2(titr,
378 . igeo(npropgi-ltitr+1,iprop),ltitr)
379 CALL ancmsg(msgid=519,
380 . msgtype=msgwarning,
381 . anmode=aninfo_blind_2,
382 . i1=idprop,
383 . c1=titr,i2=ix(nixc,nft+i),
384 . r1=ems(i),r2=mst)
385 ENDIF
386 ENDIF
387 ENDDO
388c--------------------------------------------
389 ELSEIF(igtyp == 17 ) THEN
390c--------------------------------------------
391 ippid = 2
392 ipmat = ippid + npt
393 ipang = 1
394 ipthk = ipang + npt
395 ippos = ipthk + npt
396 nthk = ippos + npt
397C
398 IF(idrape == 0) THEN
399 DO i=lft,llt
400 ems(i)=zero
401 mzi = zero
402 DO n=1,npt
403 i1 = ippid + n
404 i2 = ipmat + n
405 i3 = ipthk + n
406 i4 = ippos + n
407 thickt = stack%GEO(1 ,isubstack)
408 thkly = stack%GEO(i3 ,isubstack)*thk(i)
409 matly = stack%IGEO(i2,isubstack)
410 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth
411 ems_layer(i,n) = pm(89,matly)*thkly*area(i)*fourth
412 mzi = mzi + ems_layer(i,n)*stack%GEO(i4 ,isubstack)
413 ipid = stack%IGEO(i1, isubstack)
414 istack(i,n) = igeo(102,ipid)
415 ENDDO
416 mst = rho(i)*thk(i)*area(i)*fourth
417 dms = abs(ems(i)-mst)/mst
418 IF(dms > em02) errm = 1
419 IF(igmat <= 0) THEN
420 IF(dms > em02) THEN
421 idprop = igeo(1,iprop)
422 CALL fretitl2(titr,
423 . igeo(npropgi-ltitr+1,iprop),ltitr)
424 CALL ancmsg(msgid=519,
425 . msgtype=msgwarning,
426 . anmode=aninfo_blind_2,
427 . i1=idprop,
428 . c1=titr,i2=ix(nixc,nft+i),
429 . r1=ems(i),r2=mst)
430 ENDIF
431 ENDIF
432 !!
433 DO n=1,npt
434 i4 = ippos + n
435 IF(abs(stack%GEO(i4 ,isubstack)) < em15)
436 . stack%GEO(i4 ,isubstack)=zero
437 ENDDO
438 IF(abs(mzi) > em02) THEN
439 DO n=1,npt
440 i4 = ippos + n
441 dms = mzi/ems(i)
442 stack%GEO(i4 ,isubstack) = stack%GEO(i4 ,isubstack) - dms
443 ENDDO
444 ENDIF
445 ENDDO
446 ELSE ! idrape > 0
447 DO i=lft,llt
448 ems(i)=zero
449 mzi = zero
450 ie = indx(nft+i)
451 IF(ie == 0) THEN
452 DO n=1,npt
453 i1 = ippid + n
454 i2 = ipmat + n
455 i3 = ipthk + n
456 i4 = ippos + n
457 thickt = stack%GEO(1 ,isubstack)
458 thkly = stack%GEO(i3 ,isubstack)*thk(i)
459 matly = stack%IGEO(i2,isubstack)
460 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth
461 ems_layer(i,n) = pm(89,matly)*thkly*area(i)*fourth
462 mzi = mzi + ems_layer(i,n)*stack%GEO(i4 ,isubstack)
463 ipid = stack%IGEO(i1, isubstack)
464 istack(i,n) = igeo(102,ipid)
465 ENDDO
466 ELSE
467 thick_drape = drape(ie)%THICK
468 scale = thk(i)/thick_drape
469 DO n=1,npt
470 i1 = ippid + n
471 i2 = ipmat + n
472 i3 = ipthk + n
473 i4 = ippos + n
474 ip = drape(ie)%INDX_PLY(n)
475 IF(ip > 0 ) THEN
476 drape_ply => drape(ie)%DRAPE_PLY(ip)
477 thinning = drape_ply%RDRAPE(1,1) ! one slice by layer
478 thickt = stack%GEO(1 ,isubstack)
479 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
480 thkly = thkly*thinning
481 !! SCALE = THK(I)/THICK_DRAPE
482 matly = stack%IGEO(i2,isubstack)
483 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth*scale
484 ems_layer(i,n) = pm(89,matly)*thkly*area(i)*fourth *scale
485 ratio_thkly(n) = thkly/thick_drape
486 IF (n == 1) THEN
487 posly(n) = -half + half*ratio_thkly(n)
488 ELSE
489 posly(n) = posly(n-1)
490 . + half*(ratio_thkly(n)+ratio_thkly(n-1))
491 ENDIF ! IF (N == 1)
492
493 mzi = mzi + ems_layer(i,n)*posly(n)
494 ipid = stack%IGEO(i1, isubstack)
495 istack(i,n) = igeo(102,ipid)
496 ELSE ! IP=0
497 thickt = stack%GEO(1 ,isubstack)
498 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
499 matly = stack%IGEO(i2,isubstack)
500 !! SCALE = THK(I)/THICK_DRAPE
501 ems(i) = ems(i) + pm(89,matly)*thkly*area(i)*fourth*scale
502 ems_layer(i,n) = pm(89,matly)*thkly*area(i)*fourth*scale
503 ratio_thkly(n) = thkly/thick_drape
504 IF (n == 1) THEN
505 posly(n) = -half + half*ratio_thkly(n)
506 ELSE
507 posly(n) = posly(n-1)
508 . + half*(ratio_thkly(n)+ratio_thkly(n-1))
509 ENDIF ! IF (N == 1)
510 mzi = mzi + ems_layer(i,n)*posly(n)
511 ipid = stack%IGEO(i1, isubstack)
512 istack(i,n) = igeo(102,ipid)
513 ENDIF
514 ENDDO
515 ENDIF ! IE
516 mst = rho(i)*thk(i)*area(i)*fourth
517 dms = abs(ems(i)-mst)/mst
518 IF(dms > em02) errm = 1
519 IF(igmat <= 0) THEN
520 IF(dms > em02) THEN
521 idprop = igeo(1,iprop)
522 CALL fretitl2(titr,
523 . igeo(npropgi-ltitr+1,iprop),ltitr)
524 CALL ancmsg(msgid=519,
525 . msgtype=msgwarning,
526 . anmode=aninfo_blind_2,
527 . i1=idprop,
528 . c1=titr,i2=ix(nixc,nft+i),
529 . r1=ems(i),r2=mst)
530 ENDIF
531 ENDIF
532C correction la position de la surface neutre au CDG
533 DO n=1,npt
534 i4 = ippos + n
535 IF(abs(stack%GEO(i4 ,isubstack)) < em15)
536 . stack%GEO(i4 ,isubstack)=zero
537 ENDDO
538 IF(abs(mzi) > em02) THEN
539 DO n=1,npt
540 i4 = ippos + n
541 dms = mzi/ems(i)
542 stack%GEO(i4 ,isubstack) = stack%GEO(i4 ,isubstack) - dms
543 ENDDO
544 ENDIF
545 ENDDO
546 ENDIF ! IDRAPE
547c--------------------------------------------
548 ELSEIF (igtyp == 51 .OR. igtyp == 52 ) THEN
549c--------------------------------------------
550
551 ipang = 1
552 ippid = 2
553 ipmat = ippid + nlay ! NPT
554 ipthk = ipang + nlay ! NPT
555 ippos = ipthk + nlay ! NPT
556 nthk = ippos + nlay ! NPT
557C
558 ipos = igeo(99,iprop)
559 thickt = stack%GEO(1,isubstack)
560 zshift = geo(199,iprop)
561 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
562 IF(idrape == 0 ) THEN
563 DO i=lft,llt
564 ems(i)=zero
565 mzi = zero
566 ipt_all = 0
567 scale = thk(i)/thickt
568 DO ilay=1,nlay ! NPT
569 nptt = elbuf_str%BUFLY(ilay)%NPTT
570 i1 = ippid + ilay
571 i2 = ipmat + ilay
572 i3 = ipthk + ilay
573 i4 = ippos + ilay
574 ipid = stack%IGEO(i1,isubstack) ! or ipid = ix(nixc-1,nft+i)
575 iint = igeo(47,iprop)
576 thickt = stack%GEO(1,isubstack)
577 IF(iint == 1) THEN
578 thkly = stack%GEO(i3,isubstack)*thickt
579 laypos = stack%GEO(i4,isubstack)
580 ipid_ly = stack%IGEO(i1,isubstack)
581 matly = stack%IGEO(i2,isubstack)
582 DO it=1,nptt
583 ipt = ipt_all + it
584 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
585 ratio_thkly(ipt) = thk_it(it)/thickt
586 IF (ipt == 1) THEN
587 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
588 ELSE
589 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
590 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
591 ENDIF ! IF (IPT == 1)
592 pos_it(it) = posly(ipt)*thk(i)
593 ENDDO
594 ELSEIF(iint == 2) THEN
595 thkly = stack%GEO(i3,isubstack)*thickt
596 laypos = stack%GEO(i4,isubstack)
597 ipid_ly = stack%IGEO(i1,isubstack)
598 matly = stack%IGEO(i2,isubstack)
599 DO it=1,nptt
600 ipt = ipt_all + it
601 thk_it(it) = half*thkly*w_gauss(it,nptt) ! Gauss distribution of NPTT through layer
602 ratio_thkly(ipt) = thk_it(it)/thickt
603 IF (ipt == 1 ) THEN
604 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
605 ELSE
606 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
607 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
608 ENDIF ! IF (IPT == 1)
609 pos_it(it) = posly(ipt)*thk(i)
610 ENDDO
611 ENDIF ! iint
612 !!SCALE = THK(I)/THICKT
613 DO it=1,nptt
614 thk_it(it) = thk_it(it)*scale
615 ems_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
616 ems(i) = ems(i) + ems_nptt
617 mzi = mzi + ems_nptt * pos_it(it)/thk(i)
618 ENDDO
619 ipt_all = ipt_all + nptt
620 ENDDO ! NLAY
621 !!
622 mst = rho(i)*thk(i)*area(i)*fourth
623 dms = abs(ems(i)-mst)/mst
624 IF (dms > em02) errm = 1
625 IF(igtyp == 51 .AND. igmat < 0) THEN
626 IF (dms > em02) THEN
627 idprop = igeo(1,iprop)
628 CALL fretitl2(titr,
629 . igeo(npropgi-ltitr+1,iprop),ltitr)
630 CALL ancmsg(msgid=519,
631 . msgtype=msgwarning,
632 . anmode=aninfo_blind_2,
633 . i1=idprop,
634 . c1=titr,i2=ix(nixc,nft+i),
635 . r1=ems(i),r2=mst)
636 ENDIF
637 ENDIF
638 ENDDO ! lft,llt
639C------
640 IF( jthe > 0 ) THEN
641 DO i=lft,llt
642 emscp(i)=rhocp(i)*thk(i)*area(i)*fourth
643 ENDDO
644 ENDIF
645 ELSE ! idrape > 0
646 DO i=lft,llt
647 ems(i)=zero
648 mzi = zero
649 ipt_all = 0
650 ie = indx(nft+i)
651 IF(ie == 0) THEN
652 scale = thk(i)/thickt
653 DO ilay=1,nlay ! NPT
654 nptt = elbuf_str%BUFLY(ilay)%NPTT
655 i1 = ippid + ilay
656 i2 = ipmat + ilay
657 i3 = ipthk + ilay
658 i4 = ippos + ilay
659 ipid = stack%IGEO(i1,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
660 iint = igeo(47,iprop)
661 thickt = stack%GEO(1,isubstack)
662 IF(iint == 1) THEN
663 thkly = stack%GEO(i3,isubstack)*thickt
664 laypos = stack%GEO(i4,isubstack)
665 ipid_ly = stack%IGEO(i1,isubstack)
666 matly = stack%IGEO(i2,isubstack)
667 DO it=1,nptt
668 ipt = ipt_all + it
669 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
670 ratio_thkly(ipt) = thk_it(it)/thickt
671 IF (ipt == 1) THEN
672 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
673 ELSE
674 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
675 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
676 ENDIF ! IF (IPT == 1)
677 pos_it(it) = posly(ipt)*thk(i)
678 ENDDO
679 ELSEIF(iint == 2) THEN
680 thkly = stack%GEO(i3,isubstack)*thickt
681 laypos = stack%GEO(i4,isubstack)
682 ipid_ly = stack%IGEO(i1,isubstack)
683 matly = stack%IGEO(i2,isubstack)
684 DO it=1,nptt
685 ipt = ipt_all + it
686 thk_it(it) = half*thkly*w_gauss(it,nptt) ! Gauss distribution of NPTT through layer
687 ratio_thkly(ipt) = thk_it(it)/thickt
688 IF (ipt == 1 ) THEN
689 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
690 ELSE
691 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
692 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
693 ENDIF ! IF (IPT == 1)
694 pos_it(it) = posly(ipt)*thk(i)
695 ENDDO
696 ENDIF ! iint
697 !! SCALE = THK(I)/THICK
698 DO it=1,nptt
699 thk_it(it) = scale*thk_it(it)
700 ems_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
701 ems(i) = ems(i) + ems_nptt
702 mzi = mzi + ems_nptt * pos_it(it)/thk(i)
703 ENDDO
704 ipt_all = ipt_all + nptt
705 ENDDO ! NLAY
706 ELSE ! IE
707 thick_drape = drape(ie)%THICK
708 scale = thk(i)/thick_drape
709 DO ilay=1,nlay ! NPT
710 ip = drape(ie)%INDX_PLY(ilay)
711 nptt = elbuf_str%BUFLY(ilay)%NPTT ! = NSLICE
712 IF(ip > 0 ) THEN
713 drape_ply => drape(ie)%DRAPE_PLY(ip)
714 nslice = drape_ply%NSLICE
715 i1 = ippid + ilay
716 i2 = ipmat + ilay
717 i3 = ipthk + ilay
718 i4 = ippos + ilay
719 ipid = stack%IGEO(i1,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
720 iint = igeo(47,iprop)
721 thickt = stack%GEO(1,isubstack)
722 IF(iint == 1) THEN
723 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
724 ipid_ly = stack%IGEO(i1,isubstack)
725 matly = stack%IGEO(i2,isubstack)
726 DO it=1,nptt
727 ipt = ipt_all + it
728 thinning = drape_ply%RDRAPE(it,1)
729 thk_it(it) = thkly*thinning/nptt
730 ratio_thkly(ipt) = thk_it(it)/thick_drape
731 IF (ipt == 1) THEN
732 posly(ipt) = zshift + half*ratio_thkly(ipt)
733 ELSE
734 posly(ipt) = posly(ipt-1)
735 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
736 ENDIF ! IF (IPT == 1)
737 pos_it(it) = posly(ipt)*thk(i)
738 ENDDO ! nptt
739 ELSEIF(iint == 2) THEN
740 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
741 ipid_ly = stack%IGEO(i1,isubstack)
742 matly = stack%IGEO(i2,isubstack)
743 DO it=1,nptt
744 ipt = ipt_all + it
745 thinning = drape_ply%RDRAPE(it,1)
746 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
747 ratio_thkly(ipt) = thk_it(it)/thick_drape
748 IF (ipt == 1) THEN
749 posly(ipt) = zshift + half*ratio_thkly(ipt)
750 ELSE
751 posly(ipt) = posly(ipt-1)
752 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
753 ENDIF ! IF (IPT == 1)
754 pos_it(it) = posly(ipt)*thk(i)
755 ENDDO ! nptt
756 ENDIF ! iint
757 !! SCALE = THK(I)/THICK_DRAPE
758 DO it=1,nptt
759 thk_it(it)= scale*thk_it(it)
760 ems_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
761 ems(i) = ems(i) + ems_nptt
762 mzi = mzi + ems_nptt * pos_it(it)/thk(i)
763 ENDDO
764 ELSE ! Ip not draped
765 nptt = elbuf_str%BUFLY(ilay)%NPTT
766 i1 = ippid + ilay
767 i2 = ipmat + ilay
768 i3 = ipthk + ilay
769 i4 = ippos + ilay
770 ipid = stack%IGEO(i1,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
771 iint = igeo(47,iprop)
772 thickt = stack%GEO(1,isubstack)
773 IF(iint == 1) THEN
774 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
775 ipid_ly = stack%IGEO(i1,isubstack)
776 matly = stack%IGEO(i2,isubstack)
777 DO it=1,nptt
778 ipt = ipt_all + it
779 thk_it(it) = thkly/nptt
780 ratio_thkly(ipt) = thk_it(it)/thick_drape
781 IF (ipt == 1) THEN
782 posly(ipt) = zshift + half*ratio_thkly(ipt)
783 ELSE
784 posly(ipt) = posly(ipt-1)
785 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
786 ENDIF ! IF (IPT == 1)
787 pos_it(it) = posly(ipt)*thk(i)
788 ENDDO ! nptt
789 ELSEIF(iint == 2) THEN
790 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
791 ipid_ly = stack%IGEO(i1,isubstack)
792 matly = stack%IGEO(i2,isubstack)
793 DO it=1,nptt
794 ipt = ipt_all + it
795 thk_it(it) = half*thkly*w_gauss(it,nptt)
796 ratio_thkly(ipt) = thk_it(it)/thk(i)/thick_drape
797 IF (ipt == 1) THEN
798 posly(ipt) = zshift + half*ratio_thkly(ipt)
799 ELSE
800 posly(ipt) = posly(ipt-1)
801 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
802 ENDIF ! IF (IPT == 1)
803 pos_it(it) = posly(ipt)*thk(i)
804 ENDDO ! nptt
805 ENDIF ! iint
806 !! SCALE = THK(I)/THICK_DRAPE
807 DO it=1,nptt
808 thk_it(it ) = scale*thk_it(it)
809 ems_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
810 ems(i) = ems(i) + ems_nptt
811 mzi = mzi + ems_nptt * pos_it(it)/thk(i)
812 ENDDO
813 ENDIF ! IP
814 ipt_all = ipt_all + nptt
815 ENDDO ! NLAY
816 ENDIF ! IE
817
818 mst = rho(i)*thk(i)*area(i)*fourth
819 dms = abs(ems(i)-mst)/mst
820 IF (dms > em02) errm = 1
821 IF(igtyp == 51 .AND. igmat < 0) THEN
822 IF (dms > em02) THEN
823 idprop = igeo(1,iprop)
824 CALL fretitl2(titr,
825 . igeo(npropgi-ltitr+1,iprop),ltitr)
826 CALL ancmsg(msgid=519,
827 . msgtype=msgwarning,
828 . anmode=aninfo_blind_2,
829 . i1=idprop,
830 . c1=titr,i2=ix(nixc,nft+i),
831 . r1=ems(i),r2=mst)
832 ENDIF
833 ENDIF
834C
835 ENDDO
836C------
837 IF( jthe > 0 ) THEN
838 DO i=lft,llt
839 emscp(i)=rhocp(i)*thk(i)*area(i)*fourth
840 ENDDO
841 ENDIF
842 ENDIF ! idrape
843c--------------------------------------------
844 ELSE ! IGTYP
845c--------------------------------------------
846 DO i=lft,llt
847 ems(i)=rho(i)*thk(i)*area(i)*fourth
848 ENDDO
849C
850 IF( jthe > 0 ) THEN
851 DO i=lft,llt
852 emscp(i) = rhocp(i)*thk(i)*area(i)*fourth
853 ENDDO
854 ENDIF
855C
856 ENDIF ! IGTYP
857C----------------------------------------------
858C INITIALISATION DES MASSES NODALES
859c elt mass (parith/on spmd + adaptive meshing)
860C----------------------------------------------
861 IF(jthe == 0 ) THEN
862 DO i=lft,llt
863 msc(i)=ems(i)
864 ENDDO
865 ELSE
866 DO i=lft,llt
867 msc(i)=ems(i)
868 mcps(i) = emscp(i)
869 ENDDO
870 ENDIF
871c
872 IF(ishxfem_ply > 0 )THEN
873 ippid = 2
874 ipmat = ippid + npt
875 ipang = 1
876 ipthk = ipang + npt
877 ippos = ipthk + npt
878 nthk = ippos + npt
879 DO i=lft,llt
880 ii = nft + i
881 DO j=1,npt
882 ip = istack(i,j)
883 ms_layerc(ii,ip)= ems_layer(i,j)
884C
885 i4 = ippos + j
886 thickt = stack%GEO(1 ,isubstack)
887 islv = igeo(32,iprop)
888 IF(islv == 0)THEN
889 zi_layerc(ii,ip) = thickt*stack%GEO(i4 ,isubstack)
890 zply(ip) = zi_layerc(ii,ip)
891 ENDIF
892cc ZJ = ZI_LAYERC(II,IP)
893 zj = thickt*stack%GEO(i4 ,isubstack)
894 msz2c(ii) = msz2c(ii) + ems_layer(i,j)*zj*zj
895 ENDDO
896 ENDDO
897 ENDIF
898C
899 IF (jthe > 0 .or. nintemp > 0) THEN
900 IF (nintemp > 0 ) THEN
901 DO i= lft,llt
902 IF(temp(ix1(i))== zero) temp(ix1(i)) = pm(79,imid)
903 IF(temp(ix2(i))== zero) temp(ix2(i)) = pm(79,imid)
904 IF(temp(ix3(i))== zero) temp(ix3(i)) = pm(79,imid)
905 IF(temp(ix4(i))== zero) temp(ix4(i)) = pm(79,imid)
906 ENDDO
907 ELSE
908 DO i=lft,llt
909 temp(ix1(i))= pm(79,imid)
910 temp(ix2(i))= pm(79,imid)
911 temp(ix3(i))= pm(79,imid)
912 temp(ix4(i))= pm(79,imid)
913 ENDDO
914 ENDIF
915 ENDIF
916C----------------------------------------------
917C INERTIES ELEMENTS /4
918C----------------------------------------------
919 IF(iner_9_12 /= zero)THEN
920 fac=iner_9_12
921 ELSEIF(ihbe>=11)THEN
922 fac=twelve
923 ELSE
924 fac=nine
925 ENDIF
926 IF (igtyp == 11 .OR. igtyp == 16) THEN
927 ipthk = 300
928 ippos = 400
929 DO i=lft,llt
930 xi(i)=zero
931 DO n=1,npt
932 i1=ipthk+n
933 i3=ippos+n
934 thickt= geo(200,iprop)
935 thkly = geo(i1,iprop)*thk(i)
936 laypos = geo(i3,iprop)*thk(i)
937 i2=ipmat+n
938 matly = igeo(i2,iprop)
939 msl = pm(89,matly)*thkly*area(i)*fourth
940 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
941 xi(i) = xi(i) + xil
942 ENDDO
943 ENDDO
944 ELSEIF(igtyp == 17) THEN
945 ippid = 2
946 ipmat = ippid + npt
947 ipang = 1
948 ipthk = ipang + npt
949 ippos = ipthk + npt
950 nthk = ippos + npt
951 ipos = igeo(99,iprop)
952 thickt = stack%GEO(1,isubstack)
953 zshift = geo(199,iprop)
954 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
955C
956 IF(idrape == 0) THEN
957 DO i=lft,llt
958 xi(i)=zero
959 DO n=1,npt
960 i1 = ippid + n
961 i2 = ipmat + n
962 i3 = ipthk + n
963 i4 = ippos + n
964 thickt= stack%GEO(1,isubstack)
965 thkly = stack%GEO(i3 ,isubstack )*thk(i)
966 laypos = stack%GEO(i4,isubstack)*thk(i)
967 matly = stack%IGEO(i2 ,isubstack)
968 msl = pm(89,matly)*thkly*area(i)*fourth
969 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
970 xi(i) = xi(i) + xil
971 ENDDO
972 ENDDO
973 ELSE ! IDRAPE > 0
974 DO i=lft,llt
975 xi(i)=zero
976 ie = indx(nft + i)
977 IF(ie == 0) THEN
978 DO n=1,npt
979 i1 = ippid + n
980 i2 = ipmat + n
981 i3 = ipthk + n
982 i4 = ippos + n
983 thickt= stack%GEO(1,isubstack)
984 thkly = stack%GEO(i3 ,isubstack )*thk(i)
985 laypos = stack%GEO(i4,isubstack)*thk(i)
986 matly = stack%IGEO(i2 ,isubstack)
987 msl = pm(89,matly)*thkly*area(i)*fourth
988 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
989 xi(i) = xi(i) + xil
990 ENDDO
991 ELSE
992 thick_drape = drape(ie)%THICK
993 scale = thk(i)/thick_drape
994 DO n=1,npt
995 ip = drape(ie)%INDX_PLY(n)
996 IF(ip > 0) THEN
997 drape_ply => drape(ie)%DRAPE_PLY(ip)
998 i1 = ippid + n
999 i2 = ipmat + n
1000 i3 = ipthk + n
1001 i4 = ippos + n
1002 thinning = drape_ply%RDRAPE(1,1)
1003 thickt= stack%GEO(1,isubstack)
1004 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1005 thkly = thkly*thinning
1006 matly = stack%IGEO(i2,isubstack)
1007 ratio_thkly(n) = thkly/thick_drape
1008 IF (n == 1) THEN
1009 posly(n) = zshift + half*ratio_thkly(n)
1010 ELSE
1011 posly(n) = posly(n-1)
1012 . + half*(ratio_thkly(n)+ratio_thkly(n-1))
1013 ENDIF ! IF (N == 1)
1014 laypos = posly(n)*thk(i)
1015 thkly = thkly*scale
1016 msl = pm(89,matly)*thkly*area(i)*fourth
1017 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
1018 xi(i) = xi(i) + xil
1019 ELSE !IP ==0
1020 i1 = ippid + n
1021 i2 = ipmat + n
1022 i3 = ipthk + n
1023 i4 = ippos + n
1024 thickt= stack%GEO(1,isubstack)
1025 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1026 matly = stack%IGEO(i2,isubstack)
1027 ratio_thkly(n) = thkly/thick_drape
1028 IF (n == 1) THEN
1029 posly(n) = zshift + half*ratio_thkly(n)
1030 ELSE
1031 posly(n) = posly(n-1)
1032 . + half*(ratio_thkly(n)+ratio_thkly(n-1))
1033 ENDIF ! IF (N == 1)
1034 laypos = posly(n)*thk(i)
1035 thkly = thkly*scale
1036 msl = pm(89,matly)*thkly*area(i)*fourth
1037 xil = msl*(area(i)/fac+thkly*thkly*one_over_12+laypos*laypos)
1038 xi(i) = xi(i) + xil
1039 ENDIF
1040 ENDDO
1041 ENDIF
1042 ENDDO
1043 ENDIF
1044 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
1045 ipang = 1
1046 ippid = 2
1047 ipmat = ippid + nlay ! NPT
1048 ipthk = ipang + nlay ! NPT
1049 ippos = ipthk + nlay ! NPT
1050 nthk = ippos + nlay ! NPT
1051 ipos = igeo(99,iprop)
1052 thickt = stack%GEO(1,isubstack)
1053 zshift = geo(199,iprop)
1054 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
1055C
1056 IF(idrape == 0) THEN
1057 DO i=lft,llt
1058 scale = thk(i)/thickt
1059 xi(i)=zero
1060 ipt_all = 0
1061 DO ilay=1, nlay ! NPT
1062 nptt = elbuf_str%BUFLY(ilay)%NPTT
1063 i1 = ippid + ilay
1064 i2 = ipmat + ilay
1065 i3 = ipthk + ilay
1066 i4 = ippos + ilay
1067 ipid_ly = stack%IGEO(i1,isubstack)
1068 matly = stack%IGEO(i2,isubstack)
1069 iint = igeo(47,iprop)
1070 IF(iint == 1) THEN
1071 thickt = stack%GEO(1,isubstack)
1072 thkly = stack%GEO(i3,isubstack)*thickt
1073 laypos = stack%GEO(i4 ,isubstack)*thickt
1074 DO it=1,nptt
1075 ipt = ipt_all + it
1076 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
1077 ratio_thkly(ipt) = thk_it(it)/thickt
1078 IF (ipt == 1 ) THEN
1079 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
1080 ELSE
1081 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
1082 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
1083 ENDIF ! IF (IPT == 1)
1084 pos_it(it) = posly(ipt)*thk(i)
1085 ENDDO
1086 ELSEIF(iint == 2) THEN
1087 thickt = stack%GEO(1,isubstack)
1088 thkly = stack%GEO(i3,isubstack)*thickt
1089 laypos = stack%GEO(i4 ,isubstack)*thickt
1090 DO it=1,nptt
1091 ipt = ipt_all + it
1092 thk_it(it) = half*thkly*w_gauss(it,nptt)
1093 ratio_thkly(ipt) = thk_it(it)/thickt
1094 IF (ipt == 1) THEN
1095 posly(ipt) = zshift + half*ratio_thkly(ipt)
1096 ELSE
1097 posly(ipt) = posly(ipt-1)
1098 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1099 ENDIF ! IF (IPT 1)
1100 pos_it(it) = posly(ipt)*thk(i)
1101 ENDDO ! nptt
1102 ENDIF ! iint
1103C
1104 !!SCALE = THK(I)/THICKT
1105 DO it=1,nptt
1106 thk_it(it) = scale*thk_it(it)
1107 ms_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
1108 xil = ms_nptt*(area(i)/fac+thk_it(it)*thk_it(it)*one_over_12+
1109 . pos_it(it)*pos_it(it))
1110 xi(i) = xi(i) + xil
1111 ENDDO
1112 ipt_all = ipt_all + nptt
1113 ENDDO ! NLAY
1114 ENDDO ! LFT,LLT
1115C
1116 ELSE ! IDRAPE > 0
1117 DO i=lft,llt
1118 xi(i)=zero
1119 ipt_all = 0
1120 ie = indx(nft + i)
1121 IF(ie == 0) THEN
1122 scale = thk(i)/thickt
1123 DO ilay=1, nlay ! NPT
1124 nptt = elbuf_str%BUFLY(ilay)%NPTT
1125 i1 = ippid + ilay
1126 i2 = ipmat + ilay
1127 i3 = ipthk + ilay
1128 i4 = ippos + ilay
1129 ipid_ly = stack%IGEO(i1,isubstack)
1130 matly = stack%IGEO(i2,isubstack)
1131 iint = igeo(47,iprop)
1132 IF(iint == 1) THEN
1133 thickt = stack%GEO(1,isubstack)
1134 thkly = stack%GEO(i3,isubstack)*thickt
1135 laypos = stack%GEO(i4 ,isubstack)*thickt
1136 DO it=1,nptt
1137 ipt = ipt_all + it
1138 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
1139 ratio_thkly(ipt) = thk_it(it)/thickt
1140 IF (ipt == 1 ) THEN
1141 posly(ipt) = zshift + half*ratio_thkly(ipt) ! integr. point "IT" position ratio
1142 ELSE
1143 posly(ipt) = posly(ipt-1) ! integr. point "IT" position ratio
1144 . + half*(ratio_thkly(ipt) + ratio_thkly(ipt-1))
1145 ENDIF ! IF (IPT == 1)
1146 pos_it(it) = posly(ipt)*thk(i)
1147 ENDDO
1148 ELSEIF(iint == 2) THEN
1149 thickt = stack%GEO(1,isubstack)
1150 thkly = stack%GEO(i3,isubstack)*thickt
1151 laypos = stack%GEO(i4 ,isubstack)*thickt
1152 DO it=1,nptt
1153 ipt = ipt_all + it
1154 thk_it(it) = half*thkly*w_gauss(it,nptt)
1155 ratio_thkly(ipt) = thk_it(it)/thickt
1156 IF (ipt == 1) THEN
1157 posly(ipt) = zshift + half*ratio_thkly(ipt)
1158 ELSE
1159 posly(ipt) = posly(ipt-1)
1160 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1161 ENDIF ! IF (IPT 1)
1162 pos_it(it) = posly(ipt)*thk(i)
1163 ENDDO ! nptt
1164 ENDIF ! iint
1165C
1166 !!SCALE = THK(I)/THICKT
1167 DO it=1,nptt
1168 thk_it(it) = scale*thk_it(it)
1169 ms_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
1170 xil = ms_nptt*(area(i)/fac+thk_it(it)*thk_it(it)*one_over_12+
1171 . pos_it(it)*pos_it(it))
1172 xi(i) = xi(i) + xil
1173 ENDDO
1174 ipt_all = ipt_all + nptt
1175 ENDDO ! NLAY
1176 ELSE ! IE
1177 thick_drape = drape(ie)%THICK
1178 scale = thk(i)/thick_drape
1179 DO ilay=1, nlay ! NPT
1180 nptt = elbuf_str%BUFLY(ilay)%NPTT
1181 ip = drape(ie)%INDX_PLY(ilay)
1182 IF(ip > 0) THEN
1183 drape_ply => drape(ie)%DRAPE_PLY(ip)
1184 nslice = drape_ply%NSLICE ! NPTT
1185 i1 = ippid + ilay
1186 i2 = ipmat + ilay
1187 i3 = ipthk + ilay
1188 i4 = ippos + ilay
1189 ipid_ly = stack%IGEO(i1,isubstack)
1190 matly = stack%IGEO(i2,isubstack)
1191 iint = igeo(47,iprop)
1192 IF(iint == 1) THEN
1193 thickt = stack%GEO(1,isubstack)
1194 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1195 ipid_ly = stack%IGEO(i1,isubstack)
1196 matly = stack%IGEO(i2,isubstack)
1197 DO it=1,nptt
1198 ipt = ipt_all + it
1199 thinning = drape_ply%RDRAPE(it,1)
1200 thk_it(it) = thkly*thinning/nptt
1201 ratio_thkly(ipt) = thk_it(it)/thick_drape
1202 IF (ipt == 1) THEN
1203 posly(ipt) = zshift + half*ratio_thkly(ipt)
1204 ELSE
1205 posly(ipt) = posly(ipt-1)
1206 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1207 ENDIF ! IF (IPT == 1)
1208 pos_it(it) = posly(ipt)*thk(i)
1209 ENDDO ! nptt
1210 ELSEIF(iint == 2) THEN
1211 thickt = stack%GEO(1,isubstack)
1212 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1213 DO it=1,nptt
1214 ipt = ipt_all + it
1215 thinning = drape_ply%RDRAPE(it,1)
1216 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
1217 ratio_thkly(ipt) = thk_it(it)/thick_drape
1218 IF (ipt == 1) THEN
1219 posly(ipt) = zshift + half*ratio_thkly(ipt)
1220 ELSE
1221 posly(ipt) = posly(ipt-1)
1222 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1223 ENDIF ! IF (IPT 1)
1224 pos_it(it) = posly(ipt)*thk(i)
1225 ENDDO ! nptt
1226 ENDIF ! iint
1227 !!SCALE = THK(I)/THICK_DRAPE
1228 DO it=1,nptt
1229 thk_it(it) = scale * thk_it(it)
1230 ms_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
1231 xil = ms_nptt*(area(i)/fac+thk_it(it)*thk_it(it)*one_over_12+
1232 . pos_it(it)*pos_it(it))
1233 xi(i) = xi(i) + xil
1234 ENDDO
1235 ELSE ! IP == 0
1236 nptt = elbuf_str%BUFLY(ilay)%NPTT
1237 i1 = ippid + ilay
1238 i2 = ipmat + ilay
1239 i3 = ipthk + ilay
1240 i4 = ippos + ilay
1241 ipid_ly = stack%IGEO(i1,isubstack)
1242 matly = stack%IGEO(i2,isubstack)
1243 iint = igeo(47,iprop)
1244 IF(iint == 1) THEN
1245 thickt = stack%GEO(1,isubstack)
1246 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1247 ipid_ly = stack%IGEO(i1,isubstack)
1248 matly = stack%IGEO(i2,isubstack)
1249 DO it=1,nptt
1250 ipt = ipt_all + it
1251 thk_it(it) = thkly/nptt
1252 ratio_thkly(ipt) = thk_it(it)/thick_drape
1253 IF (ipt == 1) THEN
1254 posly(ipt) = zshift + half*ratio_thkly(ipt)
1255 ELSE
1256 posly(ipt) = posly(ipt-1)
1257 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1258 ENDIF ! IF (IPT == 1)
1259 pos_it(it) = posly(ipt)*thk(i)
1260 ENDDO ! nptt
1261 ELSEIF(iint == 2) THEN
1262 thickt = stack%GEO(1,isubstack)
1263 thkly = stack%GEO(i3 ,isubstack)*thickt ! initial
1264 DO it=1,nptt
1265 ipt = ipt_all + it
1266 thk_it(it) = half*thkly*w_gauss(it,nptt)
1267 ratio_thkly(ipt) = thk_it(it)/thick_drape
1268 IF (ipt == 1) THEN
1269 posly(ipt) = zshift + half*ratio_thkly(ipt)
1270 ELSE
1271 posly(ipt) = posly(ipt-1)
1272 . + half*(ratio_thkly(ipt)+ratio_thkly(ipt-1))
1273 ENDIF ! IF (IPT 1)
1274 pos_it(it) = posly(ipt)*thk(i)
1275 ENDDO ! nptt
1276 ENDIF ! iint
1277 !!SCALE = THK(I)/THICK_DRAPE
1278 DO it=1,nptt
1279 thk_it(it) = scale*thk_it(it)
1280 ms_nptt = pm(89,matly)*thk_it(it)*area(i)*fourth
1281 xil = ms_nptt*(area(i)/fac+thk_it(it)*thk_it(it)*one_over_12+
1282 . pos_it(it)*pos_it(it))
1283 xi(i) = xi(i) + xil
1284 ENDDO
1285 ENDIF ! IP
1286 ipt_all = ipt_all + nptt
1287 ENDDO ! NLAY
1288 ENDIF ! IE
1289 ENDDO ! LFT,LLT
1290 ENDIF ! idrape
1291 ELSE
1292 zshift = geo(199,iprop)
1293 f_offset= zshift*zshift
1294 DO i=lft,llt
1295 xi(i)=ems(i)*(area(i)/fac+thk(i)*thk(i)*(one_over_12+f_offset))
1296 ENDDO
1297 ENDIF
1298C--- check stability condition
1299 IF (igtyp == 11 .OR. igtyp == 16 .OR. igtyp == 17
1300 . .OR. igtyp == 51)THEN
1301 IF(igmat <= 0) THEN !
1302 DO i=lft,llt
1303 ins = ems(i)*(area(i)/fac+thk(i)*thk(i)*one_over_12)
1304 dms = abs(xi(i)-ins)/ins
1305 IF(dms > em02) THEN
1306 idprop = igeo(1,iprop)
1307 CALL fretitl2(titr,
1308 . igeo(npropgi-ltitr+1,iprop),ltitr)
1309 CALL ancmsg(msgid=520,
1310 . msgtype=msgwarning,
1311 . anmode=aninfo_blind_2,
1312 . i1=idprop,
1313 . c1=titr,i2=ix(nixc,nft+i))
1314 ENDIF
1315 ENDDO
1316 ENDIF
1317 ENDIF
1318C----------------------------------------------
1319C-----------------------
1320 DO i=lft,llt
1321 inc(i)=xi(i)
1322 ENDDO
1323C-----------------------
1324 IF(nadmesh==0)THEN
1325 IF (nxref == 0) THEN
1326 DO i=lft,llt
1327 i1 = ix1(i)
1328 i2 = ix2(i)
1329 i3 = ix3(i)
1330 i4 = ix4(i)
1331 ip = ipart(i)
1332C
1333
1334 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1335 partsav(2,ip)=partsav(2,ip) + ems(i)*
1336 . (x(1,i1)+x(1,i2)+x(1,i3)+x(1,i4))
1337 partsav(3,ip)=partsav(3,ip) + ems(i)*
1338 . (x(2,i1)+x(2,i2)+x(2,i3)+x(2,i4))
1339 partsav(4,ip)=partsav(4,ip) + ems(i)*
1340 . (x(3,i1)+x(3,i2)+x(3,i3)+x(3,i4))
1341 xx = (x(1,i1)*x(1,i1)+x(1,i2)*x(1,i2)
1342 . +x(1,i3)*x(1,i3)+x(1,i4)*x(1,i4))
1343 xy = (x(1,i1)*x(2,i1)+x(1,i2)*x(2,i2)
1344 . +x(1,i3)*x(2,i3)+x(1,i4)*x(2,i4))
1345 yy = (x(2,i1)*x(2,i1)+x(2,i2)*x(2,i2)
1346 . +x(2,i3)*x(2,i3)+x(2,i4)*x(2,i4))
1347 yz = (x(2,i1)*x(3,i1)+x(2,i2)*x(3,i2)
1348 . +x(2,i3)*x(3,i3)+x(2,i4)*x(3,i4))
1349 zz = (x(3,i1)*x(3,i1)+x(3,i2)*x(3,i2)
1350 . +x(3,i3)*x(3,i3)+x(3,i4)*x(3,i4))
1351 zx = (x(3,i1)*x(1,i1)+x(3,i2)*x(1,i2)
1352 . +x(3,i3)*x(1,i3)+x(3,i4)*x(1,i4))
1353C
1354 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i) * (yy+zz)
1355 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i) * (zz+xx)
1356 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i) * (xx+yy)
1357 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1358 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1359 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1360C
1361 partsav(11,ip)=partsav(11,ip) + ems(i)*
1362 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1363 partsav(12,ip)=partsav(12,ip) + ems(i)*
1364 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1365 partsav(13,ip)=partsav(13,ip) + ems(i)*
1366 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1367C
1368 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1369 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1370 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1371 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1372 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1373 END DO
1374 ELSE
1375 DO i=lft,llt
1376 i1 = ix1(i)
1377 i2 = ix2(i)
1378 i3 = ix3(i)
1379 i4 = ix4(i)
1380 ip = ipart(i)
1381C
1382 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1383 partsav(2,ip)=partsav(2,ip) + ems(i)*
1384 . (xrefc(1,1,i)+xrefc(2,1,i)+xrefc(3,1,i)+xrefc(4,1,i))
1385 partsav(3,ip)=partsav(3,ip) + ems(i)*
1386 . (xrefc(1,2,i)+xrefc(2,2,i)+xrefc(3,2,i)+xrefc(4,2,i))
1387 partsav(4,ip)=partsav(4,ip) + ems(i)*
1388 . (xrefc(1,3,i)+xrefc(2,3,i)+xrefc(3,3,i)+xrefc(4,3,i))
1389 xx = xrefc(1,1,i)**2 + xrefc(2,1,i)**2
1390 . + xrefc(3,1,i)**2 + xrefc(4,1,i)
1391 xy = xrefc(1,1,i)*xrefc(1,2,i) + xrefc(2,1,i)*xrefc(2,2,i)
1392 . + xrefc(3,1,i)*xrefc(3,2,i) + xrefc(4,1,i)*xrefc(4,2,i)
1393 yy = xrefc(1,2,i)**2 + xrefc(2,2,i)**2
1394 . + xrefc(3,2,i)**2 + xrefc(4,2,i)
1395 yz = xrefc(1,2,i)*xrefc(1,3,i) + xrefc(2,2,i)*xrefc(2,3,i)
1396 . + xrefc(3,2,i)*xrefc(3,3,i) + xrefc(4,2,i)*xrefc(4,3,i)
1397 zz = xrefc(1,3,i)**2 + xrefc(2,3,i)**2
1398 . + xrefc(3,3,i)**2 + xrefc(4,3,i)
1399 zx = xrefc(1,3,i)*xrefc(1,1,i) + xrefc(2,3,i)*xrefc(2,1,i)
1400 . + xrefc(3,3,i)*xrefc(3,1,i) + xrefc(4,3,i)*xrefc(4,1,i)
1401C
1402 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1403 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1404 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1405 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1406 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1407 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1408C
1409 partsav(11,ip)=partsav(11,ip) + ems(i)*
1410 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1411 partsav(12,ip)=partsav(12,ip) + ems(i)*
1412 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1413 partsav(13,ip)=partsav(13,ip) + ems(i)*
1414 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1415C
1416 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1417 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1418 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1419 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1420 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1421 END DO
1422 ENDIF
1423 ELSEIF(istatcnd==0) THEN
1424 IF (nxref == 0) THEN
1425 DO i=lft,llt
1426 n=nft+i
1427 IF(sh4tree(3,n) >= 0)THEN
1428 i1 = ix1(i)
1429 i2 = ix2(i)
1430 i3 = ix3(i)
1431 i4 = ix4(i)
1432 ip = ipart(i)
1433C
1434 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1435 partsav(2,ip)=partsav(2,ip) + ems(i)*
1436 . (x(1,i1)+x(1,i2)+x(1,i3)+x(1,i4))
1437 partsav(3,ip)=partsav(3,ip) + ems(i)*
1438 . (x(2,i1)+x(2,i2)+x(2,i3)+x(2,i4))
1439 partsav(4,ip)=partsav(4,ip) + ems(i)*
1440 . (x(3,i1)+x(3,i2)+x(3,i3)+x(3,i4))
1441 xx = (x(1,i1)*x(1,i1)+x(1,i2)*x(1,i2)
1442 . +x(1,i3)*x(1,i3)+x(1,i4)*x(1,i4))
1443 xy = (x(1,i1)*x(2,i1)+x(1,i2)*x(2,i2)
1444 . +x(1,i3)*x(2,i3)+x(1,i4)*x(2,i4))
1445 yy = (x(2,i1)*x(2,i1)+x(2,i2)*x(2,i2)
1446 . +x(2,i3)*x(2,i3)+x(2,i4)*x(2,i4))
1447 yz = (x(2,i1)*x(3,i1)+x(2,i2)*x(3,i2)
1448 . +x(2,i3)*x(3,i3)+x(2,i4)*x(3,i4))
1449 zz = (x(3,i1)*x(3,i1)+x(3,i2)*x(3,i2)
1450 . +x(3,i3)*x(3,i3)+x(3,i4)*x(3,i4))
1451 zx = (x(3,i1)*x(1,i1)+x(3,i2)*x(1,i2)
1452 . +x(3,i3)*x(1,i3)+x(3,i4)*x(1,i4))
1453C
1454 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1455 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1456 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1457 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1458 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1459 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1460C
1461 partsav(11,ip)=partsav(11,ip) + ems(i)*
1462 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1463 partsav(12,ip)=partsav(12,ip) + ems(i)*
1464 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1465 partsav(13,ip)=partsav(13,ip) + ems(i)*
1466 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1467C
1468 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1469 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1470 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1471 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1472 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1473 END IF
1474 END DO
1475 ELSE
1476 DO i=lft,llt
1477 n=nft+i
1478 IF(sh4tree(3,n) >= 0)THEN
1479 i1 = ix1(i)
1480 i2 = ix2(i)
1481 i3 = ix3(i)
1482 i4 = ix4(i)
1483 ip = ipart(i)
1484C
1485 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1486 partsav(2,ip)=partsav(2,ip) + ems(i)*
1487 . (xrefc(1,1,i)+xrefc(2,1,i)+xrefc(3,1,i)+xrefc(4,1,i))
1488 partsav(3,ip)=partsav(3,ip) + ems(i)*
1489 . (xrefc(1,2,i)+xrefc(2,2,i)+xrefc(3,2,i)+xrefc(4,2,i))
1490 partsav(4,ip)=partsav(4,ip) + ems(i)*
1491 . (xrefc(1,3,i)+xrefc(2,3,i)+xrefc(3,3,i)+xrefc(4,3,i))
1492 xx = xrefc(1,1,i)**2 + xrefc(2,1,i)**2
1493 . + xrefc(3,1,i)**2 + xrefc(4,1,i)
1494 xy = xrefc(1,1,i)*xrefc(1,2,i) + xrefc(2,1,i)*xrefc(2,2,i)
1495 . + xrefc(3,1,i)*xrefc(3,2,i) + xrefc(4,1,i)*xrefc(4,2,i)
1496 yy = xrefc(1,2,i)**2 + xrefc(2,2,i)**2
1497 . + xrefc(3,2,i)**2 + xrefc(4,2,i)
1498 yz = xrefc(1,2,i)*xrefc(1,3,i) + xrefc(2,2,i)*xrefc(2,3,i)
1499 . + xrefc(3,2,i)*xrefc(3,3,i) + xrefc(4,2,i)*xrefc(4,3,i)
1500 zz = xrefc(1,3,i)**2 + xrefc(2,3,i)**2
1501 . + xrefc(3,3,i)**2 + xrefc(4,3,i)
1502 zx = xrefc(1,3,i)*xrefc(1,1,i) + xrefc(2,3,i)*xrefc(2,1,i)
1503 . + xrefc(3,3,i)*xrefc(3,1,i) + xrefc(4,3,i)*xrefc(4,1,i)
1504C
1505 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1506 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1507 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1508 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1509 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1510 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1511C
1512 partsav(11,ip)=partsav(11,ip) + ems(i)*
1513 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1514 partsav(12,ip)=partsav(12,ip) + ems(i)*
1515 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1516 partsav(13,ip)=partsav(13,ip) + ems(i)*
1517 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1518C
1519 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1520 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1521 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1522 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1523 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1524 END IF
1525 END DO
1526 ENDIF
1527 ELSE
1528 IF (nxref == 0) THEN
1529 DO i=lft,llt
1530 n=nft+i
1531 IF(sh4tree(3,n) == 0 .OR. sh4tree(3,n) == -1)THEN
1532 i1 = ix1(i)
1533 i2 = ix2(i)
1534 i3 = ix3(i)
1535 i4 = ix4(i)
1536 ip = ipart(i)
1537C
1538 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1539 partsav(2,ip)=partsav(2,ip) + ems(i)*
1540 . (x(1,i1)+x(1,i2)+x(1,i3)+x(1,i4))
1541 partsav(3,ip)=partsav(3,ip) + ems(i)*
1542 . (x(2,i1)+x(2,i2)+x(2,i3)+x(2,i4))
1543 partsav(4,ip)=partsav(4,ip) + ems(i)*
1544 . (x(3,i1)+x(3,i2)+x(3,i3)+x(3,i4))
1545 xx = (x(1,i1)*x(1,i1)+x(1,i2)*x(1,i2)
1546 . +x(1,i3)*x(1,i3)+x(1,i4)*x(1,i4))
1547 xy = (x(1,i1)*x(2,i1)+x(1,i2)*x(2,i2)
1548 . +x(1,i3)*x(2,i3)+x(1,i4)*x(2,i4))
1549 yy = (x(2,i1)*x(2,i1)+x(2,i2)*x(2,i2)
1550 . +x(2,i3)*x(2,i3)+x(2,i4)*x(2,i4))
1551 yz = (x(2,i1)*x(3,i1)+x(2,i2)*x(3,i2)
1552 . +x(2,i3)*x(3,i3)+x(2,i4)*x(3,i4))
1553 zz = (x(3,i1)*x(3,i1)+x(3,i2)*x(3,i2)
1554 . +x(3,i3)*x(3,i3)+x(3,i4)*x(3,i4))
1555 zx = (x(3,i1)*x(1,i1)+x(3,i2)*x(1,i2)
1556 . +x(3,i3)*x(1,i3)+x(3,i4)*x(1,i4))
1557C
1558 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1559 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1560 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1561 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1562 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1563 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1564C
1565 partsav(11,ip)=partsav(11,ip) + ems(i)*
1566 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1567 partsav(12,ip)=partsav(12,ip) + ems(i)*
1568 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1569 partsav(13,ip)=partsav(13,ip) + ems(i)*
1570 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1571C
1572 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1573 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1574 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1575 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1576 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1577 END IF
1578 END DO
1579 ELSE
1580 DO i=lft,llt
1581 n=nft+i
1582 IF(sh4tree(3,n) == 0 .OR. sh4tree(3,n) == -1)THEN
1583 i1 = ix1(i)
1584 i2 = ix2(i)
1585 i3 = ix3(i)
1586 i4 = ix4(i)
1587 ip = ipart(i)
1588C
1589 partsav(1,ip)=partsav(1,ip) + four*ems(i)
1590 partsav(2,ip)=partsav(2,ip) + ems(i)*
1591 . (xrefc(1,1,i)+xrefc(2,1,i)+xrefc(3,1,i)+xrefc(4,1,i))
1592 partsav(3,ip)=partsav(3,ip) + ems(i)*
1593 . (xrefc(1,2,i)+xrefc(2,2,i)+xrefc(3,2,i)+xrefc(4,2,i))
1594 partsav(4,ip)=partsav(4,ip) + ems(i)*
1595 . (xrefc(1,3,i)+xrefc(2,3,i)+xrefc(3,3,i)+xrefc(4,3,i))
1596 xx = xrefc(1,1,i)**2 + xrefc(2,1,i)**2
1597 . + xrefc(3,1,i)**2 + xrefc(4,1,i)
1598 xy = xrefc(1,1,i)*xrefc(1,2,i) + xrefc(2,1,i)*xrefc(2,2,i)
1599 . + xrefc(3,1,i)*xrefc(3,2,i) + xrefc(4,1,i)*xrefc(4,2,i)
1600 yy = xrefc(1,2,i)**2 + xrefc(2,2,i)**2
1601 . + xrefc(3,2,i)**2 + xrefc(4,2,i)
1602 yz = xrefc(1,2,i)*xrefc(1,3,i) + xrefc(2,2,i)*xrefc(2,3,i)
1603 . + xrefc(3,2,i)*xrefc(3,3,i) + xrefc(4,2,i)*xrefc(4,3,i)
1604 zz = xrefc(1,3,i)**2 + xrefc(2,3,i)**2
1605 . + xrefc(3,3,i)**2 + xrefc(4,3,i)
1606 zx = xrefc(1,3,i)*xrefc(1,1,i) + xrefc(2,3,i)*xrefc(2,1,i)
1607 . + xrefc(3,3,i)*xrefc(3,1,i) + xrefc(4,3,i)*xrefc(4,1,i)
1608C
1609 partsav(5,ip) =partsav(5,ip) + four*xi(i) + ems(i)*(yy+zz)
1610 partsav(6,ip) =partsav(6,ip) + four*xi(i) + ems(i)*(zz+xx)
1611 partsav(7,ip) =partsav(7,ip) + four*xi(i) + ems(i)*(xx+yy)
1612 partsav(8,ip) =partsav(8,ip) - ems(i) * xy
1613 partsav(9,ip) =partsav(9,ip) - ems(i) * yz
1614 partsav(10,ip)=partsav(10,ip) - ems(i) * zx
1615C
1616 partsav(11,ip)=partsav(11,ip) + ems(i)*
1617 . (v(1,i1)+v(1,i2)+v(1,i3)+v(1,i4))
1618 partsav(12,ip)=partsav(12,ip) + ems(i)*
1619 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
1620 partsav(13,ip)=partsav(13,ip) + ems(i)*
1621 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
1622C
1623 partsav(14,ip)=partsav(14,ip) + half * ems(i) *
1624 . (v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
1625 . +v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
1626 . +v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
1627 . +v(1,i4)*v(1,i4)+v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
1628 END IF
1629 END DO
1630 ENDIF
1631 END IF
1632
1633C----------------------------------------------
1634C INITIALISATION DES RIGIDITES NODALES POUR INTERFACES
1635C----------------------------------------------
1636 IF (igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
1637 IF(i7stifs/=0)THEN
1638 DO i=lft,llt
1639 et= stack%PM(2,isubstack)*thk(i)
1640 stc(i)=et
1641 nshnod(ix1(i))=nshnod(ix1(i))+1
1642 nshnod(ix2(i))=nshnod(ix2(i))+1
1643 nshnod(ix3(i))=nshnod(ix3(i))+1
1644 nshnod(ix4(i))=nshnod(ix4(i))+1
1645 ENDDO
1646 ENDIF
1647 ELSEIF(igtyp == 11 .AND. igmat > 0) THEN
1648 ipgmat = 700
1649 IF(i7stifs/=0)THEN
1650 DO i=lft,llt
1651 et=geo(ipgmat +2 ,iprop)*thk(i)
1652 stc(i)=et
1653 nshnod(ix1(i))=nshnod(ix1(i))+1
1654 nshnod(ix2(i))=nshnod(ix2(i))+1
1655 nshnod(ix3(i))=nshnod(ix3(i))+1
1656 nshnod(ix4(i))=nshnod(ix4(i))+1
1657 ENDDO
1658 ENDIF
1659 ELSE ! igmat = 0
1660 IF(i7stifs/=0)THEN
1661 DO i=lft,llt
1662 et=pm(20,imid)*thk(i)
1663 stc(i)=et
1664 nshnod(ix1(i))=nshnod(ix1(i))+1
1665 nshnod(ix2(i))=nshnod(ix2(i))+1
1666 nshnod(ix3(i))=nshnod(ix3(i))+1
1667 nshnod(ix4(i))=nshnod(ix4(i))+1
1668 ENDDO
1669 ENDIF
1670 ENDIF ! igmat
1671 DEALLOCATE(posly,ratio_thkly,thk_it,pos_it )
1672C-----------
1673 RETURN
1674 END
cbainit3
subroutine cbainit3(elbuf_str, ixc, pm, x, geo, xmas, in, nvc, dtelem, igrsh4n, xrefc, nel, ithk, ihbe, igrsh3n, thke, isigsh, sigsh, stifn, stifr, partsav, v, ipart, msc, inc, skew, i8mi, nsigsh, igeo, ipm, iuser, etnod, nshnod, stc, ptshel, bufmat, sh4tree, mcp, mcps, temp, ms_layer, zi_layer, itag, itagel, iparg, ms_layerc, zi_layerc, part_area, cpt_eltens, msz2c, zply, itagn, itage, ixfem, npf, tf, xfem_str, isubstack, stack, rnoise, drape, sh4ang, iddlevel, geo_stack, igeo_stack, strc, perturb, iyldini, ele_area, nloc_dmg, ng, group_param, idrape, drapeg, mat_param, fail_fractal, fail_brokmann, glob_therm)
Definition cbainit3.F:81
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
inirig_mat
subroutine inirig_mat(ixc, ixs, ixtg, ixs10, x, v, pm, geo, ms, in, ptg, msc, mss, mstg, inc, intg, thkc, thkt, partsav, iparts, ipartc, ipartt, veul, dtelem, ihbe, isolnod, nvc, i8mi, msnf, mssf, igeo, etnod, nshnod, stc, sttg, wma, sh4tree, sh3tree, mcp, mcpc, temp, mcps, mssx, mcpsx, ins, stifn, stifr, connec, irig_node, numel, nindx, xrefc, xreftg, xrefs, mssa, sh3trim, isubstack, bufmat, ipm, stack, rnoise, strc, strtg, perturb, nel, group_param, igtyp, defaults, glob_therm)
Definition inirig_mat.F:69
initia
subroutine initia(iparg, elbuf, ms, in, v, x, ixs, ixq, ixc, ixt, ixp, ixr, detonators, geo, pm, rby, npby, lpby, npc, npts, pld, veul, ale_connectivity, skew, fill, ipart, itab, sensors, skvol, ixtg, thk, nloc_dmg, group_param_tab, glob_therm, igrnod, igrsurf, bufsf, vr, bufmat, xlas, las, dtelem, mss, msq, msc, mst, msp, msr, mstg, ptg, inc, nod2eltg, knod2eltg, inp, inr, intg, index, itri, kxx, ixx, xelemwa, iwa, nod2elq, knod2elq, nod2els, knod2els, kxsp, ixsp, nod2sp, ispcond, icode, iskew, iskn, ispsym, xframe, isptag, spbuf, mssx, nsigi, npbyl, lpbyl, rbyl, msnf, mssf, nsigsh, igeo, ipm, nsigs, nsigsph, vns, vnsx, stc, stt, stp, str, sttg, stur, bns, bnsx, volnod, bvolnod, etnod, nshnod, stifint, fxbdep, fxbvit, fxbacc, fxbipm, fxbrpm, fxbelm, fxbsig, fxbmod, ins, ptshel, ptsh3n, ptsol, ptquad, wma, ptsph, fxbnod, mbufel, mdepl, fxani, numel, nsigrs, sh4tree, sh3tree, mcp, temp, imerge2, iadmerge2, slnrbm, nslnrbm, rmstifn, rmstifr, ms_layer, zi_layer, itag, itagel, mcpc, mcptg, xrefc, xreftg, xrefs, mssa, msrt, irbe2, lrbe2, inivol, kvol, nbsubmat, ixs10, ixs16, ixs20, totaddmas, ipmas, stifn, msz2, itagn, sitage, itage, ixr_kj, elbuf_tab, nom_opt, ptr_nopt_rbe2, ptr_nopt_adm, ptr_nopt_fun, sol2sph, irst, sh3trim, xfem_tab, kxig3d, ixig3d, msig3d, knot, nctrlmax, wige, stack, rnoise, drape, sh4ang, sh3ang, geo_stack, igeo_stack, stifintr, strc, strp, strr, strtg, perturb, itagnd, nativ_sms, iloadp, facload, ptspri, nsigbeam, ptbeam, nsigtruss, pttruss, multi_fvm, sigi, sigsh, sigsp, sigsph, sigrs, sigbeam, sigtruss, strsglob, straglob, orthoglob, isigsh, iyldini, ksigsh3, fail_ini, iusolyld, iuser, iddlevel, inimap1d, inimap2d, func2d, fvm_inivel, tagprt_sms, igrbric, igrquad, igrsh4n, igrsh3n, igrpart, totmas, knotlocpc, knotlocel, vnige, bnige, fxbglm, fxbcpm, fxbcps, fxblm, fxbfls, fxbdls, fxb_matrix, fxb_matrix_add, fxb_last_adress, ptr_nopt_fxb, r_skew, knod2el1d, nod2el1d, ebcs_tab, rby_iniaxis, alea, knod2elc, nod2elc, dr, slrbody, drapeg, ipari, intbuf_tab, interfaces, mat_param, npreload_a, preload_a, fail_fractal, fail_brokmann, defaults, ndamp_freq_range, dampr, ibeam_vector, rbeam_vector, ikine)
Definition initia.F:188
max
#define max(a, b)
Definition macros.h:21
drape_mod
Definition drape_mod.F:54
message_mod
Definition message_mod.F:1249
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003
stack_mod
Definition stack_mod.F:72
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:889
fretitl2
subroutine fretitl2(titr, iasc, l)
Definition freform.F:804
starter
program starter
Definition starter.F:39