Download ccgqme.eso

Back to the list

   1 : C CCGQME    SOURCE    GOUNAND   07/07/05    21:15:13     5784           
   2 :       SUBROUTINE CCGQME(LCOF,NOMLOI,
   3 :      $     FC,
   4 :      $     IMPR,IRET)
   5 :       IMPLICIT REAL*8 (A-H,O-Z)
   6 :       IMPLICIT INTEGER (I-N)
   7 : C***********************************************************************
   8 : C NOM         : CCGQME
   9 : C DESCRIPTION : Lois de comportement aux points de Gauss :
  10 : C               Qualité du maillage : alignement et isotropie
  11 : C               cf. Huang
  12 : C     
  13 : C LANGAGE     : ESOPE
  14 : C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  15 : C               mél : gounand@semt2.smts.cea.fr
  16 : C***********************************************************************
  17 : C APPELES          :
  18 : C APPELE PAR       :
  19 : C***********************************************************************
  20 : C ENTREES            :
  21 : C ENTREES/SORTIES    :
  22 : C SORTIES            : -
  23 : C TRAVAIL            :
  24 : C***********************************************************************
  25 : C VERSION    : v1, 11/05/07, version initiale
  26 : C HISTORIQUE : v1, 11/05/07, création
  27 : C HISTORIQUE :
  28 : C***********************************************************************
  29 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
  30 : C en cas de modification de ce sous-programme afin de faciliter
  31 : C la maintenance !
  32 : C***********************************************************************
  33 : -INC CCOPTIO
  34 : -INC CCREEL
  35 : CBEGININCLUDE SMCHAEL
  36 :       SEGMENT MCHAEL
  37 :       POINTEUR IMACHE(N1).MELEME
  38 :       POINTEUR ICHEVA(N1).MCHEVA
  39 :       ENDSEGMENT
  40 :       SEGMENT MCHEVA
  41 :       REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM)
  42 :       ENDSEGMENT
  43 :       SEGMENT LCHEVA
  44 :       POINTEUR LISCHE(NBCHE).MCHEVA
  45 :       ENDSEGMENT
  46 : CENDINCLUDE SMCHAEL
  47 :       INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
  48 :       POINTEUR FC.MCHEVA
  49 :       POINTEUR LCOF.LCHEVA
  50 :       POINTEUR JMAJAC.MCHEVA
  51 :       POINTEUR JMIJAC.MCHEVA
  52 :       POINTEUR JDTJAC.MCHEVA
  53 :       POINTEUR JMAREG.MCHEVA
  54 :       POINTEUR JMET.MCHEVA
  55 :       POINTEUR JTHE.MCHEVA
  56 :       POINTEUR JGAM.MCHEVA
  57 :       CHARACTER*8 NOMLOI
  58 :       INTEGER ICOF
  59 : *     
  60 : -INC TMXMAT
  61 : * Objets temporaires
  62 :       POINTEUR JAC.MXMAT,JT.MXMAT
  63 :       POINTEUR G.MXMAT,IG.MXMAT,H.MXMAT,HIG.MXMAT,GIH.MXMAT
  64 :       POINTEUR ME.MXMAT,JTM.MXMAT,MJ.MXMAT
  65 : *     
  66 :       SEGMENT MCOF
  67 :       POINTEUR COEF(IDIM,IDIM).MCHEVA
  68 :       ENDSEGMENT
  69 :       POINTEUR MET.MCOF
  70 : *     
  71 :       LOGICAL LBID
  72 :       INTEGER LAXSP
  73 :       REAL*8 DEUPI,XR
  74 :       REAL*8 XL,XM
  75 : *     
  76 :       INTEGER IMPR,IRET
  77 : *     
  78 : * Executable statements
  79 : *     
  80 :       IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgqme'
  81 : C      IF (.NOT.(IDIM.EQ.1)) THEN
  82 : C         WRITE(IOIMP,*) 'IDIM=',IDIM,' ?'
  83 : C         GOTO 9999
  84 : C      ENDIF
  85 :       NLFC=FC.VELCHE(/6)
  86 :       NPFC=FC.VELCHE(/5)
  87 :       ICOF=0 
  88 : *     
  89 : *     Récupération des coefficients de la metrique
  90 : *     
  91 :       SEGINI MET
  92 :       DO IIDIM=1,IDIM
  93 :          ICOF=ICOF+1
  94 :          JMET=LCOF.LISCHE(ICOF)
  95 :          IF (ICOF.EQ.1) THEN 
  96 :             NLJM=JMET.VELCHE(/6)
  97 :             NPJM=JMET.VELCHE(/5)
  98 :          ELSE
  99 :             NLJM2=JMET.VELCHE(/6)
 100 :             NPJM2=JMET.VELCHE(/5)
 101 :             IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN 
 102 :                WRITE(IOIMP,*) 'Erreur grave dims JMET'
 103 :                GOTO 9999
 104 :             ENDIF
 105 :          ENDIF
 106 :          MET.COEF(IIDIM,IIDIM)=JMET
 107 :       ENDDO
 108 :       DO IIDIM=1,IDIM
 109 :          NJ=IDIM-IIDIM
 110 :          IF (NJ.GE.1) THEN
 111 :             DO JIDIM=IIDIM+1,IDIM
 112 :                ICOF=ICOF+1
 113 :                JMET=LCOF.LISCHE(ICOF)
 114 :                NLJM2=JMET.VELCHE(/6)
 115 :                NPJM2=JMET.VELCHE(/5)
 116 :                IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN 
 117 :                   WRITE(IOIMP,*) 'Erreur grave dims JMET2'
 118 :                   GOTO 9999
 119 :                ENDIF
 120 :                MET.COEF(IIDIM,JIDIM)=JMET
 121 :             ENDDO
 122 :          ENDIF
 123 :       ENDDO
 124 : *     
 125 :       ICOF=ICOF+1
 126 :       JMAJAC=LCOF.LISCHE(ICOF)
 127 :       NLJA=JMAJAC.VELCHE(/6)
 128 :       NPJA=JMAJAC.VELCHE(/5)
 129 :       IREF=JMAJAC.VELCHE(/4)
 130 :       IREL=JMAJAC.VELCHE(/3)
 131 : *     
 132 :       IF (IREL.NE.IDIM) THEN
 133 :          WRITE(IOIMP,*) 'Erreur dims JMAJAC'
 134 :          GOTO 9999
 135 :       ENDIF
 136 : *     
 137 :       ICOF=ICOF+1
 138 :       ICOF=ICOF+1
 139 :       ICOF=ICOF+1
 140 :       JMAREG=LCOF.LISCHE(ICOF)
 141 :       NLJR=JMAREG.VELCHE(/6)
 142 :       NPJR=JMAREG.VELCHE(/5)
 143 :       I1  =JMAREG.VELCHE(/4)
 144 :       I2  =JMAREG.VELCHE(/3)
 145 :       IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF))
 146 :      $     THEN
 147 :          WRITE(IOIMP,*) 'Erreur dims JMAREG'
 148 :          GOTO 9999
 149 :       ENDIF
 150 : *     
 151 : * Objets temporaires et à préconditionner
 152 : *     
 153 :       LDIM1=IREL
 154 :       LDIM2=IREF
 155 :       SEGINI,JAC
 156 :       SEGINI,MJ
 157 :       LDIM1=IREF
 158 :       LDIM2=IREL
 159 :       SEGINI,JT
 160 :       SEGINI,JTM
 161 :       LDIM1=IREF
 162 :       LDIM2=IREF
 163 :       SEGINI,G
 164 :       SEGINI,IG
 165 :       SEGINI,H
 166 :       SEGINI,HIG
 167 :       SEGINI,GIH
 168 :       LDIM1=IREL
 169 :       LDIM2=IREL
 170 :       SEGINI,ME
 171 : *     
 172 : * Calcul de la métrique des éléments réguliers
 173 : *     
 173 :       CALL MAMA(JMAREG.VELCHE,IREF,IREF,
 174 :      $     'JTJ     ',H.XMAT,IREF,IREF,
 175 :      $     IMPR,IRET)
 176 :       IF (IRET.NE.0) GOTO 9999
 177 : *      SEGPRT,H
 178 : *     
 179 :       DO ILFC=1,NLFC
 180 :          IF (NLJM.EQ.1) THEN
 181 :             ILJM=1
 182 :          ELSE
 183 :             ILJM=ILFC
 184 :          ENDIF
 185 :          IF (NLJA.EQ.1) THEN
 186 :             ILJA=1
 187 :          ELSE
 188 :             ILJA=ILFC
 189 :          ENDIF
 190 :          DO IPFC=1,NPFC
 191 :             IF (NPJM.EQ.1) THEN
 192 :                IPJM=1
 193 :             ELSE
 194 :                IPJM=IPFC
 195 :             ENDIF
 196 :             IF (NPJA.EQ.1) THEN
 197 :                IPJA=1
 198 :             ELSE
 199 :                IPJA=IPFC
 200 :             ENDIF
 201 : *     
 202 : * Copie des coefficients de la métrique
 203 : *     
 204 :             DO IIDIM=1,IDIM
 205 :                JMET=MET.COEF(IIDIM,IIDIM)
 206 :                ME.XMAT(IIDIM,IIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
 207 :             ENDDO
 208 :             DO IIDIM=1,IIDIM
 209 :                NJ=IDIM-IIDIM
 210 :                IF (NJ.GE.1) THEN
 211 :                   DO JIDIM=IIDIM+1,IDIM
 212 :                      JMET=MET.COEF(IIDIM,JIDIM)
 213 :                      ME.XMAT(IIDIM,JIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
 214 :                      ME.XMAT(JIDIM,IIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
 215 :                   ENDDO
 216 :                ENDIF
 217 :             ENDDO
 218 : *     
 219 : * Copie du jacobien
 220 : *     
 220 :             CALL MAMA(JMAJAC.VELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF,
 221 :      $           'COPIE   ',
 222 :      $           JAC.XMAT,IREL,IREF,
 223 :      $           IMPR,IRET)
 224 :             IF (IRET.NE.0) GOTO 9999
 225 : *            SEGPRT,JAC
 226 : *     
 227 : * Calcul de la métrique G
 228 : *     
 229 : *     Calcul de Jt
 229 :             CALL MAMA(JAC.XMAT,IREL,IREF,
 230 :      $           'TRANSPOS',JT.XMAT,IREF,IREL,
 231 :      $           IMPR,IRET)
 232 :             IF (IRET.NE.0) GOTO 9999
 233 : *     Calcul de MJ
 233 :             CALL MAMAMA(ME.XMAT,IREL,IREL,JAC.XMAT,IREL,IREF,
 234 :      $           'FOIS    ',MJ.XMAT,IREL,IREF,
 235 :      $           IMPR,IRET)
 236 :             IF (IRET.NE.0) GOTO 9999
 237 : *     Calcul de JTM
 237 :             CALL MAMA(MJ.XMAT,IREL,IREF,
 238 :      $           'TRANSPOS',JTM.XMAT,IREF,IREL,
 239 :      $           IMPR,IRET)
 240 :             IF (IRET.NE.0) GOTO 9999
 241 : *     Calcul de G=JtMJ
 241 :             CALL MAMAMA(JT.XMAT,IREF,IREL,MJ.XMAT,IREL,IREF,
 242 :      $           'FOIS    ',G.XMAT,IREF,IREF,
 243 :      $           IMPR,IRET)
 244 :             IF (IRET.NE.0) GOTO 9999
 245 : * Calcul de l'inverse, du déterminant et trace de l'inverse de g
 246 :             LBID=.FALSE.
 246 :             CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,LBID,IMPR,IRET)
 247 :             IF (IRET.NE.0) GOTO 9999
 248 : * Calcul de hg-1, de sa trace et de son déterminant
 248 :             CALL MAMAMA(H.XMAT,IREF,IREF,IG.XMAT,IREF,IREF,
 249 :      $           'FOIS    ',HIG.XMAT,IREF,IREF,IMPR,IRET)
 250 :             IF (IRET.NE.0) GOTO 9999
 251 :             LBID=.FALSE.
 251 :             CALL GEOLI2(IREF,1,1,HIG.XMAT,GIH.XMAT,DETHIG,
 252 :      $           LBID,IMPR,IRET)
 253 :             IF (IRET.NE.0) GOTO 9999
 254 :             XM=1.D0/SQRT(DETHIG)
 254 :             CALL MARE(GIH.XMAT,IREF,IREF,'TRACE   ',
 255 :      $           XL,IMPR,IRET)
 256 :             IF (IRET.NE.0) GOTO 9999
 257 : *     
 258 : * Calcul des qualités de maillage
 259 : *     
 260 :             IF (NOMLOI.EQ.'QEQU    ') THEN
 261 :                CONTRI=XM
 262 :             ELSEIF (NOMLOI.EQ.'QALI') THEN
 263 :                IF (IREF.EQ.1) THEN
 264 :                   CONTRI=1.D0
 265 :                ELSE
 266 :                   XIREF=DBLE(IREF)
 267 :                   XNUM=XL
 268 :                   XDEN=XIREF*(XM**(2.D0/XIREF))
 269 :                   XEXP=XIREF/(2.D0*(XIREF-1.D0))
 270 :                   CONTRI=(XNUM/XDEN)**XEXP
 271 : C                  SEGPRT,GIH
 272 : C                  WRITE(IOIMP,*) 'TRGIH =',XL
 273 : C                  WRITE(IOIMP,*) 'DETGIH=',(1.D0/DETHIG)
 274 : C                  WRITE(IOIMP,*) 'XNUM=',XNUM
 275 : C                  WRITE(IOIMP,*) 'XDEN=',XDEN
 276 :                ENDIF
 277 :             ELSE
 278 :                WRITE(IOIMP,*) 'Erreur grave'
 279 :                GOTO 9999
 280 :             ENDIF
 281 :             FC.VELCHE(1,1,1,1,IPFC,ILFC)=
 282 :      $           FC.VELCHE(1,1,1,1,IPFC,ILFC)+
 283 :      $           CONTRI
 284 :          ENDDO
 285 :       ENDDO
 286 :       SEGSUP,ME
 287 :       SEGSUP,GIH
 288 :       SEGSUP,HIG
 289 :       SEGSUP,H
 290 :       SEGSUP,IG
 291 :       SEGSUP,G
 292 :       SEGSUP,JTM
 293 :       SEGSUP,JT
 294 :       SEGSUP,MJ
 295 :       SEGSUP,JAC
 296 :       SEGSUP,MET
 297 : *     
 298 : * Normal termination
 299 : *     
 300 :          IRET=0
 301 :          RETURN
 302 : *     
 303 : * Format handling
 304 : *     
 305 : *     
 306 : * Error handling
 307 : *     
 308 :  9999    CONTINUE
 309 :          IRET=1
 310 :          WRITE(IOIMP,*) 'An error was detected in subroutine ccgqme'
 311 :          RETURN
 312 : *     
 313 : * End of subroutine CCGQME
 314 : *     
 315 :          END
 316 : 
 317 :  
 318 :  
 319 :  
 320 :  

© Cast3M 2003 - All rights reserved.
Disclaimer