Sources Esope | Cast3M

hhoelq
C HHOELQ    SOURCE    OF166741  26/02/23    21:15:02     12480          
 
C===== FORMULATION HHO =================================================
C=                  Sous-programme HHOELQ                              =
C=                  ---------------------                              =
C= Ce sous programme est l'equivalent des sous programmes PRQUOI et    =
C= ELQUOI pour la formulation HHO. (Il s'en inspire fortement.)        =
C=                                                                     =
C= Entree :                                                            =
C= --------                                                            =
C=   IPMODL      pointeur sur une zone elementaire d'un modele         =
C=               SEGMENT ACTIF en entree et en sortie                  =
C=                                                                     =
C= Sortie :                                                            =
C= --------                                                            =
C=   Le pointeur IPMODL mis a jour                                     =
C=======================================================================
 
      SUBROUTINE HHOELQ (IPMODL)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
 
-INC CCHHOPA
-INC CCHHOPR
 
-INC SMMODEL
-INC SMELEME
-INC SMINTE
-INC SMLENTI
 
C=======================================================================
C= INTHHO : SAVE local des pointeurs sur les segments d'integration    =
C= NBMODE : Nombre de MODEs de calcul (= nb valeurs possibles IFOMOD)  =
C= NINFOS : Dimension (JG) du tableau infele dans IMODEL = 16          =
C=======================================================================
C IFOMOD :  -1     0     1     2     3     4     5     6
C          PLAN  AXIS  FOUR  TRID  UNID  UNID  UNID  FREQ
C                                  PLAN  AXIS  SPHE
C NBMODE = IFOMOD + 2 mais on ne considere pas IFOMOD = 6
C NTYNTE : 1 --> INTEGRATION aux NOEUDS
C          2 --> INTEGRATION aux Points d'INTEGRATION
C          3 --> INTEGRATION au centre de GRAVITE (champ constant par element)
C NOFACE : pour ordre faces (+1 car ordre 0 accepte)
C NNFACE : nombre de noeuds decrivant une face (2 en 2D, 3 a HHO_MAX_EDGE en 3D)
C NOCELL : pour ordre cellules (+1 car ordre 0 accepte)
C NNCELL : nombre de noeuds decrivant une cellule (3 a HHO_MAX_EDGE en 2D)
C          pour le 3D voir la remarque ci-dessous :
C Il faudra revoir ce stockage quand on aura du 3D
C en 2D --> on a des polygones definis par leur nombre de noeuds, ce 
C           qui les rend unique !
C=======================================================================
 
      PARAMETER ( NTYNTE = 3 , NBMODE = 7 )
      PARAMETER ( NOCELL = 6 , NNCELL = NFAMAX )
      PARAMETER ( NOFACE = 6 , NNFACE = NFAMAX )
      PARAMETER ( NINTEG = NNCELL*NOCELL*NTYNTE*NBMODE )
      PARAMETER ( NINFOS = 16 )
 
      CHARACTER*(16) motHHO
      LOGICAL b_z
 
      INTEGER INTHHO(NNCELL,NOCELL,NTYNTE,NBMODE)
      SAVE INTHHO
      DATA INTHHO / NINTEG*0 /
 
C On se place sur la bonne "tranche" du tableau INTHHO :
      IBMODE = IFOMOD + 2
 
      IF (IBMODE.LT.1 .OR. IBMODE.GT.NBMODE) THEN
        write(IOIMP,*)
        write(IOIMP,*) 'HHOELQ : valeur de NBMODE incorrecte !'
        write(IOIMP,*) '         IBMODE = ',IBMODE,NBMODE,IFOMOD
        CALL ERREUR(5)
        RETURN
      ENDIF
 
      CALL HHONOB(IPMODL, nobHHO)
      IF (nobHHO.LE.0) THEN
        write(ioimp,*) 'HHOELQ : IPMODL incorrect (not HHO)',nobHHO
        CALL ERREUR(5)
        RETURN
      END IF
 
      imodel = IPMODL
C  Recuperation des donnees de infele en entree
      MELE   = imodel.NEFMOD
      MFR    = NUMMFR(MELE)
      meleme = imodel.IMAMOD
      IELE   = meleme.ITYPEL
 
      mlenti = imodel.ivamod(nobhho+1)
c*      segact,mlenti
 
C  Ordre des faces et des cellules :
      n_o_face = mlenti.lect(2)
      n_o_cell = mlenti.lect(4)
 
      if (iimpi.eq.1972)
     &  write(ioimp,*) 'HHOELQ : ordre face/cell',n_o_face,n_o_cell
 
      IF (n_o_face .GT. NOFACE) THEN
        write(ioimp,*) 'HHOELQ : NOFACE incorrect'
        CALL ERREUR(5)
        RETURN
      END IF
      IF (n_o_cell .GT. NOCELL) THEN
        write(ioimp,*) 'HHOELQ : NOCELL incorrect'
        CALL ERREUR(5)
        RETURN
      END IF
 
C=---------------------------------------------------------------------=
C= REMPLISSAGE DU TABLEAU infele du IMODEL                             =
C=---------------------------------------------------------------------=
C=     La dimension de infele est fixee par NINFOS                     =
C=     infele( 1)  numero de l'element fini                            =
C=     infele( 2)  nombre de points d'integration en contraintes       =
C=                 multicouche                                         =
C=     infele( 3)  nombre de points d'integration pour la masse        =
C=     infele( 4)  nombre de points d'integration pour SIGMA           =
C=                 BSIGMA et KSIGMA                                    =
C=     infele( 5)  nombre de caracteristiques                          =
C=     infele( 6)  nombre de points d'integration pour la RIGIDITE     =
C=     infele( 7)  longueur d'un tableau de travail pour l'element     =
C=     infele( 8)  nombre de fonctions de forme                        =
C=     infele( 9)  nbre de d.d.l. dans la matrice de RIGIDITE          =
C=     infele(10)  taille de la matrice de Hooke                       =
C=     infele(11)  pointeur sur le segment d'integration               =
C=     infele(12)  pointeur sur le 2nd segment d'integration des       =
C=                 elements homogeneises ou fluides (COQ6 ou COQ8)     =
C=  infele(13) : numero de la FORMULATION - MFR                        =
C=      = HHO_MFR_ELEMENT (include CHHOPA) pour formulation HHO        =
C=  infele(14) : numero de l'element GEOMETRIQUE - IELE                =
C=     infele(15)  nombre maximal de d.d.l. par noeud                  =
C=     infele(16)  nombre de composantes de contraintes ou de deform.  =
C=---------------------------------------------------------------------=
 
C  Remplissage de infele(1) : numero de l'element fini
C  --------------------------
      imodel.infele(1) = MELE
 
C  Remplissage de infele(8) :
C  --------------------------
c faces      imodel.infele(8) = -1 * imodel.INFMOD( 9)
      imodel.infele(8) = -1 * imodel.INFMOD(12)
 
C  Remplissage de infele(13) et infele(14) :
C  -----------------------------------------
      imodel.infele(13) = MFR
      imodel.infele(14) = IELE
 
C  Remplissage de infele(10), infele(15) et infele(16) :
C  -----------------------------------------------------
      IF (IFOUR.EQ.2) THEN
        imodel.infele(10) = 6
        imodel.infele(16) = 6
      ELSE
        imodel.infele(10) = 4
        imodel.infele(16) = 4
      ENDIF
      idifo = mlenti.lect(9)
      imodel.infele(15) = idifo * MAX(mlenti.lect(3),mlenti.lect(5))
 
      imodel.infele(2) = 0
 
      NBPINT = mlenti.lect(8)
      imodel.infele(3) = NBPINT
      imodel.infele(4) = NBPINT
      imodel.infele(6) = NBPINT
C= infele(5) = Nombre de caracteristiques : utilite ?
      imodel.infele(5) = 4
C= infele(7) = Taille du segment de travail a recuperer !
      imodel.infele(7) = 100
 
      imodel.infele(9) = mlenti.lect(11)
 
      imodel.infele(12) = -9
 
      nb_faces = mlenti.lect(7)
 
C  Remplissage de infele(11) : INTEGRATION DE L'ELEMENT FINI
C  ---------------------------
C On remplit tous les segments, dans infmod(3 ...), dont le type est : =
C  1 = Champ aux noeuds
C  2 = Point de Gauss, centre de gravite et champ CONSTANT
C  3 = Point de Gauss pour la rigidite
C  4 = Point de Gauss pour la masse
C  5 = Point de Gauss - calcul des contraintes
 
      iin = 0
      DO itg = 1, 3
C  Si le segment d'integration a deja ete rempli : iin est non nul ici !
C  Si le segment n'a pas deja ete rempli : iin = 0
        iin = INTHHO(nb_faces,n_o_cell+1,itg,IBMODE)
        IF (iin.EQ.0) THEN
          minte = 0
**          nbno = nb_faces
          nbno = 0
          IF (itg.EQ.1) THEN
            nbpgau = nb_faces
            SEGINI,minte
          ELSE IF (itg.EQ.2) THEN
            nbpgau = 1
            SEGINI,minte
          ELSE IF (itg.EQ.3) THEN
            nbpgau = NBPINT
            SEGINI,minte
          ENDIF
          DO i = 1, nbpgau
            minte.POIGAU(i) = 1.D0 / REAL(nbpgau)
          END DO
          IPT1 = minte
          IF (IPT1.NE.0) CALL SAVSEG(IPT1)
          iin = IPT1
          INTHHO(nb_faces,n_o_cell+1,itg,IBMODE) = iin
        END IF
        imodel.infmod(2+itg) = iin
      END DO
 
      DO itg = 4,5
        imodel.infmod(2+itg) = imodel.infmod(2+3)
      END DO
 
C- On met dans infele(11) le dernier segment rempli (=infmod(2+3))
      imodel.infele(11) = iin
      MINTE = iin
      IF (MINTE.GT.0) SEGACT,MINTE
 
c      RETURN
      END