OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
c3inmas.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com04_c.inc"
#include "param_c.inc"
#include "remesh_c.inc"
#include "scr03_c.inc"
#include "scr12_c.inc"
#include "scr17_c.inc"
#include "vect01_c.inc"
#include "drape_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine c3inmas (x, xreftg, ixtg, geo, pm, ms, tiner, thke, partsav, v, ipart, mstg, intg, ptg, igeo, imat, iprop, area, etnod, nshnod, sttg, sh3tree, mcp, mcptg, temp, sh3trim, isubstack, nlay, elbuf_str, stack, thki, rnoise, drape, perturb, ix1, ix2, ix3, nintemp, x2, x3, y3, idrape, indx)

Function/Subroutine Documentation

◆ c3inmas()

subroutine c3inmas ( x,
xreftg,
integer, dimension(nixtg,*) ixtg,
geo,
pm,
ms,
tiner,
thke,
partsav,
v,
integer, dimension(*) ipart,
mstg,
intg,
ptg,
integer, dimension(npropgi,*) igeo,
integer imat,
integer iprop,
area,
etnod,
integer, dimension(*) nshnod,
sttg,
integer, dimension(ksh3tree,*) sh3tree,
mcp,
mcptg,
temp,
integer, dimension(*) sh3trim,
integer isubstack,
integer nlay,
type(elbuf_struct_), target elbuf_str,
type (stack_ply) stack,
thki,
rnoise,
type (drape_), dimension(numelc_drape + numeltg_drape), target drape,
integer, dimension(nperturb) perturb,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, intent(in) nintemp,
x2,
x3,
y3,
integer idrape,
integer, dimension (stdrape) indx )

Definition at line 38 of file c3inmas.F.

