Télécharger rayn.eso

Retour à la liste

Numérotation des lignes :

  1. C RAYN SOURCE PV 17/09/29 21:15:46 9578
  2. SUBROUTINE RAYN
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5.  
  6. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  7. C C
  8. C But : Calculer unbe matrice de conductivité due aux échanges C
  9. C par rayonnement C
  10. C C
  11. C Syntaxe : CND1 = RAYN MO1 CHTEM CHRAY (CONSSB) ; C
  12. C C
  13. C MO1 : modèle C
  14. C CHTEM : temperature moyenne (MCHAML) C
  15. C CHRAY : matrice de rayonnement (MCHAML) C
  16. C CONSSB : constante de Stefan (REEL) C
  17. C C
  18. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  19.  
  20. -INC CCOPTIO
  21. -INC CCREEL
  22. *-
  23. -INC SMCOORD
  24. -INC SMELEME
  25. -INC SMCHAML
  26. -INC SMRIGID
  27. -INC SMMODEL
  28. -INC SMINTE
  29.  
  30. POINTEUR CHAMTM.MCHELM,CHAMR.MCHELM
  31. POINTEUR MELBMA.MELEME,MELCOQ.MELEME
  32. POINTEUR MAPBMA.MELEME,MAPCOQ.MELEME
  33. POINTEUR MAPCOM.MELEME
  34.  
  35. SEGMENT ,INFOEL
  36. LOGICAL KCOQ(N1),KQUAD(N1)
  37. ENDSEGMENT
  38.  
  39. SEGMENT,GL2LOC
  40. INTEGER IEQUIV(3,IECART)
  41. ENDSEGMENT
  42.  
  43. SEGMENT/INFO/(INFELL(JG)),INFO1.INFO,INFO2.INFO
  44. SEGMENT LESINF
  45. INTEGER LINFO(NBSZEL)
  46. ENDSEGMENT
  47.  
  48. SEGMENT LESAIJ
  49. REAL*8 AIJ(NBNNMX,NBELMX,NBSZEL)
  50. ENDSEGMENT
  51.  
  52. DIMENSION T3(2)
  53. C ... Attention !!! 20 correspond ici au nombre maxi de noeuds par
  54. C élément
  55. C ... XEL - coordonnées d'origine, XCH - coordonnées choisies pour
  56. C VPAST (ou VPAST2),
  57. C XTR - coordonnées transformées par VCORL1 ...
  58. DIMENSION XEL(3,20),XCH(3,3),XTR(3,20)
  59. C ... BPSS - matrice de transformation ...
  60. DIMENSION BPSS(3,3)
  61. DIMENSION SHP(6,20)
  62.  
  63. C-------- Fin des déclarations -------------------------------
  64.  
  65. CHARACTER*8 TYPE
  66. DATA INTTYP /3/
  67. DATA CONSSB /5.67D-8/
  68.  
  69. C----------- Fin des DATA ------------------------------------
  70.  
  71. C ... Lecture ...
  72.  
  73. TYPE='MMODEL '
  74. CALL LIROBJ(TYPE,IRET,1,IRETOU)
  75. IF(IRETOU.EQ.0) THEN
  76. MOTERR(1:8)=TYPE
  77. CALL ERREUR(37)
  78. RETURN
  79. ENDIF
  80. IF(IRETOU.EQ.1) MMODEL=IRET
  81. SEGACT MMODEL
  82.  
  83. TYPE='MCHAML '
  84. CALL LIROBJ(TYPE,IRET,1,IRETOU)
  85. IF(IRETOU.EQ.0) THEN
  86. MOTERR(1:8)=TYPE
  87. CALL ERREUR(37)
  88. RETURN
  89. ENDIF
  90. IF(IRETOU.EQ.1) MCHEL1=IRET
  91.  
  92. TYPE='MCHAML '
  93. CALL LIROBJ(TYPE,IRET,1,IRETOU)
  94. IF(IRETOU.EQ.0) THEN
  95. MOTERR(1:8)=TYPE
  96. CALL ERREUR(37)
  97. RETURN
  98. ENDIF
  99. IF(IRETOU.EQ.1) MCHEL2=IRET
  100.  
  101. SEGACT, MCHEL1, MCHEL2
  102. IF((MCHEL1.TITCHE).EQ.'MATRICE DE RAYONNEMENT') THEN
  103. CHAMR = MCHEL1
  104. CHAMTM= MCHEL2
  105. ELSEIF((MCHEL2.TITCHE).EQ.'MATRICE DE RAYONNEMENT') THEN
  106. CHAMR = MCHEL2
  107. CHAMTM= MCHEL1
  108. ELSE
  109. CALL ERREUR(21)
  110. RETURN
  111. ENDIF
  112. SEGDES, MCHEL1, MCHEL2
  113.  
  114. CALL LIRREE(XVAL,0,IRETOU)
  115. IF(IRETOU.EQ.1) CONSSB=XVAL
  116.  
  117. C ... Verification que les supports du MMODEL et des
  118. C MCHAML sont bien les mêmes ...
  119.  
  120. CALL RAYN1(MMODEL,CHAMR)
  121. CALL RAYN1(MMODEL,CHAMTM)
  122.  
  123. C CHAMTM = pointeur vers le MCHAML de température moyenne
  124. C CHAMR = pointeur vers le MCHAML de matrices de rayonnement
  125. C MMODEL = pointeur vers le modèle
  126.  
  127. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  128. C ... Recherche du support géométrique de la rigidité ... C
  129. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  130.  
  131. SEGACT,MMODEL
  132.  
  133. C ... On créé un MELEME contenant tous les maillages-supports
  134. C du CHAMTM ...
  135. NBSZEL=KMODEL(/1)
  136. N1=NBSZEL
  137. SEGINI,INFOEL
  138. DO 1000 I=1,NBSZEL
  139. IMODEL=KMODEL(I)
  140. SEGACT,IMODEL
  141. IF(FORMOD(1).NE.'THERMIQUE'.or.matmod(2).NE.'RAYONNEMENT')THEN
  142. CALL ERREUR(21)
  143. RETURN
  144. ENDIF
  145. MELEME=IMAMOD
  146. SEGACT,MELEME
  147.  
  148. IF (IDIM.EQ.3) THEN
  149. C
  150. IF ((NEFMOD.EQ.4).OR.(NEFMOD.EQ.8)) THEN
  151. C TRI3 ou QUA4
  152. KQUAD(I)=.FALSE.
  153. KCOQ(I) =.FALSE.
  154. ELSEIF ((NEFMOD.EQ.27).OR.(NEFMOD.EQ.49)) THEN
  155. C COQ3 ou COQ4
  156. KQUAD(I)=.FALSE.
  157. KCOQ(I) =.TRUE.
  158. ELSEIF ((NEFMOD.EQ.6).OR.(NEFMOD.EQ.10)) THEN
  159. C TRI6 ou QUA8
  160. KQUAD(I)=.TRUE.
  161. KCOQ(I) =.FALSE.
  162. ELSEIF ((NEFMOD.EQ.56).OR.(NEFMOD.EQ.41)) THEN
  163. C COQ6 ou COQ8
  164. KQUAD(I)=.TRUE.
  165. KCOQ(I) =.TRUE.
  166. ELSE
  167. CALL ERREUR(16)
  168. RETURN
  169. ENDIF
  170. C
  171. ELSEIF (IDIM.EQ.2) THEN
  172. C
  173. IF (NEFMOD.EQ.2) THEN
  174. C SEG2
  175. KQUAD(I)=.FALSE.
  176. KCOQ(I) =.FALSE.
  177. ELSEIF (NEFMOD.EQ.44) THEN
  178. C COQ2
  179. KQUAD(I)=.FALSE.
  180. KCOQ(I) =.TRUE.
  181. ELSEIF (NEFMOD.EQ.3) THEN
  182. C SEG3
  183. KQUAD(I)=.TRUE.
  184. KCOQ(I) =.FALSE.
  185. ELSE
  186. CALL ERREUR(16)
  187. RETURN
  188. ENDIF
  189. C
  190. ELSE
  191. CALL ERREUR(14)
  192. RETURN
  193. ENDIF
  194.  
  195. C SEGDES,IMODEL
  196. 1000 CONTINUE
  197.  
  198. NBMCOQ=0
  199. DO 1010 I=1,NBSZEL
  200. IF(KCOQ(I)) NBMCOQ=NBMCOQ+1
  201. 1010 CONTINUE
  202. NBMBMA=NBSZEL-NBMCOQ
  203.  
  204. NBNN=0
  205. NBELEM=0
  206. NBREF=0
  207. NBSOUS=NBMBMA
  208. IF(NBMBMA.GT.0) THEN
  209. SEGINI MELBMA
  210. J=0
  211. DO 1020 I=1,NBSZEL
  212. IF(.NOT.KCOQ(I)) THEN
  213. J=J+1
  214. IMODEL=KMODEL(I)
  215. MELBMA.LISOUS(J)=IMAMOD
  216. ENDIF
  217. 1020 CONTINUE
  218. ELSE
  219. MELBMA=0
  220. ENDIF
  221. IF(IIMPI.GE.4) write(*,*) ' MELBMA',MELBMA
  222.  
  223. NBNN=0
  224. NBELEM=0
  225. NBREF=0
  226. NBSOUS=NBMCOQ
  227. IF(NBMCOQ.GT.0) THEN
  228. SEGINI MELCOQ
  229. J=0
  230. DO 1030 I=1,NBSZEL
  231. IF(KCOQ(I)) THEN
  232. J=J+1
  233. IMODEL=KMODEL(I)
  234. MELCOQ.LISOUS(J)=IMAMOD
  235. ENDIF
  236. 1030 CONTINUE
  237. ELSE
  238. MELCOQ=0
  239. ENDIF
  240. IF(IIMPI.GE.4) write(*,*) ' MELCOQ',MELCOQ
  241.  
  242. C ... On les transforme en maillages composés de POI1 ...
  243. IF(NBMBMA.GT.0) THEN
  244. MAPBMA=MELBMA
  245. CALL CHANGE(MAPBMA,1)
  246.  
  247. c CALL ECROBJ('MAILLAGE',MAPBMA)
  248. c CALL NBNO
  249. c CALL LIRENT(NBNBMA,1,IRETOU)
  250. c IF(IERR.NE.0) GOTO 9999
  251. NBNBMA=MAPBMA.NUM(/2)
  252. ELSE
  253. MAPBMA=0
  254. NBNBMA=0
  255. ENDIF
  256.  
  257. IF(NBMCOQ.GT.0) THEN
  258. MAPCOQ=MELCOQ
  259. CALL CHANGE(MAPCOQ,1)
  260.  
  261. c CALL ECROBJ('MAILLAGE',MAPCOQ)
  262. c CALL NBNO
  263. c CALL LIRENT(NBNCOQ,1,IRETOU)
  264. c IF(IERR.NE.0) GOTO 9999
  265. NBNCOQ=MAPCOQ.NUM(/2)
  266. ELSE
  267. MAPCOQ=0
  268. NBNCOQ=0
  269. ENDIF
  270.  
  271. C ... On détruit les maillages chapeaux ...
  272. SEGSUP,MELBMA,MELCOQ
  273.  
  274. C ... On va vérifier si les deux ensembles de POI1 sont distincts ...
  275. C ... En calculant le nombre de noeuds de la différence symétrique
  276. C des deux maillages ...
  277.  
  278. IF(NBNBMA*NBNCOQ.NE.0) THEN
  279. CALL ECROBJ('MAILLAGE',MAPBMA)
  280. CALL ECROBJ('MAILLAGE',MAPCOQ)
  281. CALL PRDIFF
  282. CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
  283. IF(IIMPI.GE.4) write(*,*) ' IPT1 ',IPT1
  284. IF(IERR.NE.0) GOTO 9999
  285.  
  286. C!! segprt,IPT1
  287.  
  288. c CALL ECROBJ('MAILLAGE',IPT1)
  289. c CALL NBNO
  290. c CALL LIRENT(NBNDS,1,IRETOU)
  291. c IF(IERR.NE.0) GOTO 9999
  292. SEGACT IPT1
  293. NBNDS=IPT1.NUM(/2)
  294. IF(IIMPI.GE.4) THEN
  295. write(*,*) ' NBNDS ',NBNDS
  296. write(*,*) 'NBNBMA+NBNCOQ ',NBNBMA+NBNCOQ
  297. ENDIF
  298.  
  299. SEGSUP,IPT1
  300.  
  301. C ... Et en le comparant avec la somme des nombres de noeuds ...
  302. IF(NBNBMA+NBNCOQ.NE.NBNDS) THEN
  303.  
  304. CALL ECROBJ('MAILLAGE',MAPBMA)
  305. CALL ECROBJ('MAILLAGE',MAPCOQ)
  306. CALL INTERS
  307. CALL LIROBJ('MAILLAGE',MAPCOM,1,IRETOU)
  308. IF(IERR.NE.0) GOTO 9999
  309.  
  310. c CALL ECROBJ('MAILLAGE',MAPCOM)
  311. c CALL NBNO
  312. c CALL LIRENT(NBNCOM,1,IRETOU)
  313. c IF(IERR.NE.0) GOTO 9999
  314.  
  315. CALL CHANGE(MAPCOM,1)
  316. NBNCOM=MAPCOM.NUM(/2)
  317.  
  318. C ... Maintenant on va enlever les noeuds du MAPCOM des deux autres
  319. C maillages ...
  320.  
  321. CALL ECROBJ('MAILLAGE',MAPBMA)
  322. CALL ECROBJ('MAILLAGE',MAPCOM)
  323. CALL PRDIFF
  324. CALL LIROBJ('MAILLAGE',MAPBMA,1,IRETOU)
  325. IF(IERR.NE.0) GOTO 9999
  326. NBNBMA=NBNBMA-NBNCOM
  327.  
  328. CALL ECROBJ('MAILLAGE',MAPCOQ)
  329. CALL ECROBJ('MAILLAGE',MAPCOM)
  330. CALL PRDIFF
  331. CALL LIROBJ('MAILLAGE',MAPCOQ,1,IRETOU)
  332. IF(IERR.NE.0) GOTO 9999
  333. NBNCOQ=NBNCOQ-NBNCOM
  334. ELSE
  335. NBNCOM=0
  336. MAPCOM=0
  337. ENDIF
  338.  
  339. ELSE
  340. NBNCOM=0
  341. MAPCOM=0
  342. ENDIF
  343.  
  344. IF (IIMPI.GE.4) THEN
  345. write(*,*) 'NBNBMA = ',NBNBMA
  346. write(*,*) 'NBNCOQ = ',NBNCOQ
  347. write(*,*) 'NBNCOM = ',NBNCOM
  348. ENDIF
  349.  
  350. C ... On créé le maillage - support de la rigidité ...
  351. C ... En même temps on va chercher les N° mini et maxi des noeuds ...
  352. NBNN=NBNBMA+NBNCOQ+NBNCOM
  353. NBELEM=1
  354. NBSOUS=0
  355. NBREF=NBSZEL
  356. SEGINI,MELEME
  357. ITYPEL=28
  358.  
  359. IF(NBNBMA.GT.0) THEN
  360. SEGACT,MAPBMA
  361. NOMINI=MAPBMA.NUM(1,1)
  362. NOMAXI=MAPBMA.NUM(1,1)
  363. DO 1100 I=1,NBNBMA
  364. NUM(I,1)=MAPBMA.NUM(1,I)
  365. IF(MAPBMA.NUM(1,I).LT.NOMINI) NOMINI=MAPBMA.NUM(1,I)
  366. IF(MAPBMA.NUM(1,I).GT.NOMAXI) NOMAXI=MAPBMA.NUM(1,I)
  367. 1100 CONTINUE
  368. SEGDES,MAPBMA
  369. ENDIF
  370.  
  371. IF(NBNCOQ.GT.0) THEN
  372. SEGACT,MAPCOQ
  373. IF(NBNBMA.EQ.0) THEN
  374. NOMINI=MAPCOQ.NUM(1,1)
  375. NOMAXI=MAPCOQ.NUM(1,1)
  376. ENDIF
  377. DO 1110 I=1,NBNCOQ
  378. NUM(I+NBNBMA,1)=MAPCOQ.NUM(1,I)
  379. IF(MAPCOQ.NUM(1,I).LT.NOMINI) NOMINI=MAPCOQ.NUM(1,I)
  380. IF(MAPCOQ.NUM(1,I).GT.NOMAXI) NOMAXI=MAPCOQ.NUM(1,I)
  381. 1110 CONTINUE
  382. SEGDES,MAPCOQ
  383. ENDIF
  384.  
  385. IF(NBNCOM.GT.0) THEN
  386. SEGACT,MAPCOM
  387. IF(NBNBMA+NBNCOQ.EQ.0) THEN
  388. NOMINI=MAPCOM.NUM(1,1)
  389. NOMAXI=MAPCOM.NUM(1,1)
  390. ENDIF
  391. DO 1120 I=1,NBNCOM
  392. NUM(I+NBNBMA+NBNCOQ,1)=MAPCOM.NUM(1,I)
  393. IF(MAPCOM.NUM(1,I).LT.NOMINI) NOMINI=MAPCOM.NUM(1,I)
  394. IF(MAPCOM.NUM(1,I).GT.NOMAXI) NOMAXI=MAPCOM.NUM(1,I)
  395. 1120 CONTINUE
  396. SEGDES,MAPCOM
  397. ENDIF
  398.  
  399. DO 1200 I=1,NBREF
  400. IMODEL=KMODEL(I)
  401. LISREF(I)=IMAMOD
  402. c-inutile SEGDES,IMODEL
  403. C ... On profite de cette boucle pour activer les maillages ...
  404. IPT1=IMAMOD
  405. SEGACT,IPT1
  406. 1200 CONTINUE
  407. ICOLOR(1)=7
  408.  
  409. cdebug segprt,meleme
  410. SEGDES,MELEME
  411.  
  412. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  413. C ... Préparation de la structure du MRIGID ... C
  414. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  415.  
  416. NRIGE=7
  417. NRIGEL=1
  418. SEGINI,MRIGID
  419. MTYMAT='RIGIDITE'
  420. COERIG(1)=1.D0
  421. IRIGEL(1,1)=MELEME
  422. C ... Les noeuds des massifs ont un seul DDL : T
  423. C coques ont deux DDL : TINF TSUP
  424. C communs ont trois DDL : T TINF et TSUP
  425. NLIGRP=NBNBMA+2*NBNCOQ+3*NBNCOM
  426. NLIGRD=NLIGRP
  427. SEGINI,DESCR
  428.  
  429. DO 2000 I=1,NBNBMA
  430. LISINC(I)='T '
  431. LISDUA(I)='Q '
  432. NOELEP(I)=I
  433. NOELED(I)=I
  434. 2000 CONTINUE
  435.  
  436. DO 2010 I=1,NBNCOQ
  437. LISINC(NBNBMA+2*I-1)='TSUP'
  438. LISINC(NBNBMA+2*I )='TINF'
  439. LISDUA(NBNBMA+2*I-1)='QSUP'
  440. LISDUA(NBNBMA+2*I )='QINF'
  441. NOELEP(NBNBMA+2*I-1)=NBNBMA+I
  442. NOELEP(NBNBMA+2*I )=NBNBMA+I
  443. NOELED(NBNBMA+2*I-1)=NBNBMA+I
  444. NOELED(NBNBMA+2*I )=NBNBMA+I
  445. 2010 CONTINUE
  446.  
  447. DO 2020 I=1,NBNCOM
  448. LISINC(2*NBNCOQ+NBNBMA+3*I-2)='T'
  449. LISINC(2*NBNCOQ+NBNBMA+3*I-1)='TSUP'
  450. LISINC(2*NBNCOQ+NBNBMA+3*I )='TINF'
  451. LISDUA(2*NBNCOQ+NBNBMA+3*I-2)='Q'
  452. LISDUA(2*NBNCOQ+NBNBMA+3*I-1)='QSUP'
  453. LISDUA(2*NBNCOQ+NBNBMA+3*I )='QINF'
  454. NOELEP(2*NBNCOQ+NBNBMA+3*I-2)=2*NBNCOQ+NBNBMA+I
  455. NOELEP(2*NBNCOQ+NBNBMA+3*I-1)=2*NBNCOQ+NBNBMA+I
  456. NOELEP(2*NBNCOQ+NBNBMA+3*I )=2*NBNCOQ+NBNBMA+I
  457. NOELED(2*NBNCOQ+NBNBMA+3*I-2)=2*NBNCOQ+NBNBMA+I
  458. NOELED(2*NBNCOQ+NBNBMA+3*I-1)=2*NBNCOQ+NBNBMA+I
  459. NOELED(2*NBNCOQ+NBNBMA+3*I )=2*NBNCOQ+NBNBMA+I
  460. 2020 CONTINUE
  461.  
  462. IRIGEL(3,1)=DESCR
  463.  
  464. cdebug segprt,descr
  465. SEGDES,DESCR
  466.  
  467. NELRIG=1
  468. * SEGINI,IMATRI
  469. SEGINI,XMATRI
  470. C ... On remplira XMATRI plus tard ...
  471. * IMATTT(1)=XMATRI
  472. * SEGDES,IMATRI
  473.  
  474. IRIGEL(4,1)=xMATRI
  475. IRIGEL(5,1)=NIFOUR
  476. IRIGEL(6,1)=0
  477. C ... On est dans le cas asymétrique ...
  478. IRIGEL(7,1)=2
  479. xmatri.symre=2
  480. IFORIG=IFOMOD
  481. SEGDES,MRIGID
  482.  
  483. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  484. C ... Calcul proprement dit de la matrice de rigidité ... C
  485. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  486.  
  487. C ... On commence par préparer le segment GL2LOC ...
  488. IECART=NOMAXI-NOMINI+1
  489. SEGINI,GL2LOC
  490.  
  491. C ... et le remplir ...
  492. IF(NBNBMA.GT.0) THEN
  493. SEGACT,MAPBMA
  494. DO 3000 I=1,NBNBMA
  495. IEQUIV(1,MAPBMA.NUM(1,I)-NOMINI+1)=I
  496. 3000 CONTINUE
  497. SEGDES,MAPBMA
  498. ENDIF
  499.  
  500. IF(NBNCOQ.GT.0) THEN
  501. SEGACT,MAPCOQ
  502. DO 3010 I=1,NBNCOQ
  503. IEQUIV(2,MAPCOQ.NUM(1,I)-NOMINI+1)=NBNBMA+2*I-1
  504. IEQUIV(3,MAPCOQ.NUM(1,I)-NOMINI+1)=NBNBMA+2*I
  505. 3010 CONTINUE
  506. SEGDES,MAPCOQ
  507. ENDIF
  508.  
  509. IF(NBNCOM.GT.0) THEN
  510. SEGACT,MAPCOM
  511. DO 3020 I=1,NBNCOQ
  512. IEQUIV(1,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I-2
  513. IEQUIV(2,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I-1
  514. IEQUIV(3,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I
  515. 3020 CONTINUE
  516. SEGDES,MAPCOM
  517. ENDIF
  518.  
  519. cdebug segprt,gl2loc
  520.  
  521. C ... On détruit les MAPOIN ? dont on n'aura plus besoin ? vraiment ?...
  522. *** IF(NBNBMA.GT.0) SEGSUP,MAPBMA
  523. IF(NBNCOQ.GT.0) SEGSUP,MAPCOQ
  524. IF(NBNCOM.GT.0) SEGSUP,MAPCOM
  525.  
  526. C ... On commence par extraire un segment INFO pour chaque zone
  527. C élémentaire. On profite de cette boucle pour trouver le
  528. C maximum des NBNN sur toutes les zones élémentaires ...
  529.  
  530. SEGINI,LESINF
  531. DO 3900 I=1,NBSZEL
  532.  
  533. IMODEL=KMODEL(I)
  534. MELE=NEFMOD
  535. segact imodel
  536. if( 2+inttyp.gt.infmod(/1) ) then
  537. CALL ELQUOI(MELE,0,INTTYP,ikk,IMODEL)
  538. info = ikk
  539. else
  540. jg=16
  541. segini info
  542. do iou=1,16
  543. infell(iou)=infele(iou)
  544. enddo
  545. infell(11)=infmod(2+inttyp)
  546. endif
  547. LINFO(I)=info
  548. MELEME=IMAMOD
  549. NBNN=NUM(/1)
  550. NBEL=NUM(/2)
  551. IF(I.EQ.1) THEN
  552. NBNNMX = NBNN
  553. NBELMX = NBEL
  554. ENDIF
  555. IF(NBNN.GT.NBNNMX) NBNNMX = NBNN
  556. IF(NBEL.GT.NBELMX) NBELMX = NBEL
  557.  
  558. 3900 CONTINUE
  559.  
  560. C ... NBNNMX et NBELMX servent à dimensionner le segment LESAIJ
  561. C qui contiendra les termes suivants :
  562. C
  563. C ik /
  564. C a = | N sur l'élément N° k de la zone élémentaire N° i,
  565. C j / j j étant l'indice du noeud de l'élément (et de
  566. C la fonction de forme associée)
  567. C
  568.  
  569. SEGINI,LESAIJ
  570. IDIM1=IDIM-1
  571. DO 3910 I=1,NBSZEL
  572.  
  573. IMODEL=KMODEL(I)
  574. MELEME=IMAMOD
  575. NBNN=NUM(/1)
  576. NBEL=NUM(/2)
  577.  
  578. INFO=LINFO(I)
  579. MINTE=INFELL(11)
  580. SEGACT,MINTE
  581. cdebug write(*,*) '-----------------------------'
  582. cdebug write(*,*) 'ZONE N° ',I
  583. cdebug segprt,minte
  584. NBNGAU=POIGAU(/1)
  585.  
  586. DO 3915 IEL=1,NBEL
  587.  
  588. CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XEL)
  589.  
  590. C ... SEG2 ...
  591. IF(ITYPEL.EQ.2) THEN
  592. XCH(1,1)=XEL(1,1)
  593. XCH(2,1)=XEL(2,1)
  594. XCH(1,2)=XEL(1,2)
  595. XCH(2,2)=XEL(2,2)
  596. C ... SEG3 ...
  597. ELSE IF(ITYPEL.EQ.3) THEN
  598. XCH(1,1)=XEL(1,1)
  599. XCH(2,1)=XEL(2,1)
  600. XCH(1,2)=XEL(1,3)
  601. XCH(2,2)=XEL(2,3)
  602. C ... TRI3 ou QUA4 ...
  603. ELSE IF(ITYPEL.EQ.4 .OR. ITYPEL.EQ.8) THEN
  604. XCH(1,1)=XEL(1,1)
  605. XCH(2,1)=XEL(2,1)
  606. XCH(3,1)=XEL(3,1)
  607. XCH(1,2)=XEL(1,2)
  608. XCH(2,2)=XEL(2,2)
  609. XCH(3,2)=XEL(3,2)
  610. XCH(1,3)=XEL(1,3)
  611. XCH(2,3)=XEL(2,3)
  612. XCH(3,3)=XEL(3,3)
  613. C ... TRI6 ou QUA8 ...
  614. ELSE IF(ITYPEL.EQ.6 .OR. ITYPEL.EQ.10) THEN
  615. XCH(1,1)=XEL(1,1)
  616. XCH(2,1)=XEL(2,1)
  617. XCH(3,1)=XEL(3,1)
  618. XCH(1,2)=XEL(1,3)
  619. XCH(2,2)=XEL(2,3)
  620. XCH(3,2)=XEL(3,3)
  621. XCH(1,3)=XEL(1,5)
  622. XCH(2,3)=XEL(2,5)
  623. XCH(3,3)=XEL(3,5)
  624. ELSE
  625. CALL ERREUR(16)
  626. RETURN
  627. ENDIF
  628.  
  629. IF(IDIM.EQ.2)THEN
  630. CALL VPAST2(XCH,BPSS)
  631. ELSE IF(IDIM.EQ.3) THEN
  632. CALL VPAST(XCH,BPSS)
  633. ENDIF
  634.  
  635. CALL VCORL1(XEL,XTR,BPSS,NBNN)
  636.  
  637. DO 3920 IPG=1,NBNGAU
  638.  
  639. DO 3930 INN=1,NBNN
  640. SHP(1,INN)=SHPTOT(1,INN,IPG)
  641. SHP(2,INN)=SHPTOT(2,INN,IPG)
  642. SHP(3,INN)=SHPTOT(3,INN,IPG)
  643. 3930 CONTINUE
  644.  
  645. CALL JACOBI(XTR,SHP,IDIM1,NBNN,DJAC)
  646. cdebug write(*,*) 'Zone ',I,', Element ',IEL,', PG ',IPG,
  647. cdebug & '-> J = ',DJAC
  648.  
  649. IF(IFOMOD.NE.0) THEN
  650.  
  651. DO 3935 INN=1,NBNN
  652. AIJ(INN,IEL,I) = AIJ(INN,IEL,I) +
  653. & SHPTOT(1,INN,IPG)*POIGAU(IPG)*DJAC
  654. 3935 CONTINUE
  655.  
  656. ELSE
  657. C ... axisymetrique : calcul du rayon au point d'integration
  658. C
  659. CALL DISTRR(XEL,SHP,NBNN,RR)
  660. RR=RR*2.D0*XPI
  661.  
  662. DO 3936 INN=1,NBNN
  663. AIJ(INN,IEL,I) = AIJ(INN,IEL,I) +
  664. & SHPTOT(1,INN,IPG)*POIGAU(IPG)*RR*DJAC
  665. 3936 CONTINUE
  666. ENDIF
  667.  
  668. 3920 CONTINUE
  669.  
  670. 3915 CONTINUE
  671.  
  672. 3910 CONTINUE
  673.  
  674. cdebug segprt,lesaij
  675.  
  676. C ... Double boucle sur tous les éléments concernés ...
  677. C ... Tout ce qui se termine par 1 concerne la boucle extérieure :
  678. C I1, IMODE1, IPT1, INFO1, MINTE1 ...
  679.  
  680. SEGACT,CHAMTM
  681. CALL TCHAMR(CHAMR,1)
  682.  
  683. BILAN = 0.D0
  684. DO 4000 I1=1,NBSZEL
  685.  
  686. IMODE1=KMODEL(I1)
  687. IPT1=IMODE1.IMAMOD
  688. NBNN1=IPT1.NUM(/1)
  689.  
  690. C ... On va chercher la température moyenne de l'élément IEL1 ...
  691. C ... On commence par trouver les bons numéros de composantes,
  692. C puis on active les MELVAL associés ...
  693.  
  694. SEGACT,CHAMTM
  695. MCHAM1=CHAMTM.ICHAML(I1)
  696. C write(*,*) ' I1 MCHAM1 = ',I1,MCHAM1
  697. SEGACT,MCHAM1
  698. NC=MCHAM1.IELVAL(/1)
  699. IF(KCOQ(I1)) THEN
  700. NCTINF=0
  701. NCTSUP=0
  702. DO 4005 I=1,NC
  703. IF(MCHAM1.NOMCHE(I).EQ.'TINF ') NCTINF=I
  704. IF(MCHAM1.NOMCHE(I).EQ.'TSUP ') NCTSUP=I
  705. 4005 CONTINUE
  706.  
  707. IF(NCTINF.EQ.0) GOTO 9999
  708. MELVA1=MCHAM1.IELVAL(NCTINF)
  709. SEGACT,MELVA1
  710.  
  711. IF(NCTSUP.EQ.0) GOTO 9999
  712. MELVA2=MCHAM1.IELVAL(NCTSUP)
  713. SEGACT,MELVA2
  714. ELSE
  715. NCT=0
  716. DO 4006 I=1,NC
  717. IF(MCHAM1.NOMCHE(I).EQ.'T ') NCT=I
  718. 4006 CONTINUE
  719.  
  720. IF(NCT.EQ.0) GOTO 9999
  721. MELVA1=MCHAM1.IELVAL(NCT)
  722. SEGACT,MELVA1
  723. ENDIF
  724. SEGDES,MCHAM1
  725. SEGDES,CHAMTM
  726.  
  727. DO 4010 IEL1=1,IPT1.NUM(/2)
  728.  
  729. C ... On trouve la (les) température(s) correspondants à
  730. C l'élément IEL1 ...
  731.  
  732. IF(KCOQ(I1)) THEN
  733. IF(MELVA1.VELCHE(/2).EQ.1) THEN
  734. T=MELVA1.VELCHE(1,1)
  735. ELSE
  736. T=MELVA1.VELCHE(1,IEL1)
  737. ENDIF
  738. T3(2)=T*T*T
  739.  
  740. IF(MELVA2.VELCHE(/2).EQ.1) THEN
  741. T=MELVA2.VELCHE(1,1)
  742. ELSE
  743. T=MELVA2.VELCHE(1,IEL1)
  744. ENDIF
  745. T3(1)=T*T*T
  746. ELSE
  747. IF(MELVA1.VELCHE(/2).EQ.1) THEN
  748. T=MELVA1.VELCHE(1,1)
  749. ELSE
  750. T=MELVA1.VELCHE(1,IEL1)
  751. ENDIF
  752. T3(1)=T*T*T
  753. ENDIF
  754.  
  755. C ... Deuxième niveau de la boucle sur les éléments ...
  756.  
  757. DO 4020 I2=1,NBSZEL
  758.  
  759. IMODE2=KMODEL(I2)
  760. IPT2=IMODE2.IMAMOD
  761. NBNN2=IPT2.NUM(/1)
  762.  
  763. DO 4030 IEL2=1,IPT2.NUM(/2)
  764.  
  765. C ... On va chercher les surfaces des éléments IEL1 et IEL2 ...
  766.  
  767. MCHAML=CHAMR.ICHAML(I2)
  768. SEGACT,MCHAML
  769. NC2=IELVAL(/1)
  770. C ... Ici on suppose que la surface est toujours dernière ...
  771. MELVAL=IELVAL(NC2)
  772. SURF2=VELCHE(1,IEL2)
  773.  
  774. MCHAML=CHAMR.ICHAML(I1)
  775. SEGACT,MCHAML
  776. NC1=IELVAL(/1)
  777. C ... Ici on suppose que la surface est toujours dernière ...
  778. MELVAL=IELVAL(NC1)
  779. SURF1=VELCHE(1,IEL1)
  780.  
  781. C ... On a supposé que dans le champ de facteurs de forme les
  782. C composantes sont toujours dans l'ordre suivant :
  783. C MIDL SURF dans le cas des solides
  784. C SUPE INFE SURF dans le cas des coques
  785. C
  786. DO 4040 IC1=1,NC1-1
  787. MELVAL=IELVAL(IC1)
  788. MCHELM=IELCHE(1,IEL1)
  789. MCHAML=ICHAML(I2)
  790. DO 4050 IC2=1,NC2-1
  791. MELVAL=IELVAL(IC2)
  792. COEFF=VELCHE(1,IEL2)
  793. C write(*,*)' iel1 iel2 coeff ',IEL1,IEL2,COEFF
  794.  
  795. DO 4060 INN1=1,NBNN1
  796. NGL1=IPT1.NUM(INN1,IEL1)
  797. NDDL1=IC1
  798. IF(KCOQ(I1)) NDDL1=NDDL1+1
  799. NLOC1 = IEQUIV(NDDL1,NGL1-NOMINI+1)
  800. DO 4070 INN2=1,NBNN2
  801.  
  802. NGL2=IPT2.NUM(INN2,IEL2)
  803. NDDL2=IC2
  804. IF(KCOQ(I2)) NDDL2=NDDL2+1
  805. NLOC2 = IEQUIV(NDDL2,NGL2-NOMINI+1)
  806.  
  807. C write(*,*)' in1 in2 coeff ',INN1,INN2,COEFF
  808. TERME = COEFF*AIJ(INN1,IEL1,I1)
  809. & *AIJ(INN2,IEL2,I2)*T3(IC1)*CONSSB/(SURF1*SURF2)
  810. RE(NLOC2,NLOC1,1)=RE(NLOC2,NLOC1,1) + TERME
  811. BILAN = BILAN + TERME
  812.  
  813. cdebug write(*,*) 'On rajoute ',TERME,' au terme ',NLOC1,NLOC2
  814.  
  815. 4070 CONTINUE
  816.  
  817. 4060 CONTINUE
  818.  
  819. 4050 CONTINUE
  820.  
  821. 4040 CONTINUE
  822.  
  823. 4030 CONTINUE
  824.  
  825. 4020 CONTINUE
  826.  
  827. 4010 CONTINUE
  828. SEGDES,MELVA1
  829. IF(KCOQ(I1)) SEGDES,MELVA2
  830.  
  831. 4000 CONTINUE
  832. IF (IIMPI.GE.4) THEN
  833. write(*,*) 'Bilan = ',BILAN
  834. ENDIF
  835.  
  836. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  837. C ... Sortie de la rigidité calculée ... C
  838. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  839.  
  840. C!! segprt,xmatri
  841. SEGDES,XMATRI
  842. CALL ECROBJ('RIGIDITE',MRIGID)
  843.  
  844. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  845. C ... Le ménage à la fin ... C
  846. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  847.  
  848. C ... On désactive le modèle ...
  849. DO 5000 I=1,NBSZEL
  850. IMODEL=KMODEL(I)
  851. MELEME=IMAMOD
  852. SEGDES,MELEME
  853. SEGDES,IMODEL
  854. INFO=LINFO(I)
  855. SEGSUP,INFO
  856. 5000 CONTINUE
  857. SEGDES,MMODEL
  858. SEGSUP,LESINF
  859.  
  860. C ... desactivation de la matrice de rayonnement
  861. CALL TCHAMR(CHAMR,0)
  862.  
  863. C ... On supprime les segments de travail ...
  864. SEGSUP,INFOEL
  865. SEGSUP,GL2LOC
  866. SEGSUP,LESAIJ
  867.  
  868. RETURN
  869.  
  870. 9999 CONTINUE
  871. CALL ERREUR(26)
  872. RETURN
  873. END
  874.  
  875.  
  876.  
  877.  
  878.  
  879.  
  880.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales