Esope source files

   1 : C THNUMAC2  SOURCE    FANDEUR   11/07/19    21:17:49     7042           
   2 : 
   3 : C=======================================================================
   4 : C=                            T H N U M A C 2                          =
   5 : C=                            ---------------                          =
   6 : C=                      (TNUMAC dans le cas de la thermique)           =
   7 : C=  Fonction :                                                         =
   8 : C=  ----------                                                         =
   9 : C=  Calcul de la matrice de CAPACITE thermohydrique pour les elements  =
  10 : C=  finis MASSIFs a integration NUMERIQUE                              =
  11 : C=                                                                     =
  12 : C=  Parametres :  (E)=Entree  (S)=Sortie                               =
  13 : C=  ------------                                                       =
  14 : C=   NEF      (E)   Numero de l'ELEMENT FINI dans NOMTP (cf. CCHAMP)   =
  15 : C=   IMAIL    (E)   Numero du segment IMODEL dans le segment MMODEL    =
  16 : C=   IPMODE   (E)   Pointeur sur un segment IMODEL suppose ACTIF       =
  17 : C=   IPCHEM   (E)   Pointeur sur un segment MCHELM de CARACTERISTIQUES =
  18 : C=   IPRIGI  (E/S)  Pointeur sur l'objet RIGIDITE (CONDUCTIVITE)       =
  19 : C=                                                                     =
  20 : C=   Zakaria HABIBI le 30 juin 2008.                                   =
  21 : C=======================================================================
  22 : 
  23 :       SUBROUTINE THNUMAC2 (NEF,ipmail,ipinte,ipint1,IVAMAT,NMATT,
  24 :      &                     ipmatr,LRE)
  25 : 
  26 :       IMPLICIT INTEGER(I-N)
  27 :       IMPLICIT REAL*8 (A-H,O-Z)
  28 : 
  29 : -INC CCOPTIO
  30 : -INC CCREEL
  31 : 
  32 : -INC SMCHAML
  33 : -INC SMCOORD
  34 : -INC SMELEME
  35 : -INC SMINTE
  36 : -INC SMRIGID
  37 : 
  38 :       SEGMENT MMAT1
  39 :         REAL*8 VALMAT(NV1)
  40 :         REAL*8 CEL1(NBNN,NBNN),CEL2(NBNN,NBNN),CEL3(NBNN,NBNN)
  41 :         REAL*8 CEL4(NBNN,NBNN),CEL5(NBNN,NBNN),CEL6(NBNN,NBNN)
  42 :         REAL*8 CEL7(NBNN,NBNN),CEL8(NBNN,NBNN),CEL9(NBNN,NBNN)
  43 :         REAL*8 XE(3,NBNN)
  44 :         REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN),FORME(NBNN)
  45 :         REAL*8 CMAT(NDIM,NDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM)
  46 :       ENDSEGMENT
  47 : 
  48 :       SEGMENT MAXE
  49 :         REAL*8 TXR(IDIM,IDIM),XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM)
  50 :       ENDSEGMENT
  51 : 
  52 :       SEGMENT MPTVAL
  53 :         INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU)
  54 :         CHARACTER*16 TYVAL(NCOSOU)
  55 :       ENDSEGMENT
  56 : 
  57 : C INITIALISATIONS ET VERIFICATIONS
  58 : C ================================
  59 : C Recuperation d'informations sur le maillage elementaire
  60 : C =====
  61 :       MELEME = IPMAIL
  62 : c*      SEGACT,MELEME
  63 :       NBNN   = NUM(/1)
  64 :       NBELEM = NUM(/2)
  65 : C =====
  66 : C Recuperation d'informations sur le maillage elementaire
  67 : C =====
  68 :       MINTE = ipinte
  69 : c*      SEGACT,MINTE
  70 :       NBPGAU = POIGAU(/1)
  71 : C =====
  72 : C Recuperation des fonctions de forme et de leurs derivees au
  73 : C centre de gravite de l'element pour le calcul des axes locaux
  74 : C d'orthotropie ou d'anisotropie
  75 : C =====
  76 :       IF (ipint1.NE.0) THEN
  77 :         MINTE1 = ipint1
  78 : c*        SEGACT,MINTE1
  79 :         NBSH = MINTE1.SHPTOT(/2)
  80 :       ENDIF
  81 : C =====
  82 : C Initialisation des segments de travail
  83 : C =====
  84 :       MPTVAL = IVAMAT
  85 :       IF (IFOMOD.EQ.1) THEN
  86 :         NDIM=3
  87 :       ELSE
  88 :         NDIM=IDIM
  89 :       ENDIF
  90 :       NV1 = NMATT
  91 :       SEGINI,MMAT1
  92 :       MAXE = 0
  93 :       IF (ipint1.NE.0) THEN
  94 :         SEGINI,MAXE
  95 :       ENDIF
  96 : C =====
  97 : C Matrice de CAPACITE thermohydrique
  98 : C =====
  99 :       XMATRI = ipmatr
 100 : c*      SEGACT,XMATRI*MOD
 101 : c*      NLIGRP = 3*NBNN = LRE
 102 : c*      NLIGRD = 3*NBNN = LRE
 103 : 
 104 : C BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IPMAIL
 105 : C ======================================================
 106 :       DO iElt=1,NBELEM
 107 : C ===
 108 : C Mise a zero de la matrice de CAPACITE de l'element iElt
 109 : C ===
 109 :         CALL ZERO(CEL1,NBNN,NBNN)
 109 :         CALL ZERO(CEL2,NBNN,NBNN)
 109 :         CALL ZERO(CEL3,NBNN,NBNN)
 109 :         CALL ZERO(CEL4,NBNN,NBNN)
 109 :         CALL ZERO(CEL5,NBNN,NBNN)
 109 :         CALL ZERO(CEL6,NBNN,NBNN)
 109 :         CALL ZERO(CEL7,NBNN,NBNN)
 109 :         CALL ZERO(CEL8,NBNN,NBNN)
 109 :         CALL ZERO(CEL9,NBNN,NBNN)
 110 : C ===
 111 : C Recuperation des coordonnees GLOBALES des noeuds de l'element
 112 : C ===
 112 :         CALL DOXE(XCOOR,IDIM,NBNN,NUM,iElt,XE)
 113 : C ===
 114 : C Calculs des axes locaux d'orthotropie ou d'anisotropie
 115 : C ===
 116 :         IF (ipint1.NE.0) THEN
 116 :           CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR)
 117 :           IF (NBSH.EQ.-1) THEN
 117 :             CALL ERREUR(525)
 118 :             GOTO 9990
 119 :           ENDIF
 120 :         ENDIF
 121 : C ===
 122 : C Boucle sur les points de Gauss de l'element iElt
 123 : C ===
 124 :         iFois=0
 125 :         DO iGau=1,NbPGau
 126 : C - Calcul du jacobien, des fonctions de forme et de leurs
 127 : C   derivees au point de Gauss iGau
 127 :           CALL CAPA4(NEF,iGau,NBNN,XE,SHPTOT,SHP,FORME,DJAC)
 128 :           IF (IERR.NE.0) GOTO 9990
 129 :           IF (DJAC.LT.XZero) iFois=iFois+1
 130 :           DJAC=ABS(DJAC)*POIGAU(iGau)
 131 : C - Erreur si le jacobien est nul en ce point de Gauss
 132 :           IF (DJAC.LT.XPetit) THEN
 133 :             INTERR(1)=iElt
 133 :             CALL ERREUR(259)
 134 :             GOTO 9990
 135 :           ENDIF
 136 : C - Recuperation de la ou des valeurs de conductibilite au point
 137 : C   de Gauss iGau (tableau VALMAT)
 138 :           DO i=1,NMATT
 139 :             IF (IVAL(i).NE.0) THEN
 140 :               MELVAL=IVAL(i)
 141 :               IBMN=MIN(iElt,VELCHE(/2))
 142 :               IGMN=MIN(iGau,VELCHE(/1))
 143 :               VALMAT(i)=VELCHE(IGMN,IBMN)
 144 :             ELSE
 145 :               VALMAT(i)=XZERO
 146 :             ENDIF
 147 :           ENDDO
 148 : C - Cas d'un materiau ISOTROPE de conductibilite K
 149 : C   Calcul de la contribution du point de Gauss a la matrice
 150 : C   CAPACITE elementaire de cet element fini
 151 : C   Ajout du terme XK*transposee(N)*N
 152 : C   Seul cas valide en dimension 1
 153 :              XK=VALMAT(10)*DJAC
 153 :              CALL NTNST(FORME,XK,NBNN,1,CEL1)
 154 :              XK=VALMAT(11)*DJAC
 154 :              CALL NTNST(FORME,XK,NBNN,1,CEL2)
 155 :              XK=VALMAT(12)*DJAC
 155 :              CALL NTNST(FORME,XK,NBNN,1,CEL3)
 156 :              XK=VALMAT(13)*DJAC
 156 :              CALL NTNST(FORME,XK,NBNN,1,CEL4)
 157 :              XK=VALMAT(14)*DJAC
 157 :              CALL NTNST(FORME,XK,NBNN,1,CEL5)
 158 :              XK=VALMAT(15)*DJAC
 158 :              CALL NTNST(FORME,XK,NBNN,1,CEL6)
 159 :              XK=VALMAT(16)*DJAC
 159 :              CALL NTNST(FORME,XK,NBNN,1,CEL7)
 160 :              XK=VALMAT(17)*DJAC
 160 :              CALL NTNST(FORME,XK,NBNN,1,CEL8)
 161 :              XK=VALMAT(18)*DJAC
 161 :              CALL NTNST(FORME,XK,NBNN,1,CEL9)
 162 : 
 163 :         ENDDO
 164 : C =====
 165 : C Erreur si, en un point de Gauss, le jacobien change de signe
 166 : C =====
 167 :         IF (iFois.NE.0.AND.iFois.NE.NbPGau) THEN
 168 :           INTERR(1)=iElt
 168 :           CALL ERREUR(195)
 169 :           GOTO 9990
 170 :         ENDIF
 171 : C =====
 172 : C Stockage de la matrice de CAPACITE pour cet element fini
 173 : C =====  (remplissage de XMATRI)
 173 :        CALL REMPMCH
 174 :      & (CEL1,CEL2,CEL3,CEL4,CEL5,CEL6,CEL7,CEL8,CEL9,NBNN,RE(1,1,ielt))
 175 : 
 176 :       ENDDO
 177 :         
 178 : C  3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS
 179 : C ====================================================
 180 :  9990 CONTINUE
 181 :       SEGSUP,MMAT1
 182 :       IF (ipint1.GT.0) THEN
 183 :         SEGSUP,MAXE
 184 :       ENDIF
 185 : 
 186 :       RETURN
 187 :       END
 188 : 
 189 :