Télécharger rayn.eso

Retour à la liste

Numérotation des lignes :

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

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