Télécharger la plaquette Cast3M


fr en

Télécharger rayn.eso

Retour à la liste

Numérotation des lignes : 

rayn

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

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