Télécharger rayn.eso

Retour à la liste

Numérotation des lignes :

  1. C RAYN SOURCE PV 13/04/12 21:15:51 7756
  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. IFORIG=IFOMOD
  480. SEGDES,MRIGID
  481.  
  482. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  483. C ... Calcul proprement dit de la matrice de rigidité ... C
  484. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  485.  
  486. C ... On commence par préparer le segment GL2LOC ...
  487. IECART=NOMAXI-NOMINI+1
  488. SEGINI,GL2LOC
  489.  
  490. C ... et le remplir ...
  491. IF(NBNBMA.GT.0) THEN
  492. SEGACT,MAPBMA
  493. DO 3000 I=1,NBNBMA
  494. IEQUIV(1,MAPBMA.NUM(1,I)-NOMINI+1)=I
  495. 3000 CONTINUE
  496. SEGDES,MAPBMA
  497. ENDIF
  498.  
  499. IF(NBNCOQ.GT.0) THEN
  500. SEGACT,MAPCOQ
  501. DO 3010 I=1,NBNCOQ
  502. IEQUIV(2,MAPCOQ.NUM(1,I)-NOMINI+1)=NBNBMA+2*I-1
  503. IEQUIV(3,MAPCOQ.NUM(1,I)-NOMINI+1)=NBNBMA+2*I
  504. 3010 CONTINUE
  505. SEGDES,MAPCOQ
  506. ENDIF
  507.  
  508. IF(NBNCOM.GT.0) THEN
  509. SEGACT,MAPCOM
  510. DO 3020 I=1,NBNCOQ
  511. IEQUIV(1,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I-2
  512. IEQUIV(2,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I-1
  513. IEQUIV(3,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I
  514. 3020 CONTINUE
  515. SEGDES,MAPCOM
  516. ENDIF
  517.  
  518. cdebug segprt,gl2loc
  519.  
  520. C ... On détruit les MAPOIN ? dont on n'aura plus besoin ? vraiment ?...
  521. *** IF(NBNBMA.GT.0) SEGSUP,MAPBMA
  522. IF(NBNCOQ.GT.0) SEGSUP,MAPCOQ
  523. IF(NBNCOM.GT.0) SEGSUP,MAPCOM
  524.  
  525. C ... On commence par extraire un segment INFO pour chaque zone
  526. C élémentaire. On profite de cette boucle pour trouver le
  527. C maximum des NBNN sur toutes les zones élémentaires ...
  528.  
  529. SEGINI,LESINF
  530. DO 3900 I=1,NBSZEL
  531.  
  532. IMODEL=KMODEL(I)
  533. MELE=NEFMOD
  534. segact imodel
  535. if( 2+inttyp.gt.infmod(/1) ) then
  536. CALL ELQUOI(MELE,0,INTTYP,ikk,IMODEL)
  537. info = ikk
  538. else
  539. jg=16
  540. segini info
  541. do iou=1,16
  542. infell(iou)=infele(iou)
  543. enddo
  544. infell(11)=infmod(2+inttyp)
  545. endif
  546. LINFO(I)=info
  547. MELEME=IMAMOD
  548. NBNN=NUM(/1)
  549. NBEL=NUM(/2)
  550. IF(I.EQ.1) THEN
  551. NBNNMX = NBNN
  552. NBELMX = NBEL
  553. ENDIF
  554. IF(NBNN.GT.NBNNMX) NBNNMX = NBNN
  555. IF(NBEL.GT.NBELMX) NBELMX = NBEL
  556.  
  557. 3900 CONTINUE
  558.  
  559. C ... NBNNMX et NBELMX servent à dimensionner le segment LESAIJ
  560. C qui contiendra les termes suivants :
  561. C
  562. C ik /
  563. C a = | N sur l'élément N° k de la zone élémentaire N° i,
  564. C j / j j étant l'indice du noeud de l'élément (et de
  565. C la fonction de forme associée)
  566. C
  567.  
  568. SEGINI,LESAIJ
  569. IDIM1=IDIM-1
  570. DO 3910 I=1,NBSZEL
  571.  
  572. IMODEL=KMODEL(I)
  573. MELEME=IMAMOD
  574. NBNN=NUM(/1)
  575. NBEL=NUM(/2)
  576.  
  577. INFO=LINFO(I)
  578. MINTE=INFELL(11)
  579. SEGACT,MINTE
  580. cdebug write(*,*) '-----------------------------'
  581. cdebug write(*,*) 'ZONE N° ',I
  582. cdebug segprt,minte
  583. NBNGAU=POIGAU(/1)
  584.  
  585. DO 3915 IEL=1,NBEL
  586.  
  587. CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XEL)
  588.  
  589. C ... SEG2 ...
  590. IF(ITYPEL.EQ.2) THEN
  591. XCH(1,1)=XEL(1,1)
  592. XCH(2,1)=XEL(2,1)
  593. XCH(1,2)=XEL(1,2)
  594. XCH(2,2)=XEL(2,2)
  595. C ... SEG3 ...
  596. ELSE IF(ITYPEL.EQ.3) THEN
  597. XCH(1,1)=XEL(1,1)
  598. XCH(2,1)=XEL(2,1)
  599. XCH(1,2)=XEL(1,3)
  600. XCH(2,2)=XEL(2,3)
  601. C ... TRI3 ou QUA4 ...
  602. ELSE IF(ITYPEL.EQ.4 .OR. ITYPEL.EQ.8) THEN
  603. XCH(1,1)=XEL(1,1)
  604. XCH(2,1)=XEL(2,1)
  605. XCH(3,1)=XEL(3,1)
  606. XCH(1,2)=XEL(1,2)
  607. XCH(2,2)=XEL(2,2)
  608. XCH(3,2)=XEL(3,2)
  609. XCH(1,3)=XEL(1,3)
  610. XCH(2,3)=XEL(2,3)
  611. XCH(3,3)=XEL(3,3)
  612. C ... TRI6 ou QUA8 ...
  613. ELSE IF(ITYPEL.EQ.6 .OR. ITYPEL.EQ.10) THEN
  614. XCH(1,1)=XEL(1,1)
  615. XCH(2,1)=XEL(2,1)
  616. XCH(3,1)=XEL(3,1)
  617. XCH(1,2)=XEL(1,3)
  618. XCH(2,2)=XEL(2,3)
  619. XCH(3,2)=XEL(3,3)
  620. XCH(1,3)=XEL(1,5)
  621. XCH(2,3)=XEL(2,5)
  622. XCH(3,3)=XEL(3,5)
  623. ELSE
  624. CALL ERREUR(16)
  625. RETURN
  626. ENDIF
  627.  
  628. IF(IDIM.EQ.2)THEN
  629. CALL VPAST2(XCH,BPSS)
  630. ELSE IF(IDIM.EQ.3) THEN
  631. CALL VPAST(XCH,BPSS)
  632. ENDIF
  633.  
  634. CALL VCORL1(XEL,XTR,BPSS,NBNN)
  635.  
  636. DO 3920 IPG=1,NBNGAU
  637.  
  638. DO 3930 INN=1,NBNN
  639. SHP(1,INN)=SHPTOT(1,INN,IPG)
  640. SHP(2,INN)=SHPTOT(2,INN,IPG)
  641. SHP(3,INN)=SHPTOT(3,INN,IPG)
  642. 3930 CONTINUE
  643.  
  644. CALL JACOBI(XTR,SHP,IDIM1,NBNN,DJAC)
  645. cdebug write(*,*) 'Zone ',I,', Element ',IEL,', PG ',IPG,
  646. cdebug & '-> J = ',DJAC
  647.  
  648. IF(IFOMOD.NE.0) THEN
  649.  
  650. DO 3935 INN=1,NBNN
  651. AIJ(INN,IEL,I) = AIJ(INN,IEL,I) +
  652. & SHPTOT(1,INN,IPG)*POIGAU(IPG)*DJAC
  653. 3935 CONTINUE
  654.  
  655. ELSE
  656. C ... axisymetrique : calcul du rayon au point d'integration
  657. C
  658. CALL DISTRR(XEL,SHP,NBNN,RR)
  659. RR=RR*2.D0*XPI
  660.  
  661. DO 3936 INN=1,NBNN
  662. AIJ(INN,IEL,I) = AIJ(INN,IEL,I) +
  663. & SHPTOT(1,INN,IPG)*POIGAU(IPG)*RR*DJAC
  664. 3936 CONTINUE
  665. ENDIF
  666.  
  667. 3920 CONTINUE
  668.  
  669. 3915 CONTINUE
  670.  
  671. 3910 CONTINUE
  672.  
  673. cdebug segprt,lesaij
  674.  
  675. C ... Double boucle sur tous les éléments concernés ...
  676. C ... Tout ce qui se termine par 1 concerne la boucle extérieure :
  677. C I1, IMODE1, IPT1, INFO1, MINTE1 ...
  678.  
  679. SEGACT,CHAMTM
  680. CALL TCHAMR(CHAMR,1)
  681.  
  682. BILAN = 0.D0
  683. DO 4000 I1=1,NBSZEL
  684.  
  685. IMODE1=KMODEL(I1)
  686. IPT1=IMODE1.IMAMOD
  687. NBNN1=IPT1.NUM(/1)
  688.  
  689. C ... On va chercher la température moyenne de l'élément IEL1 ...
  690. C ... On commence par trouver les bons numéros de composantes,
  691. C puis on active les MELVAL associés ...
  692.  
  693. SEGACT,CHAMTM
  694. MCHAM1=CHAMTM.ICHAML(I1)
  695. C write(*,*) ' I1 MCHAM1 = ',I1,MCHAM1
  696. SEGACT,MCHAM1
  697. NC=MCHAM1.IELVAL(/1)
  698. IF(KCOQ(I1)) THEN
  699. NCTINF=0
  700. NCTSUP=0
  701. DO 4005 I=1,NC
  702. IF(MCHAM1.NOMCHE(I).EQ.'TINF ') NCTINF=I
  703. IF(MCHAM1.NOMCHE(I).EQ.'TSUP ') NCTSUP=I
  704. 4005 CONTINUE
  705.  
  706. IF(NCTINF.EQ.0) GOTO 9999
  707. MELVA1=MCHAM1.IELVAL(NCTINF)
  708. SEGACT,MELVA1
  709.  
  710. IF(NCTSUP.EQ.0) GOTO 9999
  711. MELVA2=MCHAM1.IELVAL(NCTSUP)
  712. SEGACT,MELVA2
  713. ELSE
  714. NCT=0
  715. DO 4006 I=1,NC
  716. IF(MCHAM1.NOMCHE(I).EQ.'T ') NCT=I
  717. 4006 CONTINUE
  718.  
  719. IF(NCT.EQ.0) GOTO 9999
  720. MELVA1=MCHAM1.IELVAL(NCT)
  721. SEGACT,MELVA1
  722. ENDIF
  723. SEGDES,MCHAM1
  724. SEGDES,CHAMTM
  725.  
  726. DO 4010 IEL1=1,IPT1.NUM(/2)
  727.  
  728. C ... On trouve la (les) température(s) correspondants à
  729. C l'élément IEL1 ...
  730.  
  731. IF(KCOQ(I1)) THEN
  732. IF(MELVA1.VELCHE(/2).EQ.1) THEN
  733. T=MELVA1.VELCHE(1,1)
  734. ELSE
  735. T=MELVA1.VELCHE(1,IEL1)
  736. ENDIF
  737. T3(2)=T*T*T
  738.  
  739. IF(MELVA2.VELCHE(/2).EQ.1) THEN
  740. T=MELVA2.VELCHE(1,1)
  741. ELSE
  742. T=MELVA2.VELCHE(1,IEL1)
  743. ENDIF
  744. T3(1)=T*T*T
  745. ELSE
  746. IF(MELVA1.VELCHE(/2).EQ.1) THEN
  747. T=MELVA1.VELCHE(1,1)
  748. ELSE
  749. T=MELVA1.VELCHE(1,IEL1)
  750. ENDIF
  751. T3(1)=T*T*T
  752. ENDIF
  753.  
  754. C ... Deuxième niveau de la boucle sur les éléments ...
  755.  
  756. DO 4020 I2=1,NBSZEL
  757.  
  758. IMODE2=KMODEL(I2)
  759. IPT2=IMODE2.IMAMOD
  760. NBNN2=IPT2.NUM(/1)
  761.  
  762. DO 4030 IEL2=1,IPT2.NUM(/2)
  763.  
  764. C ... On va chercher les surfaces des éléments IEL1 et IEL2 ...
  765.  
  766. MCHAML=CHAMR.ICHAML(I2)
  767. SEGACT,MCHAML
  768. NC2=IELVAL(/1)
  769. C ... Ici on suppose que la surface est toujours dernière ...
  770. MELVAL=IELVAL(NC2)
  771. SURF2=VELCHE(1,IEL2)
  772.  
  773. MCHAML=CHAMR.ICHAML(I1)
  774. SEGACT,MCHAML
  775. NC1=IELVAL(/1)
  776. C ... Ici on suppose que la surface est toujours dernière ...
  777. MELVAL=IELVAL(NC1)
  778. SURF1=VELCHE(1,IEL1)
  779.  
  780. C ... On a supposé que dans le champ de facteurs de forme les
  781. C composantes sont toujours dans l'ordre suivant :
  782. C MIDL SURF dans le cas des solides
  783. C SUPE INFE SURF dans le cas des coques
  784. C
  785. DO 4040 IC1=1,NC1-1
  786. MELVAL=IELVAL(IC1)
  787. MCHELM=IELCHE(1,IEL1)
  788. MCHAML=ICHAML(I2)
  789. DO 4050 IC2=1,NC2-1
  790. MELVAL=IELVAL(IC2)
  791. COEFF=VELCHE(1,IEL2)
  792. C write(*,*)' iel1 iel2 coeff ',IEL1,IEL2,COEFF
  793.  
  794. DO 4060 INN1=1,NBNN1
  795. NGL1=IPT1.NUM(INN1,IEL1)
  796. NDDL1=IC1
  797. IF(KCOQ(I1)) NDDL1=NDDL1+1
  798. NLOC1 = IEQUIV(NDDL1,NGL1-NOMINI+1)
  799. DO 4070 INN2=1,NBNN2
  800.  
  801. NGL2=IPT2.NUM(INN2,IEL2)
  802. NDDL2=IC2
  803. IF(KCOQ(I2)) NDDL2=NDDL2+1
  804. NLOC2 = IEQUIV(NDDL2,NGL2-NOMINI+1)
  805.  
  806. C write(*,*)' in1 in2 coeff ',INN1,INN2,COEFF
  807. TERME = COEFF*AIJ(INN1,IEL1,I1)
  808. & *AIJ(INN2,IEL2,I2)*T3(IC1)*CONSSB/(SURF1*SURF2)
  809. RE(NLOC2,NLOC1,1)=RE(NLOC2,NLOC1,1) + TERME
  810. BILAN = BILAN + TERME
  811.  
  812. cdebug write(*,*) 'On rajoute ',TERME,' au terme ',NLOC1,NLOC2
  813.  
  814. 4070 CONTINUE
  815.  
  816. 4060 CONTINUE
  817.  
  818. 4050 CONTINUE
  819.  
  820. 4040 CONTINUE
  821.  
  822. 4030 CONTINUE
  823.  
  824. 4020 CONTINUE
  825.  
  826. 4010 CONTINUE
  827. SEGDES,MELVA1
  828. IF(KCOQ(I1)) SEGDES,MELVA2
  829.  
  830. 4000 CONTINUE
  831. IF (IIMPI.GE.4) THEN
  832. write(*,*) 'Bilan = ',BILAN
  833. ENDIF
  834.  
  835. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  836. C ... Sortie de la rigidité calculée ... C
  837. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  838.  
  839. C!! segprt,xmatri
  840. SEGDES,XMATRI
  841. CALL ECROBJ('RIGIDITE',MRIGID)
  842.  
  843. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  844. C ... Le ménage à la fin ... C
  845. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  846.  
  847. C ... On désactive le modèle ...
  848. DO 5000 I=1,NBSZEL
  849. IMODEL=KMODEL(I)
  850. MELEME=IMAMOD
  851. SEGDES,MELEME
  852. SEGDES,IMODEL
  853. INFO=LINFO(I)
  854. SEGSUP,INFO
  855. 5000 CONTINUE
  856. SEGDES,MMODEL
  857. SEGSUP,LESINF
  858.  
  859. C ... desactivation de la matrice de rayonnement
  860. CALL TCHAMR(CHAMR,0)
  861.  
  862. C ... On supprime les segments de travail ...
  863. SEGSUP,INFOEL
  864. SEGSUP,GL2LOC
  865. SEGSUP,LESAIJ
  866.  
  867. RETURN
  868.  
  869. 9999 CONTINUE
  870. CALL ERREUR(26)
  871. RETURN
  872. END
  873.  
  874.  
  875.  
  876.  
  877.  
  878.  

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