46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE elbufdef_mod
50 USE message_mod
51 USE stack_mod
52 USE drape_mod
54 use element_mod , only : nixtg
55C----------------------------------------------
56C INITIALIZATION OF NODAL MASSES AND INERTIAS
57C-----------------------------------------------
58C I m p l i c i t T y p e s
59C-----------------------------------------------
60#include "implicit_f.inc"
61C-----------------------------------------------
62C G l o b a l P a r a m e t e r s
63C-----------------------------------------------
64#include "mvsiz_p.inc"
65C-----------------------------------------------
66C C o m m o n B l o c k s
67C-----------------------------------------------
68#include "com04_c.inc"
69#include "param_c.inc"
70#include "remesh_c.inc"
71#include "scr03_c.inc"
72#include "scr12_c.inc"
73#include "scr17_c.inc"
74#include "vect01_c.inc"
75#include "drape_c.inc"
76C-----------------------------------------------
77C D u m m y A r g u m e n t s
78C-----------------------------------------------
79 INTEGER IXTG(NIXTG,*),IPART(*),IMAT,IPROP,
80 . IGEO(NPROPGI,*),IGTYP, NSHNOD(*), SH3TREE(KSH3TREE,*),
81 . SH3TRIM(*),ISUBSTACK,NLAY,IDRAPE,
82 . PERTURB(NPERTURB)
83 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ)
84 INTEGER, DIMENSION (STDRAPE) :: INDX
85 INTEGER ,INTENT(IN) :: NINTEMP
87 . x(3,*), geo(npropg,*), pm(npropm,*), ms(*), tiner(*),thke(*),
88 . v(3,*),partsav(20,*),mstg(*),intg(*),ptg(3,*),
89 . etnod(*), sttg(*),mcp(*),mcptg(*),temp(*),xreftg(3,3,*),thki(*),
90 . rnoise(nperturb,*), area(mvsiz)
91 TYPE(ELBUF_STRUCT_), TARGET :: ELBUF_STR
92 TYPE (STACK_PLY) :: STACK
93 TYPE (DRAPE_), DIMENSION(NUMELC_DRAPE + NUMELTG_DRAPE), TARGET :: DRAPE
94C-----------------------------------------------
95C L o c a l V a r i a b l e s
96C-----------------------------------------------
97 INTEGER I,J,K,N,IP,I1,I2,I3,I4,IPTHK,IPMAT,IPPOS,MATLY,IDPROP,
98 . IPID,IPPID,IIGEO,IADI,IADR,IPANG,NTHK,IGMAT,ILAY,NPTT,
99 . IT,MAX_NPTT,IINT,IPID_LY,IPGMAT,IPOS,IPT,IPT_ALL,NSLICE,IE,
100 . NLAY_MAX
101 INTEGER ID
102 CHARACTER(LEN=NCHARTITLE) :: TITR
103 my_real xx,yy,zz,xy,yz,zx,zshift
104 my_real x2(mvsiz), x3(mvsiz), y3(mvsiz),
105 . rho(mvsiz), thk(mvsiz), em(mvsiz), xi(mvsiz),
106 . p1, p2, p3, aa, bb, cc, a2, b2, c2, xin, ems, msl, xil,
107 . mst,dms,ins,posly,thkly,thickt, et,emscp(mvsiz),rhocp(mvsiz),
108 . ems_nptt,ms_nptt,pos_0,thinning,scale,thick_drape
109 my_real, DIMENSION(:) , ALLOCATABLE :: thk_it, pos_it
110 my_real, DIMENSION(:,:), ALLOCATABLE :: laypos, ratio_thkly
111C-----------------------------------------------
112 TYPE (DRAPE_PLY_) , POINTER :: DRAPE_PLY
113C-----------------------------------------------
114 my_real
115 . a_gauss(9,9),w_gauss(9,9)
116C-----------------------------------------------
117 DATA a_gauss /
118 1 0. ,0. ,0. ,
119 1 0. ,0. ,0. ,
120 1 0. ,0. ,0. ,
121 2 -.577350269189626,0.577350269189626,0. ,
122 2 0. ,0. ,0. ,
123 2 0. ,0. ,0. ,
124 3 -.774596669241483,0. ,0.774596669241483,
125 3 0. ,0. ,0. ,
126 3 0. ,0. ,0. ,
127 4 -.861136311594053,-.339981043584856,0.339981043584856,
128 4 0.861136311594053,0. ,0. ,
129 4 0. ,0. ,0. ,
130 5 -.906179845938664,-.538469310105683,0. ,
131 5 0.538469310105683,0.906179845938664,0. ,
132 5 0. ,0. ,0. ,
133 6 -.932469514203152,-.661209386466265,-.238619186083197,
134 6 0.238619186083197,0.661209386466265,0.932469514203152,
135 6 0. ,0. ,0. ,
136 7 -.949107912342759,-.741531185599394,-.405845151377397,
137 7 0. ,0.405845151377397,0.741531185599394,
138 7 0.949107912342759,0. ,0. ,
139 8 -.960289856497536,-.796666477413627,-.525532409916329,
140 8 -.183434642495650,0.183434642495650,0.525532409916329,
141 8 0.796666477413627,0.960289856497536,0. ,
142 9 -.968160239507626,-.836031107326636,-.613371432700590,
143 9 -.324253423403809,0. ,0.324253423403809,
144 9 0.613371432700590,0.836031107326636,0.968160239507626/
145 DATA w_gauss /
146 1 2. ,0. ,0. ,
147 1 0. ,0. ,0. ,
148 1 0. ,0. ,0. ,
149 2 1. ,1. ,0. ,
150 2 0. ,0. ,0. ,
151 2 0. ,0. ,0. ,
152 3 0.555555555555556,0.888888888888889,0.555555555555556,
153 3 0. ,0. ,0. ,
154 3 0. ,0. ,0. ,
155 4 0.347854845137454,0.652145154862546,0.652145154862546,
156 4 0.347854845137454,0. ,0. ,
157 4 0. ,0. ,0. ,
158 5 0.236926885056189,0.478628670499366,0.568888888888889,
159 5 0.478628670499366,0.236926885056189,0. ,
160 5 0. ,0. ,0. ,
161 6 0.171324492379170,0.360761573048139,0.467913934572691,
162 6 0.467913934572691,0.360761573048139,0.171324492379170,
163 6 0. ,0. ,0. ,
164 7 0.129484966168870,0.279705391489277,0.381830050505119,
165 7 0.417959183673469,0.381830050505119,0.279705391489277,
166 7 0.129484966168870,0. ,0. ,
167 8 0.101228536290376,0.222381034453374,0.313706645877887,
168 8 0.362683783378362,0.362683783378362,0.313706645877887,
169 8 0.222381034453374,0.101228536290376,0. ,
170 9 0.081274388361574,0.180648160694857,0.260610696402935,
171 9 0.312347077040003,0.330239355001260,0.312347077040003,
172 9 0.260610696402935,0.180648160694857,0.081274388361574/
173C=======================================================================
174 igtyp = nint(geo(12,iprop))
175 igmat = igeo(98,iprop)
176 max_nptt = 1
177 IF(igtyp == 51 .OR. igtyp == 52) THEN
178 DO ilay=1,nlay ! NPT
179 max_nptt = max(max_nptt ,elbuf_str%BUFLY(ilay)%NPTT)
180 ENDDO
181 ALLOCATE(thk_it(max_nptt),pos_it(max_nptt))
182 ELSE
183 ALLOCATE(thk_it(0),pos_it(0))
184 ENDIF
185 nlay_max = max(nlay, npt)
186 ALLOCATE( laypos(nlay_max*max_nptt,mvsiz),ratio_thkly(nlay_max*max_nptt,mvsiz) )
187 IF ((igtyp == 17 .OR. igtyp == 51) .AND. igmat < 0) THEN
188 DO i=lft,llt
189 IF (ndrape == 0) THEN
190 IF (thke(i)== zero) THEN
191 thk(i) = stack%GEO(1,isubstack)
192 thke(i)= thk(i)
193 ELSE
194 thk(i)=thke(i)
195 ENDIF
196 thki(i) = stack%GEO(1,isubstack)
197 ELSEIF (ndrape > 0) THEN
198 thk(i) = thke(i)
199 thki(i)= thk(i)
200 ENDIF ! IF (NDRAPE == 0)
201c
202 IF (nperturb /= 0) THEN
203 DO j=1,nperturb
204 IF (perturb(j) == 1 .AND.
205 . rnoise(j,numelc+i+nft) /= zero) THEN
206 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
207 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
208 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
209 ENDIF
210 ENDDO
211 ENDIF
212c
213 rho(i)= pm(89,imat)
214 rhocp(i) = pm(69,imat)
215 IF(thk(i)==zero.AND.igtyp/=0)THEN
216 CALL ancmsg(msgid=495,
217 . msgtype=msgerror,
218 . anmode=aninfo_blind_1,
219 . i1=ixtg(6,i))
220 ENDIF
221 ENDDO
222 ELSEIF (igtyp == 52 .OR.
223 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
224 DO i=lft,llt
225 IF (ndrape == 0) THEN
226 IF (thke(i)== zero) THEN
227 thk(i) = stack%GEO(1,isubstack)
228 thke(i)= thk(i)
229 ELSE
230 thk(i)=thke(i)
231 ENDIF
232 thki(i) = stack%GEO(1,isubstack)
233 ELSEIF (ndrape > 0) THEN
234 thk(i) = thke(i)
235 thki(i)= thk(i)
236 ENDIF ! IF (NDRAPE == 0)
237c
238 IF (nperturb /= 0) THEN
239 DO j=1,nperturb
240 IF (perturb(j) == 1 .AND.
241 . rnoise(j,numelc+i+nft) /= zero) THEN
242 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
243 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
244 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
245 ENDIF
246 ENDDO
247 ENDIF
248c
249 rho(i) = stack%PM(1,isubstack)
250 rhocp(i) = stack%PM(15,isubstack)* rho(i)/stack%PM(14,isubstack)
251 IF (thk(i) == zero .AND. igtyp /= 0) THEN
252 CALL ancmsg(msgid=495,
253 . msgtype=msgerror,
254 . anmode=aninfo_blind_1,
255 . i1=ixtg(6,i))
256 ENDIF
257 ENDDO
258 ELSE
259 DO i=lft,llt
260 IF(thke(i)==zero)THEN
261 thk(i) =geo(1,iprop)
262 thke(i)=thk(i)
263 ELSE
264 thk(i)=thke(i)
265 ENDIF
266 thki(i) = thk(i)
267c
268 IF (nperturb /= 0) THEN
269 DO j=1,nperturb
270 IF (perturb(j) == 1 .AND.
271 . rnoise(j,numelc+i+nft) /= zero) THEN
272 thk(i) = thk(i) * rnoise(j,numelc+i+nft)
273 thke(i) = thke(i) * rnoise(j,numelc+i+nft)
274 thki(i) = thki(i) * rnoise(j,numelc+i+nft)
275 ENDIF
276 ENDDO
277 ENDIF
278c
279 rho(i)= pm(89,imat)
280 rhocp(i) = pm(69,imat)
281 IF(thk(i)==zero.AND.igtyp/=0)THEN
282 CALL ancmsg(msgid=495,
283 . msgtype=msgerror,
284 . anmode=aninfo_blind_1,
285 . i1=ixtg(6,i))
286 ENDIF
287 ENDDO
288 ENDIF
289C----------------------------------------------
290C MASSES ELEMENTAIRES
291C----------------------------------------------
292 IF (igtyp==11 .OR. igtyp==16) THEN
293 ipthk = 300
294 ippos = 400
295 ipmat = 100
296 DO i=lft,llt
297 em(i)=zero
298 DO n=1,npt
299 i1=ipthk+n
300 i2=ipmat+n
301 thickt = geo(200,iprop)
302 thkly = geo(i1,iprop)*thk(i)
303 matly = igeo(i2,iprop)
304 em(i) = em(i) + pm(89,matly)*thkly*area(i)
305 ENDDO
306C--- Test stability condition for mass
307 mst = rho(i)*thk(i)*area(i)
308 dms = abs(em(i)-mst)/mst
309C! for igty=11 ---> IGMAT=-1 (old G-mat), 1 for new g-mat
310C! for igtyp=16 --->IGMAT= 0
311C! for igtyp=17 --->IGMAT= 0
312 IF(igmat <= 0) THEN
313 IF (dms > em02) THEN
314 idprop = nint(geo(40,iprop))
315 CALL fretitl2(titr,
316 . igeo(npropgi-ltitr+1,iprop),ltitr)
317 CALL ancmsg(msgid=519,
318 . msgtype=msgwarning,
319 . anmode=aninfo_blind_2,
320 . i1=igeo(1,iprop),
321 . c1=titr,i2=ixtg(nixtg,nft+i),
322 . r1=em(i),r2=mst)
323 ENDIF
324 ENDIF
325 ENDDO
326 ELSEIF (igtyp == 17) THEN
327 ippid = 2
328 ipmat = ippid + npt
329 ipang = 1
330 ipthk = ipang + npt
331 ippos = ipthk + npt
332 nthk = ippos + npt
333 IF(idrape == 0) THEN
334 DO i=lft,llt
335 em(i)=zero
336 DO n=1,npt
337 i1=ipthk + n
338 i2=ippid + n
339 i3=ipmat + n
340 thickt = stack%GEO(1,isubstack)
341 thkly = stack%GEO(i1,isubstack)*thk(i)
342 matly = stack%IGEO(i3,isubstack)
343 em(i) = em(i) + pm(89,matly)*thkly*area(i)
344 ENDDO
345C--- Test stability condition for mass
346 mst = rho(i)*thk(i)*area(i)
347 dms = abs(em(i)-mst)/mst
348 IF(igmat < 0) THEN
349 IF (dms > em02) THEN
350 idprop = nint(geo(40,iprop))
351 CALL fretitl2(titr,
352 . igeo(npropgi-ltitr+1,iprop),ltitr)
353 CALL ancmsg(msgid=519,
354 . msgtype=msgwarning,
355 . anmode=aninfo_blind_2,
356 . i1=igeo(1,iprop),
357 . c1=titr,i2=ixtg(nixtg,nft+i),
358 . r1=em(i),r2=mst)
359 ENDIF
360 ENDIF
361 ENDDO
362 ELSE ! IDRAPE > 0
363 DO i=lft,llt
364 em(i)=zero
365 ie = indx(nft + i)
366 IF(ie == 0) THEN
367 DO n=1,npt
368 i1=ipthk + n
369 i2=ippid + n
370 i3=ipmat + n
371 thickt = stack%GEO(1,isubstack)
372 thkly = stack%GEO(i1,isubstack)*thk(i)
373 matly = stack%IGEO(i3,isubstack)
374 em(i) = em(i) + pm(89,matly)*thkly*area(i)
375 ENDDO
376 ELSE
377 thick_drape = drape(ie)%THICK
378 scale = thk(i)/thick_drape
379 DO n=1,npt
380 i1=ipthk + n
381 i2=ippid + n
382 i3=ipmat + n
383 ip = drape(ie)%INDX_PLY(n)
384 IF(ip > 0) THEN
385 drape_ply => drape(ie)%DRAPE_PLY(ip)
386 nslice = drape_ply%NSLICE
387 thinning = drape_ply%RDRAPE(1,1) ! one slice by layer
388 thickt = stack%GEO(1,isubstack)
389 thkly = stack%GEO(i1,isubstack)*thickt
390 matly = stack%IGEO(i3,isubstack)
391 thkly = thkly*thinning*scale
392 matly = stack%IGEO(i3,isubstack)
393 em(i) = em(i) + pm(89,matly)*thkly*area(i)
394 ELSE
395 thickt = stack%GEO(1,isubstack)
396 thkly = stack%GEO(i1,isubstack)*thickt*scale
397 matly = stack%IGEO(i3,isubstack)
398 matly = stack%IGEO(i3,isubstack)
399 em(i) = em(i) + pm(89,matly)*thkly*area(i)
400 ENDIF
401 ENDDO
402 ENDIF ! IE
403C--- Test stability condition for mass
404 mst = rho(i)*thk(i)*area(i)
405 dms = abs(em(i)-mst)/mst
406 IF(igmat < 0) THEN
407 IF (dms > em02) THEN
408 idprop = nint(geo(40,iprop))
409 CALL fretitl2(titr,
410 . igeo(npropgi-ltitr+1,iprop),ltitr)
411 CALL ancmsg(msgid=519,
412 . msgtype=msgwarning,
413 . anmode=aninfo_blind_2,
414 . i1=igeo(1,iprop),
415 . c1=titr,i2=ixtg(nixtg,nft+i),
416 . r1=em(i),r2=mst)
417 ENDIF
418 ENDIF
419 ENDDO
420 ENDIF ! DRAPE
421 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
422 ipang = 1
423 ippid = 2
424 ipmat = ippid + nlay ! NPT
425 ipthk = ipang + nlay ! NPT
426 ippos = ipthk + nlay ! NPT
427 nthk = ippos + nlay ! NPT
428 !
429 ipos = igeo(99,iprop)
430 thickt = stack%GEO(1,isubstack)
431 zshift = geo(199,iprop)
432 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
433 IF(idrape == 0) THEN
434 DO i=lft,llt
435 scale = thk(i)/thickt
436 em(i)=zero
437 ipt_all = 0
438 DO ilay=1,nlay ! NPT
439 nptt = elbuf_str%BUFLY(ilay)%NPTT
440 i1=ipthk + ilay
441 i2=ippid + ilay
442 i3=ipmat + ilay
443 thickt = stack%GEO(1,isubstack)
444 thkly = stack%GEO(i1,isubstack)*thickt
445 ipid_ly= stack%IGEO(i2,isubstack)
446 matly = stack%IGEO(i3,isubstack)
447 ipid = stack%IGEO(ippid,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
448 iint = igeo(47,iprop)
449 IF (iint == 1) THEN ! uniform distribution - by default
450 DO it=1,nptt
451 thk_it(it) = thkly/nptt ! equally distribution of NPTT through layer
452 ENDDO
453 ELSEIF (iint == 2) THEN ! Gauss Distribution
454 DO it=1,nptt
455 thk_it(it) = half*thkly*w_gauss(it,nptt)
456 ENDDO
457 ENDIF
458 !SCALE = THK(I)/THICKT
459 DO it=1,nptt
460 thk_it(it) = thk_it(it)*scale
461 ems_nptt = pm(89,matly)*thk_it(it)*area(i)
462 em(i) = em(i) + ems_nptt
463 ENDDO
464C----
465 ENDDO
466C--- Test stability condition for mass
467 mst = rho(i)*thk(i)*area(i)
468 dms = abs(em(i)-mst)/mst
469 IF(igtyp == 51 .AND.igmat < 0) THEN
470 IF (dms > em02) THEN
471 idprop = nint(geo(40,iprop))
472 CALL fretitl2(titr,
473 . igeo(npropgi-ltitr+1,iprop),ltitr)
474 CALL ancmsg(msgid=519,
475 . msgtype=msgwarning,
476 . anmode=aninfo_blind_2,
477 . i1=igeo(1,iprop),
478 . c1=titr,i2=ixtg(nixtg,nft+i),
479 . r1=em(i),r2=mst)
480 ENDIF
481 ENDIF
482 ENDDO
483 ELSE ! IDRAPE > 0
484 thickt = stack%GEO(1,isubstack)
485 DO i=lft,llt
486 em(i)=zero
487 ie = indx(nft + i)
488 ipt_all = 0
489 scale = thk(i)/thickt
490 IF(ie == 0) THEN
491 DO ilay=1,nlay ! NPT
492 nptt = elbuf_str%BUFLY(ilay)%NPTT
493 i1=ipthk + ilay
494 i2=ippid + ilay
495 i3=ipmat + ilay
496 ipid = stack%IGEO(i2,isubstack) ! or IPID = IX(NIXC-1,NFT+I)
497 iint = igeo(47,iprop)
498 thickt = stack%GEO(1,isubstack)
499 IF(iint == 1) THEN
500 thickt = stack%GEO(1,isubstack)
501 thkly = stack%GEO(i1,isubstack)*thickt
502 ipid_ly= stack%IGEO(i2,isubstack)
503 matly = stack%IGEO(i3,isubstack)
504 DO it=1,nptt
505 thk_it(it) = thkly/nptt ! uniform distribution of NPTT through layer
506 ENDDO
507 ELSEIF(iint == 2)THEN
508 thickt = stack%GEO(1,isubstack)
509 thkly = stack%GEO(i1,isubstack)*thickt
510 ipid_ly= stack%IGEO(i2,isubstack)
511 matly = stack%IGEO(i3,isubstack)
512 DO it=1,nptt
513 thk_it(it) = half*thkly*w_gauss(it,nptt) ! G Gauss distribution of NPTT through layer
514 ENDDO
515 ENDIF ! iint
516C
517 DO it=1,nptt
518 thk_it(it) = scale*thk_it(it)
519 ems_nptt = pm(89,matly)*thk_it(it)*area(i)
520 em(i) = em(i) + ems_nptt
521 ENDDO
522 ENDDO ! NLAY
523 ELSE ! IE > 0
524 ipt_all = 0
525 thick_drape = drape(ie)%THICK
526 scale = thk(i)/thick_drape
527 DO ilay=1,nlay ! NPT
528 nptt = elbuf_str%BUFLY(ilay)%NPTT
529 i1=ipthk + ilay
530 i2=ippid + ilay
531 i3=ipmat + ilay
532 ipid = stack%IGEO(i2,isubstack) ! or ipid = ix(nixc-1,nft+i)
533 iint = igeo(47,iprop)
534 thickt = stack%GEO(1,isubstack)
535 ip = drape(ie)%INDX_PLY(ilay)
536 IF(ip == 0) THEN
537 thickt = stack%GEO(1,isubstack)
538 thkly = stack%GEO(i1,isubstack)*thickt
539 ipid_ly= stack%IGEO(i2,isubstack)
540 matly = stack%IGEO(i3,isubstack)
541 IF(iint == 1) THEN
542 DO it=1,nptt
543 thk_it(it) = thkly/nptt
544 ENDDO ! nptt
545 ELSEIF(iint == 2)THEN
546 DO it=1,nptt
547 thk_it(it) = half*thkly*w_gauss(it,nptt) ! Gauss distribution of NPTT through laye r
548 ENDDO ! nptt
549 ENDIF ! iint
550 ELSE
551 drape_ply => drape(ie)%DRAPE_PLY(ip)
552 nslice = drape_ply%NSLICE !NPTT
553 thickt = stack%GEO(1,isubstack)
554 thkly = stack%GEO(i1,isubstack)*thickt
555 ipid_ly= stack%IGEO(i2,isubstack)
556 matly = stack%IGEO(i3,isubstack)
557 IF(iint == 1) THEN
558 DO it=1,nptt
559 thinning = drape_ply%RDRAPE(it,1)
560 thk_it(it) = thkly*thinning/nptt
561 ENDDO ! nptt
562 ELSEIF(iint == 2)THEN
563 DO it=1,nptt
564 thinning = drape_ply%RDRAPE(it,1)
565 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning ! Gauss distribution of NPTT through layer
566 ENDDO ! nptt
567 ENDIF ! iint
568 ENDIF
569 DO it=1,nptt
570 thk_it(it) = scale*thk_it(it)
571 ems_nptt = pm(89,matly)*thk_it(it)*area(i)
572 em(i) = em(i) + ems_nptt
573 ENDDO
574 ENDDO ! NLAY
575 ENDIF !IE
576C--- Test stability condition for mass
577 mst = rho(i)*thk(i)*area(i)
578 dms = abs(em(i)-mst)/mst
579 IF(igtyp == 51 .AND.igmat < 0) THEN
580 IF (dms > em02) THEN
581 idprop = nint(geo(40,iprop))
582 CALL fretitl2(titr,
583 . igeo(npropgi-ltitr+1,iprop),ltitr)
584 CALL ancmsg(msgid=519,
585 . msgtype=msgwarning,
586 . anmode=aninfo_blind_2,
587 . i1=igeo(1,iprop),
588 . c1=titr,i2=ixtg(nixtg,nft+i),
589 . r1=em(i),r2=mst)
590 ENDIF
591 ENDIF
592 ENDDO ! LFT,LLT
593 ENDIF ! IDRAPE
594 ELSE
595 DO i=lft,llt
596 em(i) = rho(i)*thk(i)*area(i)
597 ENDDO
598 DO i=lft,llt
599 emscp(i)=rhocp(i)*thk(i)*area(i)
600 ENDDO
601 ENDIF
602C----------------------------------------------
603C INERTIES ELEMENTAIRES
604C----------------------------------------------
605 IF (igtyp==11 .OR. igtyp==16) THEN
606 DO i=lft,llt
607 xi(i)=zero
608 DO n=1,npt
609 i1=ipthk+n
610 i2=ipmat+n
611 i3=ippos+n
612 thickt= geo(200,iprop)
613 thkly = geo(i1,iprop)*thk(i)
614 posly = geo(i3,iprop)*thk(i)
615 matly = igeo(i2,iprop)
616 msl = pm(89,matly)*thkly*area(i)
617 xil = msl*(area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
618 xi(i) = xi(i) + xil
619 ENDDO
620C--- Test stability condition for inertia
621 IF(igmat <= 0) THEN !
622 ins = em(i)*(area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
623 dms = abs(xi(i)-ins)/ins
624 IF(dms > em02) THEN
625 idprop = nint(geo(40,iprop))
626 CALL fretitl2(titr,
627 . igeo(npropgi-ltitr+1,iprop),ltitr)
628 CALL ancmsg(msgid=520,
629 . msgtype=msgwarning,
630 . anmode=aninfo_blind_2,
631 . i1=igeo(1,iprop),
632 . c1=titr,i2=ixtg(nixtg,nft+i))
633 ENDIF
634 ENDIF
635 ENDDO
636 ELSEIF (igtyp == 17) THEN
637 ipos = igeo(99,iprop)
638 thickt = stack%GEO(1,isubstack)
639 zshift = geo(199,iprop)
640 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
641 IF(idrape == 0 ) THEN
642 DO i=lft,llt
643 xi(i)=zero
644 DO n=1,npt
645 i1 =ippid + n
646 i2 =ipmat + n
647 i3 =ipthk + n
648 i4 =ippos + n
649 thickt= stack%GEO(1,isubstack)
650 thkly = stack%GEO(i3,isubstack)*thk(i)
651 posly = stack%GEO(i4,isubstack)*thk(i)
652 matly = stack%IGEO(i2,isubstack)
653 msl = pm(89,matly)*thkly*area(i)
654 xil = msl*(area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
655 xi(i) = xi(i) + xil
656 ENDDO
657C--- Test stability condition for inertia
658 ins = em(i)*(area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
659 dms = abs(xi(i)-ins)/ins
660
661 IF(igmat <= 0) THEN !
662 IF(dms > em02) THEN
663 idprop = nint(geo(40,iprop))
664 CALL fretitl2(titr,
665 . igeo(npropgi-ltitr+1,iprop),ltitr)
666 CALL ancmsg(msgid=520,
667 . msgtype=msgwarning,
668 . anmode=aninfo_blind_2,
669 . i1=igeo(1,iprop),
670 . c1=titr,i2=ixtg(nixtg,nft+i))
671 ENDIF
672 ENDIF
673 ENDDO
674 ELSE !IDRAPE > 0
675 DO i=lft,llt
676 xi(i)=zero
677 ie = indx(nft + i)
678 IF(ie == 0) THEN
679 DO n=1,npt
680 i1 =ippid + n
681 i2 =ipmat + n
682 i3 =ipthk + n
683 i4 =ippos + n
684 thickt= stack%GEO(1,isubstack)
685 thkly = stack%GEO(i3,isubstack)*thk(i)
686 posly = stack%GEO(i4,isubstack)*thk(i)
687 matly = stack%IGEO(i2,isubstack)
688 msl = pm(89,matly)*thkly*area(i)
689 xil = msl*(area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
690 xi(i) = xi(i) + xil
691 ENDDO
692 ELSE ! IE > 0
693 thick_drape = drape(ie)%THICK
694 scale = thk(i)/thick_drape
695 DO n=1,npt
696 ip = drape(ie)%INDX_PLY(n)
697 i1 =ippid + n
698 i2 =ipmat + n
699 i3 =ipthk + n
700 i4 =ippos + n
701 IF(ip == 0 ) THEN
702 thkly = stack%GEO(i3,isubstack)*thickt
703 matly = stack%IGEO(i2,isubstack)
704 ratio_thkly(n,i) = thkly/ thick_drape
705 ELSE
706 drape_ply => drape(ie)%DRAPE_PLY(ip)
707 thinning = drape_ply%RDRAPE(1,1)
708 thkly = stack%GEO(i3,isubstack)*thickt
709 thkly = thkly*thinning
710 matly = stack%IGEO(i2,isubstack)
711 ratio_thkly(n,i) = thkly/thick_drape
712 ENDIF
713 IF (n == 1) THEN
714 laypos(n,i) = zshift + half*ratio_thkly(n,i)
715 ELSE
716 laypos(n,i) = laypos(n-1,i) + half*(ratio_thkly(n,i)+ratio_thkly(n-1,i))
717 ENDIF ! IF (N == 1)
718 posly = laypos(n,i)*thk(i)
719 thkly = thkly*scale
720 msl = pm(89,matly)*thkly*area(i)
721 xil = msl*(area(i)/nine_over_2+thkly*thkly*one_over_12+posly*posly)
722 xi(i) = xi(i) + xil
723 ENDDO
724 ENDIF ! IE
725C--- Test stability condition for inertia
726 ins = em(i)*(area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
727 dms = abs(xi(i)-ins)/ins
728
729 IF(igmat <= 0) THEN !
730 IF(dms > em02) THEN
731 idprop = nint(geo(40,iprop))
732 CALL fretitl2(titr,
733 . igeo(npropgi-ltitr+1,iprop),ltitr)
734 CALL ancmsg(msgid=520,
735 . msgtype=msgwarning,
736 . anmode=aninfo_blind_2,
737 . i1=igeo(1,iprop),
738 . c1=titr,i2=ixtg(nixtg,nft+i))
739 ENDIF
740 ENDIF
741 ENDDO
742 ENDIF ! IDRAPE
743 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
744 ipang = 1
745 ippid = 2
746 ipmat = ippid + nlay ! NPT
747 ipthk = ipang + nlay ! NPT
748 ippos = ipthk + nlay ! NPT
749 nthk = ippos + nlay ! NPT
750 !!
751
752 ipos = igeo(99,iprop)
753 thickt = stack%GEO(1,isubstack)
754 zshift = geo(199,iprop)
755 IF(ipos == 2 ) zshift = - zshift /max(thickt,em20)
756 IF(idrape == 0) THEN
757 thickt = stack%GEO(1,isubstack)
758 DO i=lft,llt
759 xi(i)=zero
760 ipt_all = 0
761 scale = thk(i)/thickt
762 DO ilay=1,nlay ! NPT
763 nptt = elbuf_str%BUFLY(ilay)%NPTT
764 i1 =ippid + ilay
765 i2 =ipmat + ilay
766 i3 =ipthk + ilay
767 i4 =ippos + ilay
768 ipid_ly = stack%IGEO(i1,isubstack)
769 matly = stack%IGEO(i2,isubstack)
770 iint = igeo(47,iprop)
771 IF(iint == 1) THEN
772 thickt = stack%GEO(1,isubstack)
773 thkly = stack%GEO(i3,isubstack)*thickt
774 DO it=1,nptt
775 ipt = ipt_all + it
776 thk_it(it) = thkly/nptt
777 ratio_thkly(ipt,i) = thk_it(it)/thickt
778 IF (ipt == 1) THEN
779 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
780 ELSE
781 laypos(ipt,i) = laypos(ipt-1,i) +
782 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
783 ENDIF ! IF (N == 1)
784 pos_it(it) = laypos(ipt,i)*thk(i)
785 ENDDO
786 ELSEIF(iint == 2) THEN
787 thickt = stack%GEO(1,isubstack)
788 thkly = stack%GEO(i3,isubstack)*thickt
789 DO it=1,nptt
790 ipt = ipt_all + it
791 thk_it(it) = half*thkly*w_gauss(it,nptt)
792 ratio_thkly(ipt,i) = thk_it(it)/thickt
793 IF (ipt == 1) THEN
794 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
795 ELSE
796 laypos(ipt,i) = laypos(ipt-1,i) +
797 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
798 ENDIF ! IF (N == 1)
799 pos_it(it) = laypos(ipt,i)*thk(i)
800 ENDDO
801 ENDIF ! iit
802 !SCALE = THK(I)/THICKT
803 DO it=1,nptt
804 ms_nptt = pm(89,matly)*thk_it(it)*area(i)
805 thk_it(it) = thk_it(it)*scale
806 xil = ms_nptt*(area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
807 . pos_it(it)*pos_it(it))
808 xi(i) = xi(i) + xil
809 ENDDO
810 ipt_all = ipt_all + nptt
811 ENDDO ! NLAY
812C Test stability condition for inertia
813 ins = em(i)*(area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
814 dms = abs(xi(i)-ins)/ins
815 IF(igtyp == 51 .AND. igmat < 0) THEN
816 IF (dms > em02) THEN
817 idprop = nint(geo(40,iprop))
818 CALL fretitl2(titr,
819 . igeo(npropgi-ltitr+1,iprop),ltitr)
820 CALL ancmsg(msgid=520,
821 . msgtype=msgwarning,
822 . anmode=aninfo_blind_2,
823 . i1=igeo(1,iprop),
824 . c1=titr,i2=ixtg(nixtg,nft+i))
825 ENDIF
826 ENDIF
827 ENDDO ! LFT,LLT
828 ELSE ! IDRAPE > 0
829 thickt = stack%GEO(1,isubstack)
830 DO i=lft,llt
831 xi(i)=zero
832 ipt_all = 0
833 ie = indx(nft + i)
834 scale = thk(i)/thickt
835 IF(ie == 0) THEN
836 DO ilay=1,nlay ! NPT
837 nptt = elbuf_str%BUFLY(ilay)%NPTT
838 i1 =ippid + ilay
839 i2 =ipmat + ilay
840 i3 =ipthk + ilay
841 i4 =ippos + ilay
842 ipid_ly = stack%IGEO(i1,isubstack)
843 matly = stack%IGEO(i2,isubstack)
844 iint = igeo(47,iprop)
845 IF(iint == 1) THEN
846 thickt = stack%GEO(1,isubstack)
847 thkly = stack%GEO(i3,isubstack)*thickt
848 DO it=1,nptt
849 ipt = ipt_all + it
850 thk_it(it) = thkly/nptt
851 ratio_thkly(ipt,i) = thk_it(it)/thickt
852 IF (ipt == 1) THEN
853 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
854 ELSE
855 laypos(ipt,i) = laypos(ipt-1,i) +
856 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
857 ENDIF ! IF (N == 1)
858 pos_it(it) = laypos(ipt,i)*thk(i)
859 ENDDO
860 ELSEIF(iint == 2) THEN
861 thickt = stack%GEO(1,isubstack)
862 thkly = stack%GEO(i3,isubstack)*thickt
863 DO it=1,nptt
864 ipt = ipt_all + it
865 thk_it(it) = half*thkly*w_gauss(it,nptt)
866 ratio_thkly(ipt,i) = thk_it(it)/thickt
867 IF (ipt == 1) THEN
868 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
869 ELSE
870 laypos(ipt,i) = laypos(ipt-1,i) +
871 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
872 ENDIF ! IF (N == 1)
873 pos_it(it) = laypos(ipt,i)*thk(i)
874 ENDDO
875 ENDIF ! iit
876 DO it=1,nptt
877 thk_it(it) = scale*thk_it(it)
878 ms_nptt = pm(89,matly)*thk_it(it)*area(i)
879 xil = ms_nptt*(area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
880 . pos_it(it)*pos_it(it))
881 xi(i) = xi(i) + xil
882 ENDDO
883 ipt_all = ipt_all + nptt
884 ENDDO ! NLAY
885 ELSE ! IE > 0
886 ipt_all = 0
887 thick_drape = drape(ie)%THICK
888 scale = thk(i)/thick_drape
889 DO ilay=1,nlay ! NPT
890 ip = drape(ie)%INDX_PLY(ilay)
891 nptt = elbuf_str%BUFLY(ilay)%NPTT
892 i1 =ippid + ilay
893 i2 =ipmat + ilay
894 i3 =ipthk + ilay
895 i4 =ippos + ilay
896 ipid_ly = stack%IGEO(i1,isubstack)
897 matly = stack%IGEO(i2,isubstack)
898 iint = igeo(47,iprop)
899 IF(ip > 0) THEN
900 drape_ply => drape(ie)%DRAPE_PLY(ip)
901 nslice = drape_ply%NSLICE ! NPTT
902 IF(iint == 1) THEN
903 thickt = stack%GEO(1,isubstack)
904 thkly = stack%GEO(i3,isubstack)*thickt
905 DO it=1,nptt
906 ipt = ipt_all + it
907 thinning = drape_ply%RDRAPE(it,1)
908 thk_it(it) = thkly*thinning/nptt
909 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
910 IF (ipt == 1) THEN
911 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
912 ELSE
913 laypos(ipt,i) = laypos(ipt-1,i) +
914 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
915 ENDIF ! IF (N == 1)
916 pos_it(it) = laypos(ipt,i)*thk(i)
917 ENDDO
918 ELSEIF(iint == 2) THEN
919 thickt = stack%GEO(1,isubstack)
920 thkly = stack%GEO(i3,isubstack)*thickt
921 DO it=1,nptt
922 ipt = ipt_all + it
923 thinning = drape_ply%RDRAPE(it,1)
924 thk_it(it) = half*thkly*w_gauss(it,nptt)*thinning
925 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
926 IF (ipt == 1) THEN
927 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
928 ELSE
929 laypos(ipt,i) = laypos(ipt-1,i) +
930 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
931 ENDIF ! IF (N == 1)
932 pos_it(it) = laypos(ipt,i)*thk(i)
933 ENDDO
934 ENDIF ! iit
935 ELSE ! IP == 0
936 IF(iint == 1) THEN
937 thickt = stack%GEO(1,isubstack)
938 thkly = stack%GEO(i3,isubstack)*thickt
939 DO it=1,nptt
940 ipt = ipt_all + it
941 thk_it(it) = thkly/nptt
942 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
943 IF (ipt == 1) THEN
944 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
945 ELSE
946 laypos(ipt,i) = laypos(ipt-1,i) +
947 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
948 ENDIF ! IF (N == 1)
949 pos_it(it) = laypos(ipt,i)*thk(i)
950 ENDDO
951 ELSEIF(iint == 2) THEN
952 thickt = stack%GEO(1,isubstack)
953 thkly = stack%GEO(i3,isubstack)*thickt
954 DO it=1,nptt
955 ipt = ipt_all + it
956 thk_it(it) = half*thkly*w_gauss(it,nptt)
957 ratio_thkly(ipt,i) = thk_it(it)/thick_drape
958 IF (ipt == 1) THEN
959 laypos(ipt,i) = zshift + half*ratio_thkly(ipt,i)
960 ELSE
961 laypos(ipt,i) = laypos(ipt-1,i) +
962 . half*(ratio_thkly(ipt,i)+ratio_thkly(ipt-1,i))
963 ENDIF ! IF (n == 1)
964 pos_it(it) = laypos(ipt,i)*thk(i)
965 ENDDO
966 ENDIF ! iit
967 ENDIF
968 DO it=1,nptt
969 thk_it(it) = scale*thk_it(it)
970 ms_nptt = pm(89,matly)*thk_it(it)*area(i)
971 xil = ms_nptt*(area(i)/nine_over_2+thk_it(it)*thk_it(it)*one_over_12 +
972 . pos_it(it)*pos_it(it))
973 xi(i) = xi(i) + xil
974 ENDDO
975 ipt_all = ipt_all + nptt
976 ENDDO ! NLAY
977 ENDIF ! ie
978C Test stability condition for inertia
979 ins = em(i)*(area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
980 dms = abs(xi(i)-ins)/ins
981 IF(igtyp == 51 .AND. igmat < 0) THEN
982 IF (dms > em02) THEN
983 idprop = nint(geo(40,iprop))
984 CALL fretitl2(titr,
985 . igeo(npropgi-ltitr+1,iprop),ltitr)
986 CALL ancmsg(msgid=520,
987 . msgtype=msgwarning,
988 . anmode=aninfo_blind_2,
989 . i1=igeo(1,iprop),
990 . c1=titr,i2=ixtg(nixtg,nft+i))
991 ENDIF
992 ENDIF
993 ENDDO ! LFT,LLT
994 ENDIF ! IDRAPE
995 ELSE
996 DO i=lft,llt
997C--- Exact: XIN = EMS*(AREA(I)/(6.*SQRT(3.))+THK(I)*THK(I)/12.)
998C--- Compatibility with shells 4 nodes:
999 xi(i)=em(i)*(area(i)/nine_over_2+thk(i)*thk(i)*one_over_12)
1000 ENDDO
1001 ENDIF
1002C---------------------------------------------------
1003C INITIALIZATION OF NODAL MASSES AND INERTIAS
1004C---------------------------------------------------
1005C mass and inertia spmd parith/on + adaptive meshing
1006 IF(jthe == 0 ) THEN
1007 DO i=lft,llt
1008 mstg(i)=em(i)
1009 intg(i)=xi(i)
1010 END DO
1011 ELSE
1012 DO i=lft,llt
1013 mstg(i)=em(i)
1014 intg(i)=xi(i)
1015 mcptg(i) = emscp(i)
1016 END DO
1017 ENDIF
1018C
1019 DO i=lft,llt
1020 aa = x2(i)
1021 a2 = aa**2
1022 b2 = (x2(i)-x3(i))**2 + y3(i)**2
1023 bb = sqrt(b2)
1024 c2 = x3(i)**2 + y3(i)**2
1025 cc = sqrt(c2)
1026 p1 = acos((a2 + c2 - b2)/(two * aa * cc)) / pi
1027 p2 = acos((a2 + b2 - c2)/(two * aa * bb)) / pi
1028 p3 = acos((b2 + c2 - a2)/(two * bb * cc)) / pi
1029 ptg(1,i)=p1
1030 ptg(2,i)=p2
1031 ptg(3,i)=p3
1032 IF ( ( ((a2 + c2 - b2)/(two * aa * cc) <= -one) .OR.
1033 . ( (a2 + c2 - b2)/(two * aa * cc) >= one) .OR.
1034 . ( (a2 + b2 - c2)/(two * aa * bb) <= -one) .OR.
1035 . ( (a2 + b2 - c2)/(two * aa * bb) >= one) .OR.
1036 . ( (b2 + c2 - a2)/(two * bb * cc) <= -one) .OR.
1037 . ( (b2 + c2 - a2)/(two * bb * cc) >= one)) .AND.
1038 . lsh3trim /= 0 ) THEN
1039 IF (sh3trim(nft+i) /= -1) CALL ancmsg(msgid=750,
1040 . msgtype=msgerror,
1041 . anmode=aninfo,
1042 . i1=ixtg(6,nft+i))
1043 ELSEIF ( ( ((a2 + c2 - b2)/(two * aa * cc) <= -one) .OR.
1044 . ( (a2 + c2 - b2)/(two * aa * cc) >= one) .OR.
1045 . ( (a2 + b2 - c2)/(two * aa * bb) <= -one) .OR.
1046 . ( (a2 + b2 - c2)/(two * aa * bb) >= one) .OR.
1047 . ( (b2 + c2 - a2)/(two * bb * cc) <= -one) .OR.
1048 . ( (b2 + c2 - a2)/(two * bb * cc) >= one)) .AND.
1049 . lsh3trim == 0 ) THEN
1050 CALL ancmsg(msgid=750,
1051 . msgtype=msgerror,
1052 . anmode=aninfo,
1053 . i1=ixtg(6,nft+i))
1054 ENDIF
1055 END DO
1056C
1057C-----
1058 IF(nadmesh==0)THEN
1059 IF (nxref == 0) THEN
1060 DO i=lft,llt
1061 p1 = ptg(1,i)
1062 p2 = ptg(2,i)
1063 p3 = ptg(3,i)
1064 ems = em(i)
1065 xin = xi(i)
1066C
1067 i1 = ix1(i)
1068 i2 = ix2(i)
1069 i3 = ix3(i)
1070C
1071 ip=ipart(i)
1072 partsav(1,ip)=partsav(1,ip) + ems
1073 partsav(2,ip)=partsav(2,ip) + ems*
1074 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1075 partsav(3,ip)=partsav(3,ip) + ems*
1076 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1077 partsav(4,ip)=partsav(4,ip) + ems*
1078 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1079 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1080 . +x(1,i3)*x(1,i3)*p3)
1081 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1082 . +x(1,i3)*x(2,i3)*p3)
1083 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1084 . +x(2,i3)*x(2,i3)*p3)
1085 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1086 . +x(2,i3)*x(3,i3)*p3)
1087 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1088 . +x(3,i3)*x(3,i3)*p3)
1089 zx = (x(3,i1)*x(1,i1)*p1+x(3,i2)*x(1,i2)*p2
1090 . +x(3,i3)*x(1,i3)*p3)
1091 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1092 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1093 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1094 partsav(8,ip) =partsav(8,ip) - ems * xy
1095 partsav(9,ip) =partsav(9,ip) - ems * yz
1096 partsav(10,ip)=partsav(10,ip) - ems * zx
1097C
1098 partsav(11,ip)=partsav(11,ip) + ems*
1099 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1100 partsav(12,ip)=partsav(12,ip) + ems*
1101 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1102 partsav(13,ip)=partsav(13,ip) + ems*
1103 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1104 partsav(14,ip)=partsav(14,ip) + half * ems *
1105 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1106 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1107 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1108 ENDDO
1109 ELSE
1110 DO i=lft,llt
1111 p1 = ptg(1,i)
1112 p2 = ptg(2,i)
1113 p3 = ptg(3,i)
1114 ems = em(i)
1115 xin = xi(i)
1116C
1117 i1 = ix1(i)
1118 i2 = ix2(i)
1119 i3 = ix3(i)
1120 ip=ipart(i)
1121C
1122 partsav(1,ip)=partsav(1,ip) + ems
1123 partsav(2,ip)=partsav(2,ip) + ems*
1124 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)+p3*xreftg(3,1,i))
1125 partsav(3,ip)=partsav(3,ip) + ems*
1126 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1127 partsav(4,ip)=partsav(4,ip) + ems*
1128 . (p1*xreftg(1,3,i)+p2*xreftg(2,3,i)+p3*xreftg(3,3,i))
1129 xx = p1*xreftg(1,1,i)**2
1130 . + p2*xreftg(2,1,i)**2
1131 . + p3*xreftg(3,1,i)**2
1132 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1133 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1134 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1135 yy = p1*xreftg(1,2,i)**2
1136 . + p2*xreftg(2,2,i)**2
1137 . + p3*xreftg(3,2,i)**2
1138 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1139 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1140 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1141 zz = p1*xreftg(1,3,i)**2
1142 . + p2*xreftg(2,3,i)**2
1143 . + p3*xreftg(3,3,i)**2
1144 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1145 . + p2*xreftg(2,3,i)*xreftg(2,1,i)
1146 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1147C
1148 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1149 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1150 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1151 partsav(8,ip) =partsav(8,ip) - ems * xy
1152 partsav(9,ip) =partsav(9,ip) - ems * yz
1153 partsav(10,ip)=partsav(10,ip) - ems * zx
1154C
1155 partsav(11,ip)=partsav(11,ip) + ems*
1156 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1157 partsav(12,ip)=partsav(12,ip) + ems*
1158 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1159 partsav(13,ip)=partsav(13,ip) + ems*
1160 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1161 partsav(14,ip)=partsav(14,ip) + half * ems *
1162 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1163 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1164 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1165 ENDDO
1166 ENDIF
1167 ELSEIF (istatcnd==0)THEN
1168 IF (nxref == 0) THEN
1169 DO i=lft,llt
1170 n=nft+i
1171 IF(sh3tree(3,n) >= 0)THEN
1172 p1 = ptg(1,i)
1173 p2 = ptg(2,i)
1174 p3 = ptg(3,i)
1175 ems = em(i)
1176 xin = xi(i)
1177C
1178 i1 = ix1(i)
1179 i2 = ix2(i)
1180 i3 = ix3(i)
1181 ip = ipart(i)
1182C
1183 partsav(1,ip)=partsav(1,ip) + ems
1184 partsav(2,ip)=partsav(2,ip) + ems*
1185 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1186 partsav(3,ip)=partsav(3,ip) + ems*
1187 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1188 partsav(4,ip)=partsav(4,ip) + ems*
1189 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1190 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1191 . +x(1,i3)*x(1,i3)*p3)
1192 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1193 . +x(1,i3)*x(2,i3)*p3)
1194 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1195 . +x(2,i3)*x(2,i3)*p3)
1196 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1197 . +x(2,i3)*x(3,i3)*p3)
1198 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1199 . +x(3,i3)*x(3,i3)*p3)
1200 zx = (x(3,i1)*x(1,i1)*p1+x(3,i2)*x(1,i2)*p2
1201 . +x(3,i3)*x(1,i3)*p3)
1202C
1203 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1204 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1205 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1206 partsav(8,ip) =partsav(8,ip) - ems * xy
1207 partsav(9,ip) =partsav(9,ip) - ems * yz
1208 partsav(10,ip)=partsav(10,ip) - ems * zx
1209C
1210 partsav(11,ip)=partsav(11,ip) + ems*
1211 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1212 partsav(12,ip)=partsav(12,ip) + ems*
1213 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1214 partsav(13,ip)=partsav(13,ip) + ems*
1215 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1216 partsav(14,ip)=partsav(14,ip) + half * ems *
1217 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1218 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1219 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1220 END IF
1221 END DO
1222 ELSE
1223 DO i=lft,llt
1224 n=nft+i
1225 IF (sh3tree(3,n) >= 0)THEN
1226 p1 = ptg(1,i)
1227 p2 = ptg(2,i)
1228 p3 = ptg(3,i)
1229 ems = em(i)
1230 xin = xi(i)
1231C
1232 i1 = ix1(i)
1233 i2 = ix2(i)
1234 i3 = ix3(i)
1235 ip = ipart(i)
1236C
1237 partsav(1,ip)=partsav(1,ip) + ems
1238 partsav(2,ip)=partsav(2,ip) + ems*
1239 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)+p3*xreftg(3,1,i))
1240 partsav(3,ip)=partsav(3,ip) + ems*
1241 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1242 partsav(4,ip)=partsav(4,ip) + ems*
1243 . (p1*xreftg(1,3,i)+p2*xreftg(2,3,i)+p3*xreftg(3,3,i))
1244 xx = p1*xreftg(1,1,i)**2
1245 . + p2*xreftg(2,1,i)**2
1246 . + p3*xreftg(3,1,i)**2
1247 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1248 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1249 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1250 yy = p1*xreftg(1,2,i)**2
1251 . + p2*xreftg(2,2,i)**2
1252 . + p3*xreftg(3,2,i)**2
1253 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1254 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1255 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1256 zz = p1*xreftg(1,3,i)**2
1257 . + p2*xreftg(2,3,i)**2
1258 . + p3*xreftg(3,3,i)**2
1259 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1260 . + p2*xreftg(2,3,i)*xreftg(2,1,i)
1261 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1262C
1263 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1264 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1265 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1266 partsav(8,ip) =partsav(8,ip) - ems * xy
1267 partsav(9,ip) =partsav(9,ip) - ems * yz
1268 partsav(10,ip)=partsav(10,ip) - ems * zx
1269C
1270 partsav(11,ip)=partsav(11,ip) + ems*
1271 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1272 partsav(12,ip)=partsav(12,ip) + ems*
1273 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1274 partsav(13,ip)=partsav(13,ip) + ems*
1275 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1276 partsav(14,ip)=partsav(14,ip) + half * ems *
1277 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1278 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1279 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1280 END IF
1281 END DO
1282 ENDIF
1283 ELSE
1284 IF (nxref == 0) THEN
1285 DO i=lft,llt
1286 n=nft+i
1287 IF(sh3tree(3,n) == 0 .OR. sh3tree(3,n) == -1 )THEN
1288 p1 = ptg(1,i)
1289 p2 = ptg(2,i)
1290 p3 = ptg(3,i)
1291 ems = em(i)
1292 xin = xi(i)
1293C
1294 i1 = ix1(i)
1295 i2 = ix2(i)
1296 i3 = ix3(i)
1297 ip = ipart(i)
1298C
1299 partsav(1,ip)=partsav(1,ip) + ems
1300 partsav(2,ip)=partsav(2,ip) + ems*
1301 . (x(1,i1)*p1+x(1,i2)*p2+x(1,i3)*p3)
1302 partsav(3,ip)=partsav(3,ip) + ems*
1303 . (x(2,i1)*p1+x(2,i2)*p2+x(2,i3)*p3)
1304 partsav(4,ip)=partsav(4,ip) + ems*
1305 . (x(3,i1)*p1+x(3,i2)*p2+x(3,i3)*p3)
1306 xx = (x(1,i1)*x(1,i1)*p1+x(1,i2)*x(1,i2)*p2
1307 . +x(1,i3)*x(1,i3)*p3)
1308 xy = (x(1,i1)*x(2,i1)*p1+x(1,i2)*x(2,i2)*p2
1309 . +x(1,i3)*x(2,i3)*p3)
1310 yy = (x(2,i1)*x(2,i1)*p1+x(2,i2)*x(2,i2)*p2
1311 . +x(2,i3)*x(2,i3)*p3)
1312 yz = (x(2,i1)*x(3,i1)*p1+x(2,i2)*x(3,i2)*p2
1313 . +x(2,i3)*x(3,i3)*p3)
1314 zz = (x(3,i1)*x(3,i1)*p1+x(3,i2)*x(3,i2)*p2
1315 . +x(3,i3)*x(3,i3)*p3)
1316 zx = (x(3,i1)*x(1,i1)*p1+x(3,i2)*x(1,i2)*p2
1317 . +x(3,i3)*x(1,i3)*p3)
1318C
1319 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1320 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1321 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1322 partsav(8,ip) =partsav(8,ip) - ems * xy
1323 partsav(9,ip) =partsav(9,ip) - ems * yz
1324 partsav(10,ip)=partsav(10,ip) - ems * zx
1325C
1326 partsav(11,ip)=partsav(11,ip) + ems*
1327 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1328 partsav(12,ip)=partsav(12,ip) + ems*
1329 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1330 partsav(13,ip)=partsav(13,ip) + ems*
1331 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1332 partsav(14,ip)=partsav(14,ip) + half * ems *
1333 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1334 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1335 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1336 END IF
1337 END DO
1338 ELSE
1339 DO i=lft,llt
1340 n=nft+i
1341 IF(sh3tree(3,n) == 0 .OR. sh3tree(3,n) == -1 )THEN
1342 p1 = ptg(1,i)
1343 p2 = ptg(2,i)
1344 p3 = ptg(3,i)
1345 ems = em(i)
1346 xin = xi(i)
1347C
1348 i1 = ix1(i)
1349 i2 = ix2(i)
1350 i3 = ix3(i)
1351 ip = ipart(i)
1352C
1353 partsav(1,ip)=partsav(1,ip) + ems
1354 partsav(2,ip)=partsav(2,ip) + ems*
1355 . (p1*xreftg(1,1,i)+p2*xreftg(2,1,i)+p3*xreftg(3,1,i))
1356 partsav(3,ip)=partsav(3,ip) + ems*
1357 . (p1*xreftg(1,2,i)+p2*xreftg(2,2,i)+p3*xreftg(3,2,i))
1358 partsav(4,ip)=partsav(4,ip) + ems*
1359 . (p1*xreftg(1,3,i)+p2*xreftg(2,3,i)+p3*xreftg(3,3,i))
1360 xx = p1*xreftg(1,1,i)**2
1361 . + p2*xreftg(2,1,i)**2
1362 . + p3*xreftg(3,1,i)**2
1363 xy = p1*xreftg(1,1,i)*xreftg(1,2,i)
1364 . + p2*xreftg(2,1,i)*xreftg(2,2,i)
1365 . + p3*xreftg(3,1,i)*xreftg(3,2,i)
1366 yy = p1*xreftg(1,2,i)**2
1367 . + p2*xreftg(2,2,i)**2
1368 . + p3*xreftg(3,2,i)**2
1369 yz = p1*xreftg(1,2,i)*xreftg(1,3,i)
1370 . + p2*xreftg(2,2,i)*xreftg(2,3,i)
1371 . + p3*xreftg(3,2,i)*xreftg(3,3,i)
1372 zz = p1*xreftg(1,3,i)**2
1373 . + p2*xreftg(2,3,i)**2
1374 . + p3*xreftg(3,3,i)**2
1375 zx = p1*xreftg(1,3,i)*xreftg(1,1,i)
1376 . + p2*xreftg(2,3,i)*xreftg(2,1,i)
1377 . + p3*xreftg(3,3,i)*xreftg(3,1,i)
1378C
1379 partsav(5,ip) =partsav(5,ip) + xin + ems * (yy+zz)
1380 partsav(6,ip) =partsav(6,ip) + xin + ems * (zz+xx)
1381 partsav(7,ip) =partsav(7,ip) + xin + ems * (xx+yy)
1382 partsav(8,ip) =partsav(8,ip) - ems * xy
1383 partsav(9,ip) =partsav(9,ip) - ems * yz
1384 partsav(10,ip)=partsav(10,ip) - ems * zx
1385C
1386 partsav(11,ip)=partsav(11,ip) + ems*
1387 . (v(1,i1)*p1+v(1,i2)*p2+v(1,i3)*p3)
1388 partsav(12,ip)=partsav(12,ip) + ems*
1389 . (v(2,i1)*p1+v(2,i2)*p2+v(2,i3)*p3)
1390 partsav(13,ip)=partsav(13,ip) + ems*
1391 . (v(3,i1)*p1+v(3,i2)*p2+v(3,i3)*p3)
1392 partsav(14,ip)=partsav(14,ip) + half * ems *
1393 . ((v(1,i1)*v(1,i1)+v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1))*p1
1394 . +(v(1,i2)*v(1,i2)+v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2))*p2
1395 . +(v(1,i3)*v(1,i3)+v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3))*p3)
1396 END IF
1397 END DO
1398 ENDIF
1399 END IF
1400C-----
1401 IF(jthe > 0 ) THEN
1402 IF(nintemp > 0 ) THEN
1403 DO i= lft,llt
1404 IF(temp(ix1(i))== zero) temp(ix1(i)) = pm(79,imat)
1405 IF(temp(ix2(i))== zero) temp(ix2(i)) = pm(79,imat)
1406 IF(temp(ix3(i))== zero) temp(ix3(i)) = pm(79,imat)
1407 ENDDO
1408 ELSE
1409 DO i=lft,llt
1410 temp(ix1(i))= pm(79,imat)
1411 temp(ix2(i))= pm(79,imat)
1412 temp(ix3(i))= pm(79,imat)
1413 ENDDO
1414 ENDIF
1415 ENDIF
1416C----------------------------------------------
1417C INITIALIZATION OF NODAL RIGIDITIES FOR INTERFACES
1418C----------------------------------------------
1419 IF (igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))THEN
1420 IF(i7stifs/=0)THEN
1421 DO i=lft,llt
1422 et=stack%PM(2,isubstack)*thk(i)
1423 sttg(i)=et
1424 nshnod(ix1(i))=nshnod(ix1(i))+1
1425 nshnod(ix2(i))=nshnod(ix2(i))+1
1426 nshnod(ix3(i))=nshnod(ix3(i))+1
1427 ENDDO
1428 ENDIF
1429 ELSEIF(igtyp == 11 .AND. igmat > 0) THEN
1430 ipgmat = 700
1431 IF(i7stifs/=0)THEN
1432 DO i=lft,llt
1433 et=geo(ipgmat +2 ,iprop)*thk(i)
1434 sttg(i)=et
1435 nshnod(ix1(i))=nshnod(ix1(i))+1
1436 nshnod(ix2(i))=nshnod(ix2(i))+1
1437 nshnod(ix3(i))=nshnod(ix3(i))+1
1438 ENDDO
1439 ENDIF
1440 ElSE
1441 IF(i7stifs/=0)THEN
1442 DO i=lft,llt
1443 et=pm(20,imat)*thk(i)
1444 sttg(i)=et
1445 nshnod(ix1(i))=nshnod(ix1(i))+1
1446 nshnod(ix2(i))=nshnod(ix2(i))+1
1447 nshnod(ix3(i))=nshnod(ix3(i))+1
1448 ENDDO
1449 ENDIF
1450 ENDIF ! IGMAT
1451 DEALLOCATE(thk_it,pos_it, laypos,ratio_thkly )
1452C-----------
1453 RETURN
#define my_real
Definition cppsort.cpp:32
subroutine area(d1, x, x2, y, y2, eint, stif0)
#define max(a, b)
Definition macros.h:21
integer, parameter nchartitle
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:895
subroutine fretitl2(titr, iasc, l)
Definition freform.F:799