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    CB215821  24/04/12    21:17:04     11897          
  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/INFO/(INFELL(JG)),INFO1.INFO,INFO2.INFO
  45.       SEGMENT LESINF
  46.          INTEGER LINFO(NBSZEL)
  47.       ENDSEGMENT
  48.  
  49.       SEGMENT LESAIJ
  50.          REAL*8 AIJ(NBNNMX,NBELMX,NBSZEL)
  51.       ENDSEGMENT
  52.  
  53.       DIMENSION T3(2)
  54. C ... Attention !!! 20 correspond ici au nombre maxi de noeuds par
  55. C     élément
  56. C ... XEL - coordonnées d'origine, XCH - coordonnées choisies pour
  57. C     VPAST (ou VPAST2),
  58. C     XTR - coordonnées transformées par VCORL1 ...
  59.       DIMENSION XEL(3,20),XCH(3,3),XTR(3,20)
  60. C ... BPSS - matrice de transformation ...
  61.       DIMENSION BPSS(3,3)
  62.       DIMENSION SHP(6,20)
  63.  
  64. C-------- Fin des déclarations  -------------------------------
  65.  
  66.       CHARACTER*8 TYPE
  67.       DATA INTTYP /3/
  68.       DATA CONSSB /5.67D-8/
  69.  
  70. C----------- Fin des DATA  ------------------------------------
  71.  
  72.       segact mcoord
  73. C ... Lecture ...
  74.  
  75.       TYPE='MMODEL  '
  76.       CALL LIROBJ(TYPE,IRET,1,IRETOU)
  77.       IF(IRETOU.EQ.0) THEN
  78.          MOTERR(1:8)=TYPE
  79.         CALL ERREUR(37)
  80.         RETURN
  81.       ENDIF
  82.       IF(IRETOU.EQ.1) MMODEL=IRET
  83.       SEGACT MMODEL
  84.  
  85.       TYPE='MCHAML  '
  86.       CALL LIROBJ(TYPE,IRET,1,IRETOU)
  87.       IF(IRETOU.EQ.0) THEN
  88.          MOTERR(1:8)=TYPE
  89.         CALL ERREUR(37)
  90.         RETURN
  91.       ENDIF
  92.       IF(IRETOU.EQ.1) MCHEL1=IRET
  93.  
  94.       TYPE='MCHAML  '
  95.       CALL LIROBJ(TYPE,IRET,1,IRETOU)
  96.       IF(IRETOU.EQ.0) THEN
  97.          MOTERR(1:8)=TYPE
  98.         CALL ERREUR(37)
  99.         RETURN
  100.       ENDIF
  101.       IF(IRETOU.EQ.1) MCHEL2=IRET
  102.  
  103.       SEGACT, MCHEL1, MCHEL2
  104.       IF((MCHEL1.TITCHE).EQ.'MATRICE DE RAYONNEMENT') THEN
  105.         CHAMR = MCHEL1
  106.         CHAMTM= MCHEL2
  107.       ELSEIF((MCHEL2.TITCHE).EQ.'MATRICE DE RAYONNEMENT') THEN
  108.         CHAMR = MCHEL2
  109.         CHAMTM= MCHEL1
  110.       ELSE
  111.         CALL ERREUR(21)
  112.         RETURN
  113.       ENDIF
  114.  
  115.       CALL LIRREE(XVAL,0,IRETOU)
  116.       IF(IRETOU.EQ.1) CONSSB=XVAL
  117.  
  118. C ... Verification  que les supports du MMODEL et des
  119. C     MCHAML sont bien les mêmes ...
  120.  
  121.       CALL RAYN1(MMODEL,CHAMR)
  122.       CALL RAYN1(MMODEL,CHAMTM)
  123.  
  124. C     CHAMTM = pointeur vers le MCHAML de température moyenne
  125. C     CHAMR  = pointeur vers le MCHAML de matrices de rayonnement
  126. C     MMODEL = pointeur vers le modèle
  127.  
  128. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  129. C ... Recherche du support géométrique de la rigidité ... C
  130. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  131.  
  132.       SEGACT,MMODEL
  133.  
  134. C ... On créé un MELEME contenant tous les maillages-supports
  135. C     du CHAMTM ...
  136.       NBSZEL=KMODEL(/1)
  137.       N1=NBSZEL
  138.       SEGINI,INFOEL
  139.       DO 1000 I=1,NBSZEL
  140.          IMODEL=KMODEL(I)
  141.          SEGACT,IMODEL
  142.          IF(FORMOD(1).NE.'THERMIQUE'.or.matmod(2).NE.'RAYONNEMENT')THEN
  143.             CALL ERREUR(21)
  144.             RETURN
  145.          ENDIF
  146.          MELEME=IMAMOD
  147.          SEGACT,MELEME
  148.  
  149.           IF (IDIM.EQ.3) THEN
  150. C
  151.                 IF ((NEFMOD.EQ.4).OR.(NEFMOD.EQ.8)) THEN
  152. C                  TRI3 ou QUA4
  153.                    KQUAD(I)=.FALSE.
  154.                    KCOQ(I) =.FALSE.
  155.                 ELSEIF ((NEFMOD.EQ.27).OR.(NEFMOD.EQ.49)) THEN
  156. C                  COQ3 ou COQ4
  157.                    KQUAD(I)=.FALSE.
  158.                    KCOQ(I) =.TRUE.
  159.                 ELSEIF ((NEFMOD.EQ.6).OR.(NEFMOD.EQ.10)) THEN
  160. C                  TRI6 ou QUA8
  161.                    KQUAD(I)=.TRUE.
  162.                    KCOQ(I) =.FALSE.
  163.                 ELSEIF ((NEFMOD.EQ.56).OR.(NEFMOD.EQ.41)) THEN
  164. C                  COQ6 ou COQ8
  165.                    KQUAD(I)=.TRUE.
  166.                    KCOQ(I) =.TRUE.
  167.                 ELSE
  168.                    CALL ERREUR(16)
  169.                    RETURN
  170.                 ENDIF
  171. C
  172.           ELSEIF (IDIM.EQ.2) THEN
  173. C
  174.                 IF (NEFMOD.EQ.2) THEN
  175. C                  SEG2
  176.                    KQUAD(I)=.FALSE.
  177.                    KCOQ(I) =.FALSE.
  178.                 ELSEIF (NEFMOD.EQ.44) THEN
  179. C                  COQ2
  180.                    KQUAD(I)=.FALSE.
  181.                    KCOQ(I) =.TRUE.
  182.                 ELSEIF (NEFMOD.EQ.3)  THEN
  183. C                  SEG3
  184.                    KQUAD(I)=.TRUE.
  185.                    KCOQ(I) =.FALSE.
  186.                 ELSE
  187.                    CALL ERREUR(16)
  188.                    RETURN
  189.                 ENDIF
  190. C
  191.           ELSE
  192.                 CALL ERREUR(14)
  193.                 RETURN
  194.           ENDIF
  195.  
  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.       ENDIF
  369.  
  370.       IF(NBNCOQ.GT.0) THEN
  371.          SEGACT,MAPCOQ
  372.          IF(NBNBMA.EQ.0) THEN
  373.             NOMINI=MAPCOQ.NUM(1,1)
  374.             NOMAXI=MAPCOQ.NUM(1,1)
  375.          ENDIF
  376.          DO 1110 I=1,NBNCOQ
  377.             NUM(I+NBNBMA,1)=MAPCOQ.NUM(1,I)
  378.             IF(MAPCOQ.NUM(1,I).LT.NOMINI) NOMINI=MAPCOQ.NUM(1,I)
  379.             IF(MAPCOQ.NUM(1,I).GT.NOMAXI) NOMAXI=MAPCOQ.NUM(1,I)
  380.  1110    CONTINUE
  381.       ENDIF
  382.  
  383.       IF(NBNCOM.GT.0) THEN
  384.          SEGACT,MAPCOM
  385.          IF(NBNBMA+NBNCOQ.EQ.0) THEN
  386.             NOMINI=MAPCOM.NUM(1,1)
  387.             NOMAXI=MAPCOM.NUM(1,1)
  388.          ENDIF
  389.          DO 1120 I=1,NBNCOM
  390.             NUM(I+NBNBMA+NBNCOQ,1)=MAPCOM.NUM(1,I)
  391.             IF(MAPCOM.NUM(1,I).LT.NOMINI) NOMINI=MAPCOM.NUM(1,I)
  392.             IF(MAPCOM.NUM(1,I).GT.NOMAXI) NOMAXI=MAPCOM.NUM(1,I)
  393.  1120    CONTINUE
  394.       ENDIF
  395.  
  396.       DO 1200 I=1,NBREF
  397.          IMODEL=KMODEL(I)
  398.          LISREF(I)=IMAMOD
  399. C    ... On profite de cette boucle pour activer les maillages ...
  400.          IPT1=IMAMOD
  401.          SEGACT,IPT1
  402.  1200 CONTINUE
  403.       ICOLOR(1)=7
  404.  
  405.  
  406. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  407. C ... Préparation de la structure du MRIGID           ... C
  408. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  409.  
  410.       NRIGE=7
  411.       NRIGEL=1
  412.       SEGINI,MRIGID
  413.       MTYMAT='RIGIDITE'
  414.       COERIG(1)=1.D0
  415.       IRIGEL(1,1)=MELEME
  416. C ... Les noeuds des massifs ont un seul DDL : T
  417. C                    coques ont deux DDL : TINF TSUP
  418. C                    communs ont trois DDL : T TINF et TSUP
  419.       NLIGRP=NBNBMA+2*NBNCOQ+3*NBNCOM
  420.       NLIGRD=NLIGRP
  421.       SEGINI,DESCR
  422.  
  423.       DO 2000 I=1,NBNBMA
  424.          LISINC(I)='T   '
  425.          LISDUA(I)='Q   '
  426.          NOELEP(I)=I
  427.          NOELED(I)=I
  428.  2000 CONTINUE
  429.  
  430.       DO 2010 I=1,NBNCOQ
  431.          LISINC(NBNBMA+2*I-1)='TSUP'
  432.          LISINC(NBNBMA+2*I  )='TINF'
  433.          LISDUA(NBNBMA+2*I-1)='QSUP'
  434.          LISDUA(NBNBMA+2*I  )='QINF'
  435.          NOELEP(NBNBMA+2*I-1)=NBNBMA+I
  436.          NOELEP(NBNBMA+2*I  )=NBNBMA+I
  437.          NOELED(NBNBMA+2*I-1)=NBNBMA+I
  438.          NOELED(NBNBMA+2*I  )=NBNBMA+I
  439.  2010 CONTINUE
  440.  
  441.       DO 2020 I=1,NBNCOM
  442.          LISINC(2*NBNCOQ+NBNBMA+3*I-2)='T'
  443.          LISINC(2*NBNCOQ+NBNBMA+3*I-1)='TSUP'
  444.          LISINC(2*NBNCOQ+NBNBMA+3*I  )='TINF'
  445.          LISDUA(2*NBNCOQ+NBNBMA+3*I-2)='Q'
  446.          LISDUA(2*NBNCOQ+NBNBMA+3*I-1)='QSUP'
  447.          LISDUA(2*NBNCOQ+NBNBMA+3*I  )='QINF'
  448.          NOELEP(2*NBNCOQ+NBNBMA+3*I-2)=2*NBNCOQ+NBNBMA+I
  449.          NOELEP(2*NBNCOQ+NBNBMA+3*I-1)=2*NBNCOQ+NBNBMA+I
  450.          NOELEP(2*NBNCOQ+NBNBMA+3*I  )=2*NBNCOQ+NBNBMA+I
  451.          NOELED(2*NBNCOQ+NBNBMA+3*I-2)=2*NBNCOQ+NBNBMA+I
  452.          NOELED(2*NBNCOQ+NBNBMA+3*I-1)=2*NBNCOQ+NBNBMA+I
  453.          NOELED(2*NBNCOQ+NBNBMA+3*I  )=2*NBNCOQ+NBNBMA+I
  454.  2020 CONTINUE
  455.  
  456.       IRIGEL(3,1)=DESCR
  457.  
  458. cdebug      segprt,descr
  459.       SEGDES,DESCR
  460.  
  461.       NELRIG=1
  462. *      SEGINI,IMATRI
  463.       SEGINI,XMATRI
  464. C ... On remplira XMATRI plus tard ...
  465. *      IMATTT(1)=XMATRI
  466. *      SEGDES,IMATRI
  467.  
  468.       IRIGEL(4,1)=xMATRI
  469.       IRIGEL(5,1)=NIFOUR
  470.       IRIGEL(6,1)=0
  471. C ... On est dans le cas asymétrique ...
  472.       IRIGEL(7,1)=2
  473.       xmatri.symre=2
  474.       IFORIG=IFOUR
  475.       SEGDES,MRIGID
  476.  
  477. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  478. C ... Calcul proprement dit de la matrice de rigidité ... C
  479. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  480.  
  481. C ... On commence par préparer le segment GL2LOC ...
  482.       IECART=NOMAXI-NOMINI+1
  483.       SEGINI,GL2LOC
  484.  
  485. C ... et le remplir ...
  486.       IF(NBNBMA.GT.0) THEN
  487.          SEGACT,MAPBMA
  488.          DO 3000 I=1,NBNBMA
  489.             IEQUIV(1,MAPBMA.NUM(1,I)-NOMINI+1)=I
  490.  3000    CONTINUE
  491.       ENDIF
  492.  
  493.       IF(NBNCOQ.GT.0) THEN
  494.          SEGACT,MAPCOQ
  495.          DO 3010 I=1,NBNCOQ
  496.             IEQUIV(2,MAPCOQ.NUM(1,I)-NOMINI+1)=NBNBMA+2*I-1
  497.             IEQUIV(3,MAPCOQ.NUM(1,I)-NOMINI+1)=NBNBMA+2*I
  498.  3010    CONTINUE
  499.       ENDIF
  500.  
  501.       IF(NBNCOM.GT.0) THEN
  502.          SEGACT,MAPCOM
  503.          DO 3020 I=1,NBNCOQ
  504.             IEQUIV(1,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I-2
  505.             IEQUIV(2,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I-1
  506.             IEQUIV(3,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I
  507.  3020    CONTINUE
  508.       ENDIF
  509.  
  510. cdebug segprt,gl2loc
  511.  
  512. C ... On détruit les MAPOIN ? dont on n'aura plus besoin ? vraiment ?...
  513. ***   IF(NBNBMA.GT.0) SEGSUP,MAPBMA
  514.       IF(NBNCOQ.GT.0) SEGSUP,MAPCOQ
  515.       IF(NBNCOM.GT.0) SEGSUP,MAPCOM
  516.  
  517. C ... On commence par extraire un segment INFO pour chaque zone
  518. C     élémentaire. On profite de cette boucle pour trouver le
  519. C     maximum des NBNN sur toutes les zones élémentaires ...
  520.  
  521.       SEGINI,LESINF
  522.       DO 3900 I=1,NBSZEL
  523.  
  524.         IMODEL=KMODEL(I)
  525.         MELE=NEFMOD
  526.         segact imodel
  527.         if( 2+inttyp.gt.infmod(/1) ) then
  528.           CALL ELQUOI(MELE,0,INTTYP,ikk,IMODEL)
  529.          info = ikk
  530.         else
  531.           jg=16
  532.           segini info
  533.           do iou=1,16
  534.           infell(iou)=infele(iou)
  535.           enddo
  536.           infell(11)=infmod(2+inttyp)
  537.         endif
  538.         LINFO(I)=info
  539.         MELEME=IMAMOD
  540.         NBNN=NUM(/1)
  541.         NBEL=NUM(/2)
  542.         IF(I.EQ.1) THEN
  543.            NBNNMX = NBNN
  544.            NBELMX = NBEL
  545.         ENDIF
  546.         IF(NBNN.GT.NBNNMX) NBNNMX = NBNN
  547.         IF(NBEL.GT.NBELMX) NBELMX = NBEL
  548.  
  549.  3900 CONTINUE
  550.  
  551. C ... NBNNMX et NBELMX servent à dimensionner le segment LESAIJ
  552. C     qui contiendra les termes suivants :
  553. C
  554. C      ik  /
  555. C     a  = | N    sur l'élément N° k de la zone élémentaire N° i,
  556. C      j   /  j   j étant l'indice du noeud de l'élément (et de
  557. C                 la fonction de forme associée)
  558. C
  559.  
  560.       SEGINI,LESAIJ
  561.       IDIM1=IDIM-1
  562.       DO 3910 I=1,NBSZEL
  563.  
  564.         IMODEL=KMODEL(I)
  565.         MELEME=IMAMOD
  566.         NBNN=NUM(/1)
  567.         NBEL=NUM(/2)
  568.  
  569.         INFO=LINFO(I)
  570.         MINTE=INFELL(11)
  571.         SEGACT,MINTE
  572. cdebug      write(*,*) '-----------------------------'
  573. cdebug      write(*,*) 'ZONE N° ',I
  574. cdebug      segprt,minte
  575.         NBNGAU=POIGAU(/1)
  576.  
  577.         DO 3915 IEL=1,NBEL
  578.  
  579.           CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XEL)
  580.  
  581. C     ... SEG2 ...
  582.           IF(ITYPEL.EQ.2) THEN
  583.              XCH(1,1)=XEL(1,1)
  584.              XCH(2,1)=XEL(2,1)
  585.              XCH(1,2)=XEL(1,2)
  586.              XCH(2,2)=XEL(2,2)
  587. C     ... SEG3 ...
  588.           ELSE IF(ITYPEL.EQ.3) THEN
  589.              XCH(1,1)=XEL(1,1)
  590.              XCH(2,1)=XEL(2,1)
  591.              XCH(1,2)=XEL(1,3)
  592.              XCH(2,2)=XEL(2,3)
  593. C     ... TRI3 ou QUA4 ...
  594.           ELSE IF(ITYPEL.EQ.4 .OR. ITYPEL.EQ.8) THEN
  595.              XCH(1,1)=XEL(1,1)
  596.              XCH(2,1)=XEL(2,1)
  597.              XCH(3,1)=XEL(3,1)
  598.              XCH(1,2)=XEL(1,2)
  599.              XCH(2,2)=XEL(2,2)
  600.              XCH(3,2)=XEL(3,2)
  601.              XCH(1,3)=XEL(1,3)
  602.              XCH(2,3)=XEL(2,3)
  603.              XCH(3,3)=XEL(3,3)
  604. C     ... TRI6 ou QUA8 ...
  605.           ELSE IF(ITYPEL.EQ.6 .OR. ITYPEL.EQ.10) THEN
  606.              XCH(1,1)=XEL(1,1)
  607.              XCH(2,1)=XEL(2,1)
  608.              XCH(3,1)=XEL(3,1)
  609.              XCH(1,2)=XEL(1,3)
  610.              XCH(2,2)=XEL(2,3)
  611.              XCH(3,2)=XEL(3,3)
  612.              XCH(1,3)=XEL(1,5)
  613.              XCH(2,3)=XEL(2,5)
  614.              XCH(3,3)=XEL(3,5)
  615.           ELSE
  616.              CALL ERREUR(16)
  617.              RETURN
  618.           ENDIF
  619.  
  620.           IF(IDIM.EQ.2)THEN
  621.              CALL VPAST2(XCH,BPSS)
  622.           ELSE IF(IDIM.EQ.3) THEN
  623.              CALL VPAST(XCH,BPSS)
  624.           ENDIF
  625.  
  626.           CALL VCORL1(XEL,XTR,BPSS,NBNN)
  627.  
  628.           DO 3920 IPG=1,NBNGAU
  629.  
  630.             DO 3930 INN=1,NBNN
  631.               SHP(1,INN)=SHPTOT(1,INN,IPG)
  632.               SHP(2,INN)=SHPTOT(2,INN,IPG)
  633.               SHP(3,INN)=SHPTOT(3,INN,IPG)
  634.  3930       CONTINUE
  635.  
  636.             CALL JACOBI(XTR,SHP,IDIM1,NBNN,DJAC)
  637. cdebug      write(*,*) 'Zone ',I,', Element ',IEL,', PG ',IPG,
  638. cdebug     &            '-> J = ',DJAC
  639.  
  640.            IF(IFOMOD.NE.0) THEN
  641.  
  642.               DO 3935 INN=1,NBNN
  643.                 AIJ(INN,IEL,I) = AIJ(INN,IEL,I) +
  644.      &                     SHPTOT(1,INN,IPG)*POIGAU(IPG)*DJAC
  645.  3935         CONTINUE
  646.  
  647.            ELSE
  648. C      ... axisymetrique : calcul du rayon au point d'integration
  649. C
  650.               CALL DISTRR(XEL,SHP,NBNN,RR)
  651.               RR=RR*2.D0*XPI
  652.  
  653.               DO 3936 INN=1,NBNN
  654.                 AIJ(INN,IEL,I) = AIJ(INN,IEL,I) +
  655.      &                     SHPTOT(1,INN,IPG)*POIGAU(IPG)*RR*DJAC
  656.  3936         CONTINUE
  657.            ENDIF
  658.  
  659.  3920     CONTINUE
  660.  
  661.  3915   CONTINUE
  662.  
  663.  3910 CONTINUE
  664.  
  665. cdebug      segprt,lesaij
  666.  
  667. C ... Double boucle sur tous les éléments concernés ...
  668. C ... Tout ce qui se termine par 1 concerne la boucle extérieure :
  669. C      I1, IMODE1, IPT1, INFO1, MINTE1 ...
  670.  
  671.       SEGACT,CHAMTM
  672.       CALL TCHAMR(CHAMR,1)
  673.  
  674.       BILAN = 0.D0
  675.       DO 4000 I1=1,NBSZEL
  676.  
  677.         IMODE1=KMODEL(I1)
  678.         IPT1=IMODE1.IMAMOD
  679.         NBNN1=IPT1.NUM(/1)
  680.  
  681. C   ... On va chercher la température moyenne de l'élément IEL1 ...
  682. C   ... On commence par trouver les bons numéros de composantes,
  683. C       puis on active les MELVAL associés ...
  684.  
  685.         SEGACT,CHAMTM
  686.         MCHAM1=CHAMTM.ICHAML(I1)
  687. C       write(*,*) ' I1 MCHAM1 = ',I1,MCHAM1
  688.         SEGACT,MCHAM1
  689.         NC=MCHAM1.IELVAL(/1)
  690.         IF(KCOQ(I1)) THEN
  691.            NCTINF=0
  692.            NCTSUP=0
  693.            DO 4005 I=1,NC
  694.               IF(MCHAM1.NOMCHE(I).EQ.'TINF') NCTINF=I
  695.               IF(MCHAM1.NOMCHE(I).EQ.'TSUP') NCTSUP=I
  696.  4005      CONTINUE
  697.  
  698.            IF(NCTINF.EQ.0) GOTO 9999
  699.            MELVA1=MCHAM1.IELVAL(NCTINF)
  700.            SEGACT,MELVA1
  701.  
  702.            IF(NCTSUP.EQ.0) GOTO 9999
  703.            MELVA2=MCHAM1.IELVAL(NCTSUP)
  704.            SEGACT,MELVA2
  705.         ELSE
  706.            NCT=0
  707.            DO 4006 I=1,NC
  708.               IF(MCHAM1.NOMCHE(I).EQ.'T') NCT=I
  709.  4006      CONTINUE
  710.  
  711.            IF(NCT.EQ.0) GOTO 9999
  712.            MELVA1=MCHAM1.IELVAL(NCT)
  713.            SEGACT,MELVA1
  714.         ENDIF
  715.  
  716.         DO 4010 IEL1=1,IPT1.NUM(/2)
  717.  
  718. C     ... On trouve la (les) température(s) correspondants à
  719. C         l'élément IEL1 ...
  720.  
  721.           IF(KCOQ(I1)) THEN
  722.               IF(MELVA1.VELCHE(/2).EQ.1) THEN
  723.                 T=MELVA1.VELCHE(1,1)
  724.               ELSE
  725.                 T=MELVA1.VELCHE(1,IEL1)
  726.               ENDIF
  727.               T3(2)=T*T*T
  728.  
  729.               IF(MELVA2.VELCHE(/2).EQ.1) THEN
  730.                 T=MELVA2.VELCHE(1,1)
  731.               ELSE
  732.                 T=MELVA2.VELCHE(1,IEL1)
  733.               ENDIF
  734.               T3(1)=T*T*T
  735.           ELSE
  736.               IF(MELVA1.VELCHE(/2).EQ.1) THEN
  737.                 T=MELVA1.VELCHE(1,1)
  738.               ELSE
  739.                 T=MELVA1.VELCHE(1,IEL1)
  740.               ENDIF
  741.             T3(1)=T*T*T
  742.           ENDIF
  743.  
  744. C     ... Deuxième niveau de la boucle sur les éléments ...
  745.  
  746.           DO 4020 I2=1,NBSZEL
  747.  
  748.             IMODE2=KMODEL(I2)
  749.             IPT2=IMODE2.IMAMOD
  750.             NBNN2=IPT2.NUM(/1)
  751.  
  752.             DO 4030 IEL2=1,IPT2.NUM(/2)
  753.  
  754. C       ... On va chercher les surfaces des éléments IEL1 et IEL2 ...
  755.  
  756.               MCHAML=CHAMR.ICHAML(I2)
  757.               SEGACT,MCHAML
  758.               NC2=IELVAL(/1)
  759. C         ... Ici on suppose que la surface est toujours dernière ...
  760.               MELVAL=IELVAL(NC2)
  761.               SURF2=VELCHE(1,IEL2)
  762.  
  763.               MCHAML=CHAMR.ICHAML(I1)
  764.               SEGACT,MCHAML
  765.               NC1=IELVAL(/1)
  766. C         ... Ici on suppose que la surface est toujours dernière ...
  767.               MELVAL=IELVAL(NC1)
  768.               SURF1=VELCHE(1,IEL1)
  769.  
  770. C         ... On a supposé que dans le champ de facteurs de forme les
  771. C             composantes sont toujours dans l'ordre suivant :
  772. C             MIDL SURF       dans le cas des solides
  773. C             SUPE INFE SURF  dans le cas des coques
  774. C
  775.               DO 4040 IC1=1,NC1-1
  776.                  MELVAL=IELVAL(IC1)
  777.                  MCHELM=IELCHE(1,IEL1)
  778.                  MCHAML=ICHAML(I2)
  779.                  DO 4050 IC2=1,NC2-1
  780.                     MELVAL=IELVAL(IC2)
  781.                     COEFF=VELCHE(1,IEL2)
  782. C             write(*,*)' iel1 iel2 coeff ',IEL1,IEL2,COEFF
  783.  
  784.                     DO 4060 INN1=1,NBNN1
  785.                        NGL1=IPT1.NUM(INN1,IEL1)
  786.                        NDDL1=IC1
  787.                        IF(KCOQ(I1)) NDDL1=NDDL1+1
  788.                        NLOC1 = IEQUIV(NDDL1,NGL1-NOMINI+1)
  789.                        DO 4070 INN2=1,NBNN2
  790.  
  791.                           NGL2=IPT2.NUM(INN2,IEL2)
  792.                           NDDL2=IC2
  793.                           IF(KCOQ(I2)) NDDL2=NDDL2+1
  794.                           NLOC2 = IEQUIV(NDDL2,NGL2-NOMINI+1)
  795.  
  796. C             write(*,*)' in1 in2  coeff ',INN1,INN2,COEFF
  797.                           TERME = COEFF*AIJ(INN1,IEL1,I1)
  798.      &      *AIJ(INN2,IEL2,I2)*T3(IC1)*CONSSB/(SURF1*SURF2)
  799.                           RE(NLOC2,NLOC1,1)=RE(NLOC2,NLOC1,1) + TERME
  800.                           BILAN = BILAN + TERME
  801.  
  802. cdebug      write(*,*) 'On rajoute ',TERME,' au terme ',NLOC1,NLOC2
  803.  
  804.  4070                  CONTINUE
  805.  
  806.  4060               CONTINUE
  807.  
  808.  4050            CONTINUE
  809.  
  810.  4040         CONTINUE
  811.  
  812.  4030       CONTINUE
  813.  
  814.  4020     CONTINUE
  815.  
  816.  4010   CONTINUE
  817.  
  818.  4000 CONTINUE
  819.       IF (IIMPI.GE.4) THEN
  820.       write(*,*) 'Bilan = ',BILAN
  821.       ENDIF
  822.  
  823. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  824. C ... Sortie de la rigidité calculée ...                         C
  825. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  826.  
  827. C!!   segprt,xmatri
  828.       SEGDES,XMATRI
  829.       CALL ECROBJ('RIGIDITE',MRIGID)
  830.  
  831. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  832. C ... Le ménage à la fin ...                                     C
  833. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  834.  
  835. C ... On désactive le modèle ...
  836.       DO 5000 I=1,NBSZEL
  837.          IMODEL=KMODEL(I)
  838.          MELEME=IMAMOD
  839.          INFO=LINFO(I)
  840.          SEGSUP,INFO
  841.  5000 CONTINUE
  842.       SEGSUP,LESINF
  843.  
  844. C ... desactivation de la matrice de rayonnement
  845.       CALL TCHAMR(CHAMR,0)
  846.  
  847. C ... On supprime les segments de travail ...
  848.       SEGSUP,INFOEL
  849.       SEGSUP,GL2LOC
  850.       SEGSUP,LESAIJ
  851.  
  852.       RETURN
  853.  
  854.  9999 CONTINUE
  855.       CALL ERREUR(26)
  856.       END
  857.  
  858.  
  859.  
  860.  
  861.  
  862.  
  863.  

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