Télécharger la plaquette Cast3M


fr en

Télécharger hhomat.eso

Retour à la liste

Numérotation des lignes : 

hhomat

  1. C HHOMAT    SOURCE    OF166741  24/05/06    21:15:12     11082          
  2. C HHOMAT    SOURCE                                                      
  3.  
  4.       SUBROUTINE HHOMAT (IPMODL, IPCARA, iret)
  5.  
  6.       IMPLICIT INTEGER(I-N)
  7.       IMPLICIT REAL*8(A-H,O-Z)
  8.  
  9. -INC PPARAM
  10. -INC CCOPTIO
  11.  
  12. -INC CCHHOPA
  13. -INC CCHHOPR
  14.  
  15. -INC SMCOORD
  16. -INC SMELEME
  17. -INC SMMODEL
  18. -INC SMCHAML
  19. -INC SMLENTI
  20. -INC SMLREEL
  21.  
  22.       EXTERNAL LONG
  23.  
  24.       CHARACTER*(LOCHAI) chaHHO, charerr
  25.       CHARACTER*(NCONCH) conHHO
  26.       CHARACTER*(LOCOMP) mocomp
  27.  
  28.       iret = 0
  29.  
  30.       mmodel = IPMODL
  31. c*      segact,mmodel*nomod (segment actif en entree)
  32.       NSOUM = mmodel.kmodel(/1)
  33.  
  34.       mchelm = IPCARA
  35.       N1 = mchelm.IMACHE(/1)
  36. c*      segact,mchelm
  37.  
  38.       DO im = 1, NSOUM
  39.         imodel = mmodel.kmodel(im)
  40.  
  41.         IF (imodel.nefmod .NE. HHO_NUM_ELEMENT) GOTO 100
  42.         CALL HHONOB(imodel, nobHHO, iret)
  43.         IF (nobHHO.LE.0) THEN
  44.           write(ioimp,*) 'HHOMAT: nobHHO undefined'
  45.           CALL ERREUR(5)
  46.           RETURN
  47.         END IF
  48.  
  49.         iva = imodel.IVAMOD(nobHHO)
  50.         CALL QUEVAL(iva,'MOT     ',iret,lchho,r_z,chaHHO,b_z,i_z)
  51.         if (iret .NE. 0) CALL ERREUR(5)
  52.  
  53.         imaHHO = imodel.IMAMOD
  54.         CALL PLACE2(mchelm.IMACHE(1),N1,i_z,imaHHO)
  55.         IF (i_z .EQ. 0) GOTO 100
  56.         conHHO = imodel.CONMOD
  57.         IF (conHHO.NE.mchelm.CONCHE(i_z)) GOTO 100
  58.  
  59.         ISUPM = mchelm.INFCHE(i_z,6)
  60.         IF (ISUPM.NE.3) THEN
  61.           write(ioimp,*) 'HHOMAT: ISUPM /= 3 (',ISUPM,')'
  62.         END IF
  63.         IF (mchelm.INFCHE(i_z,4).NE.imodel.infmod(2+ISUPM)) THEN
  64.           write(ioimp,*) 'HHOMAT: Pointeur infche/infmod ???'
  65.         END IF
  66.  
  67.         mchaml = mchelm.ICHAML(i_z)
  68. c*        segact,mchaml
  69.         noCOMP = mchaml.IELVAL(/1)
  70.         moCOMP = 'STAB    '
  71.         CALL PLACE(mchaml.NOMCHE(1),noCOMP,i_z,moCOMP)
  72. c* pas de composante STAB, je plante tout !
  73.         IF (i_z .EQ. 0) THEN
  74.           write(ioimp,*) 'STAB component is needed for HHO formulation'
  75.           CALL ERREUR(21)
  76.           GOTO 100
  77.         END IF
  78.  
  79.         meleme = imaHHO
  80. c*        segact,meleme
  81.         NBNOE = meleme.num(/1)
  82.         NBELT = meleme.num(/2)
  83.  
  84.         mlenti = imodel.ivamod(nobhho+1)
  85. c*        segact,mlenti
  86.  
  87.         NBFAC = mlenti.lect(7)
  88. c- Le nombre de points d'integration correspond a ISUPM !
  89.         NBPIN = imodel.infele(6)
  90.  
  91.         mlent3 = imodel.ivamod(nobhho+3)
  92. c*        segact,mlent3
  93.         nb3 = mlent3.lect(/1) / 2
  94.         nbfac3 = 2 * NBFAC
  95.  
  96. c-dbg        write(ioimp,*) 'HHOMAT',im,NSOUM,nbnoe,nbelt,nbpin
  97.         IF (IDIM .NE. mlenti.lect(1)) then
  98.           write(ioimp,*) 'HHOMAT: Bizarre idim'
  99.         END IF
  100.         IF (NBNOE .NE. mlenti.lect(6)) then
  101.           write(ioimp,*) 'HHOMAT: Bizarre nbnoe'
  102.         END IF
  103.         IF (IDIM.EQ.2 .AND. NBFAC .NE. NBNOE) THEN
  104.           write(ioimp,*) 'HHOMAT: Bizarre nbfac2d'
  105.         END IF
  106.         IF (NBPIN .NE. mlenti.lect(8)) then
  107.           write(ioimp,*) 'HHOMAT: Bizarre nb.pint'
  108.         END IF
  109.         IF (nb3.ne.(NBELT*NBFAC)) THEN
  110.           write(ioimp,*) 'HHOMAT: Bizarre mlent3'
  111.         END IF
  112.  
  113.         N2 = noCOMP
  114.         moCOMP = 'PIHO    '
  115.         CALL PLACE(mchaml.NOMCHE(1),noCOMP,icpi,moCOMP)
  116.         IF (icpi .EQ. 0) THEN
  117.           N2 = N2 + 1
  118.           icpi = N2
  119.         END IF
  120.         moCOMP = 'BHHO    '
  121.         CALL PLACE(mchaml.NOMCHE(1),noCOMP,icmb,moCOMP)
  122.         IF (icmb .EQ. 0) THEN
  123.           N2 = N2 + 1
  124.           icmb = N2
  125.         END IF
  126.         moCOMP = 'SHHO    '
  127.         CALL PLACE(mchaml.NOMCHE(1),noCOMP,icms,moCOMP)
  128.         IF (icms .EQ. 0) THEN
  129.           N2 = N2 + 1
  130.           icms = N2
  131.         END IF
  132.         IF (N2.GT.noCOMP) THEN
  133.           segadj,mchaml
  134.         ELSE
  135.           segact,mchaml*MOD
  136.         END IF
  137.  
  138.         mchaml.NOMCHE(icpi) = 'PIHO    '
  139.         mchaml.TYPCHE(icpi) = 'REAL*8          '
  140.         N1PTEL = NBPIN
  141.         N1EL   = NBELT
  142.         N2PTEL = 0
  143.         N2EL   = 0
  144.         IF (icpi.GT.nocomp) THEN
  145.           SEGINI,melval
  146.           mchaml.IELVAL(icpi) = melval
  147.         ELSE
  148.           melval = mchaml.IELVAL(icpi)
  149.           SEGADJ,melval
  150.         END IF
  151.         IVPIHO = melval
  152.  
  153.         mchaml.NOMCHE(icmb) = 'BHHO    '
  154.         mchaml.TYPCHE(icmb) = 'POINTEURLISTREEL'
  155.         N1PTEL = 0
  156.         N1EL   = 0
  157.         N2PTEL = NBPIN
  158.         N2EL   = NBELT
  159.         IF (icmb.GT.nocomp) THEN
  160.           SEGINI,melval
  161.           mchaml.IELVAL(icmb) = melval
  162.         ELSE
  163.           melval = mchaml.IELVAL(icmb)
  164.           SEGADJ,melval
  165.         END IF
  166.         IVMBHO = melval
  167.  
  168.         mchaml.NOMCHE(icms) = 'SHHO    '
  169.         mchaml.TYPCHE(icms) = 'POINTEURLISTREEL'
  170.         melval = mchaml.IELVAL(icms)
  171.         N1PTEL = 0
  172.         N1EL   = 0
  173.         N2PTEL = 1
  174.         N2EL   = NBELT
  175.         IF (icms.GT.nocomp) THEN
  176.           SEGINI,melval
  177.           mchaml.IELVAL(icms) = melval
  178.         ELSE
  179.           melval = mchaml.IELVAL(icms)
  180.           SEGADJ,melval
  181.         END IF
  182.         IVMSTA = melval
  183.  
  184. C- Preparation des tableaux pour l'appel a HHOC3M et les fonctions de
  185. C- calcul associees a l'element HHO ici present
  186.         ne_0  = 9
  187. * connectivity des faces : en 2D 2 noeuds par aretes --> en 3D a faire !
  188.         ne_co = nbfac3
  189.         ie_0  = 0
  190.         ie_co = ie_0  + ne_0
  191.         ie_1  = ie_co + ne_co
  192.  
  193.         nr_0  = 0
  194.         nr_co = IDIM * nbnoe
  195.         nr_pi = NBPIN
  196.         nr_bh = mlenti.lect(14)
  197.         nr_sh = mlenti.lect(16)
  198.         ir_0  = 0
  199.         ir_co = ir_0  + nr_0
  200.         ir_1  = ir_co + nr_co
  201.  
  202.         JG = ie_1
  203.         SEGINI,mlent2
  204.         ile = JG
  205.  
  206.         mlent2.lect(ie_0 + 1) = ie_co + 1
  207.         mlent2.lect(ie_0 + 2) = ie_1
  208.         mlent2.lect(ie_0 + 3) = ir_co + 1
  209.         mlent2.lect(ie_0 + 4) = ir_1
  210.  
  211.         JG = nr_0 + nr_co + MAX(nr_pi,nr_bh,nr_sh)
  212.         SEGINI,mlree2
  213.         ilr = JG
  214.  
  215. c*        SEGACT,MCOORD*NOMOD
  216.  
  217.         DO IEL = 1, NBELT
  218.  
  219. C- Recuperation des coordonnees des noeuds de l element IEL
  220.           CALL HHOCOO(meleme.num,NBNOE,mcoord.xcoor(1), IEL,
  221.      &                mlree2.prog(ir_co+1),iret)
  222.           IF (iret.NE.0) RETURN
  223. c-dbg      write(ioimp,*) 'HHOMAT-HHOCOO',IEL,NBNOE,ir_co,ir_1
  224. c-dbg      write(ioimp,*) '   ',(meleme.num(i,IEL),i=1,NBNOE)
  225. c-dbg      write(ioimp,*) '  X',(mlree2.prog(ir_co+i),i=1,NBNOE)
  226. c-dbg      write(ioimp,*) '  Y',(mlree2.prog(ir_co+NBNOE+i),i=1,NBNOE)
  227. C- Connectivite des faces de l element IELT en tenant compte du signe :
  228. C- il faudra faire une fonction pour le 3D. pour le 2D que des SEG2
  229. C- Dans la bibliotheque, la numerotation locale des sommets d'une cellule
  230. C- debute a 0 ! Attention en 2D sur le dernier element !
  231.           in1 = (IEL-1) * NBFAC3
  232.           DO i = 1, NBFAC
  233.             j = 2 * i
  234.             in2 = in1 + j
  235.             if (mlent3.lect(in2-1).ne.2)
  236.      &        write(ioimp,*) 'HHOMAT : Bizarre mlent3 2',iel,i,in1,in2
  237.             ie2 = mlent3.lect(in2)
  238.             if (ie2.eq.0)
  239.      &        write(ioimp,*) 'HHOMAT : Bizarre mlent3 ie2',iel,i,in1,in2
  240.             IF (ie2.GT.0) THEN
  241.               mlent2.lect(ie_co+j-1) = i-1
  242.               mlent2.lect(ie_co+j  ) = MOD(i,NBFAC)
  243.             ELSE
  244.               mlent2.lect(ie_co+j-1) = MOD(i,NBFAC)
  245.               mlent2.lect(ie_co+j  ) = i-1
  246.             END IF
  247. c-dbg      write(ioimp,*) 'hhomat - iel',iel,i,in2,ie2
  248.           END DO
  249.  
  250. C- Recuperation des poids de Gauss :
  251.           mlent2.lect(ie_0 + 5) = nr_pi
  252.           mlent2.lect(ie_0 + 6) = ir_1
  253.           mlent2.lect(ie_0 + 7) = 0
  254.  
  255.           DO i = ir_1+1, ilr
  256.             mlree2.prog(i) = 0.D0
  257.           END DO
  258.  
  259.           iretc = 0
  260.           CALL HHOC3M('INTG',chaHHO(1:lchho),
  261.      &                   mlent2.lect,ile, mlree2.prog,ilr,
  262.      &                   iretc,charerr)
  263.           IF (iretc.NE.0) THEN
  264.             write(ioimp,*) 'HHOMAT -> HHOC3M -> INTG : ERROR ='
  265.             write(ioimp,*) charerr(1:LONG(charerr))
  266.             iret = 21
  267.             return
  268.           END IF
  269.  
  270.           melval = IVPIHO
  271.           DO IGA = 1, NBPIN
  272.             melval.velche(IGA,IEL) = mlree2.prog(ir_1+iga)
  273.           END DO
  274.  
  275. C- Recuperation des matrices B en chaque point de Gauss :
  276.           mlent2.lect(ie_0 + 5) = ir_1 + 1
  277.           mlent2.lect(ie_0 + 6) = ir_1 + nr_bh
  278.  
  279.           melval = IVMBHO
  280.           DO IGA = 1, NBPIN
  281.  
  282.             mlent2.lect(ie_0 + 7) = iga
  283.             DO i = ir_1+1, ilr
  284.               mlree2.prog(i) = 0.D0
  285.             END DO
  286.  
  287.             iretc = 0
  288.             CALL HHOC3M('BHHO',chaHHO(1:lchho),
  289.      &                         mlent2.lect,ile, mlree2.prog,ilr,
  290.      &                         iretc,charerr)
  291.             IF (iretc.NE.0) THEN
  292.               write(ioimp,*) 'HHOMAT -> HHOC3M -> BHHO : ERROR ='
  293.               write(ioimp,*) charerr(1:LONG(charerr))
  294.               iret = 21
  295.               return
  296.             END IF
  297.  
  298. c- La matrice BHHO(9,NBDDL) :
  299.             NLIG = 9
  300.             NCOL = mlenti.lect(11)
  301.             JG = nr_bh
  302.             SEGINI,mlreel
  303.       if (JG.ne.(NLIG*NCOL)) write(ioimp,*)'HHOMAT - Bizarre nr_bh',iga
  304. c* Matrice BHHO stockee ligne par ligne dans la bibliotheque.
  305. c* Matrice BHHO stockee colonne par colonne dans Cast3M.
  306.             DO j = 1, NCOL
  307.               jc = NLIG * (j-1)
  308.               DO i = 1, NLIG
  309.                 ic = NCOL * (i-1)
  310.                 mlreel.prog(jc+i) = mlree2.prog(ir_1+ic+j)
  311.               END DO
  312.             END DO
  313. c*            SEGDES,mlreel
  314.  
  315.             melval.ielche(IGA,IEL) = mlreel
  316.  
  317.           END DO
  318.  
  319. C- Recuperation de la matrice de stabilisation :
  320.           mlent2.lect(ie_0 + 5) = ir_1 + 1
  321.           mlent2.lect(ie_0 + 6) = ir_1 + nr_sh
  322.  
  323.           DO i = ir_1+1, ilr
  324.             mlree2.prog(i) = 0.D0
  325.           END DO
  326.  
  327.           iretc = 0
  328.           CALL HHOC3M('SHHO',chaHHO(1:lchho),
  329.      &                   mlent2.lect,ile, mlree2.prog,ilr,
  330.      &                   iretc,charerr)
  331.           IF (iretc.NE.0) THEN
  332.             write(ioimp,*) 'HHOMAT -> HHOC3M -> SHHO : ERROR ='
  333.             write(ioimp,*) charerr(1:LONG(charerr))
  334.             iret = 21
  335.             return
  336.           END IF
  337.  
  338.           NLIG = mlenti.lect(11)
  339.           NCOL = mlenti.lect(11)
  340.           JG = nr_sh
  341.           SEGINI,mlreel
  342.       if (JG.ne.(NLIG*NCOL)) write(ioimp,*)'HHOMAT - Bizarre nr_sh'
  343. c* !! matrice SHHO stockee ligne par ligne dans la bibliotheque.
  344. c* !! matrice SHHO stockee colonne par colonne dans Cast3M.
  345.           DO j = 1, NCOL
  346.             jc = NLIG * (j-1)
  347.             DO i = 1, NLIG
  348.               ic = NCOL * (i-1)
  349.               mlreel.prog(jc+i) = mlree2.prog(ir_1+ic+j)
  350.             END DO
  351.           END DO
  352. c*          SEGDES,mlreel
  353.  
  354.           melval = IVMSTA
  355.           melval.ielche(1,IEL) = mlreel
  356.  
  357.         END DO
  358.  
  359.  100    CONTINUE
  360.       END DO
  361.  
  362. c      RETURN
  363.       END
  364.  
  365.  
  366.  

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