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

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