Télécharger la plaquette Cast3M


fr en

Télécharger ccgqme.eso

Retour à la liste

Numérotation des lignes : 

ccgqme

  1. C CCGQME    SOURCE    GOUNAND   26/01/09    21:15:11     12441          
  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 PPARAM
  34. -INC CCOPTIO
  35. -INC CCREEL
  36. -INC TNLIN
  37. *-INC SMCHAEL
  38.       INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
  39.       POINTEUR FC.MCHEVA
  40.       POINTEUR LCOF.LCHEVA
  41.       POINTEUR JMAJAC.MCHEVA
  42.       POINTEUR JMIJAC.MCHEVA
  43.       POINTEUR JDTJAC.MCHEVA
  44.       POINTEUR JMAREG.MCHEVA
  45.       POINTEUR JMET.MCHEVA
  46.       POINTEUR JTHE.MCHEVA
  47.       POINTEUR JGAM.MCHEVA
  48.       CHARACTER*8 NOMLOI
  49.       INTEGER ICOF
  50. *
  51. -INC TMXMAT
  52. * Objets temporaires
  53.       POINTEUR JAC.MXMAT,JT.MXMAT
  54.       POINTEUR G.MXMAT,IH.MXMAT,H.MXMAT,HIG.MXMAT,GIH.MXMAT
  55.       POINTEUR ME.MXMAT,MJ.MXMAT
  56. *
  57.       SEGMENT MCOF
  58.       POINTEUR COEF(IDIM,IDIM).MCHEVA
  59.       ENDSEGMENT
  60.       POINTEUR MET.MCOF
  61. *
  62.       LOGICAL LBID
  63.       INTEGER LAXSP
  64.       REAL*8 DEUPI,XR
  65.       REAL*8 XL,XM
  66.       REAL*8 DETGIH(1,1),DETH(1,1)
  67. *
  68.       INTEGER IMPR,IRET
  69. *
  70. * Executable statements
  71. *
  72.       IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgqme'
  73. C      IF (.NOT.(IDIM.EQ.1)) THEN
  74. C         WRITE(IOIMP,*) 'IDIM=',IDIM,' ?'
  75. C         GOTO 9999
  76. C      ENDIF
  77.       NLFC=FC.WELCHE(/6)
  78.       NPFC=FC.WELCHE(/5)
  79.       ICOF=0
  80.       XPET=SQRT(XPETIT)
  81. *
  82. *     Récupération des coefficients de la metrique
  83. *
  84.       SEGINI MET
  85.       DO IIDIM=1,IDIM
  86.          ICOF=ICOF+1
  87.          JMET=LCOF.LISCHE(ICOF)
  88.          IF (ICOF.EQ.1) THEN
  89.             NLJM=JMET.WELCHE(/6)
  90.             NPJM=JMET.WELCHE(/5)
  91.          ELSE
  92.             NLJM2=JMET.WELCHE(/6)
  93.             NPJM2=JMET.WELCHE(/5)
  94.             IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN
  95.                WRITE(IOIMP,*) 'Erreur grave dims JMET'
  96.                GOTO 9999
  97.             ENDIF
  98.          ENDIF
  99.          MET.COEF(IIDIM,IIDIM)=JMET
  100.       ENDDO
  101.       DO IIDIM=1,IDIM
  102.          NJ=IDIM-IIDIM
  103.          IF (NJ.GE.1) THEN
  104.             DO JIDIM=IIDIM+1,IDIM
  105.                ICOF=ICOF+1
  106.                JMET=LCOF.LISCHE(ICOF)
  107.                NLJM2=JMET.WELCHE(/6)
  108.                NPJM2=JMET.WELCHE(/5)
  109.                IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN
  110.                   WRITE(IOIMP,*) 'Erreur grave dims JMET2'
  111.                   GOTO 9999
  112.                ENDIF
  113.                MET.COEF(IIDIM,JIDIM)=JMET
  114.             ENDDO
  115.          ENDIF
  116.       ENDDO
  117. *
  118.       ICOF=ICOF+1
  119.       JMAJAC=LCOF.LISCHE(ICOF)
  120.       NLJA=JMAJAC.WELCHE(/6)
  121.       NPJA=JMAJAC.WELCHE(/5)
  122.       IREF=JMAJAC.WELCHE(/4)
  123.       IREL=JMAJAC.WELCHE(/3)
  124. *
  125.       IF (IREL.NE.IDIM) THEN
  126.          WRITE(IOIMP,*) 'Erreur dims JMAJAC'
  127.          GOTO 9999
  128.       ENDIF
  129. *
  130.       ICOF=ICOF+1
  131.       ICOF=ICOF+1
  132.       ICOF=ICOF+1
  133.       JMAREG=LCOF.LISCHE(ICOF)
  134.       NLJR=JMAREG.WELCHE(/6)
  135.       NPJR=JMAREG.WELCHE(/5)
  136.       I1  =JMAREG.WELCHE(/4)
  137.       I2  =JMAREG.WELCHE(/3)
  138.       IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF))
  139.      $     THEN
  140.          WRITE(IOIMP,*) 'Erreur dims JMAREG'
  141.          GOTO 9999
  142.       ENDIF
  143. *
  144. * Objets temporaires et à préconditionner
  145. *
  146.       LDIM1=IREL
  147.       LDIM2=IREF
  148.       SEGINI,JAC
  149.       SEGINI,MJ
  150.       LDIM1=IREF
  151.       LDIM2=IREL
  152.       SEGINI,JT
  153.       LDIM1=IREF
  154.       LDIM2=IREF
  155.       SEGINI,G
  156.       SEGINI,IH
  157.       SEGINI,H
  158.       SEGINI,HIG
  159.       SEGINI,GIH
  160.       LDIM1=IREL
  161.       LDIM2=IREL
  162.       SEGINI,ME
  163. *
  164. * Calcul de la métrique des éléments réguliers
  165. *
  166.       CALL MAMA(JMAREG.WELCHE,IREF,IREF,
  167.      $     'JTJ     ',H.XMAT,IREF,IREF,
  168.      $     IMPR,IRET)
  169.       IF (IRET.NE.0) GOTO 9999
  170. *      SEGPRT,H
  171. * Calcul de l'inverse, du déterminant et trace de l'inverse de h
  172.       CALL GEOLI2(IREF,1,1,H.XMAT,IH.XMAT,DETH,
  173.      $     IMPR,IRET)
  174.       IF (IRET.NE.0) THEN
  175.          WRITE(IOIMP,*)
  176.      $        'Metrique des elements reguliers non inversible'
  177.          SEGPRT,JMAREG
  178.          SEGPRT,H
  179.          GOTO 9999
  180.       ENDIF
  181. *      SEGPRT,IH
  182. *
  183.       DO ILFC=1,NLFC
  184.          IF (NLJM.EQ.1) THEN
  185.             ILJM=1
  186.          ELSE
  187.             ILJM=ILFC
  188.          ENDIF
  189.          IF (NLJA.EQ.1) THEN
  190.             ILJA=1
  191.          ELSE
  192.             ILJA=ILFC
  193.          ENDIF
  194.          DO IPFC=1,NPFC
  195.             IF (NPJM.EQ.1) THEN
  196.                IPJM=1
  197.             ELSE
  198.                IPJM=IPFC
  199.             ENDIF
  200.             IF (NPJA.EQ.1) THEN
  201.                IPJA=1
  202.             ELSE
  203.                IPJA=IPFC
  204.             ENDIF
  205. *
  206. * Copie des coefficients de la métrique
  207. *
  208.             DO IIDIM=1,IDIM
  209.                JMET=MET.COEF(IIDIM,IIDIM)
  210.                ME.XMAT(IIDIM,IIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM)
  211.             ENDDO
  212.             DO IIDIM=1,IIDIM
  213.                NJ=IDIM-IIDIM
  214.                IF (NJ.GE.1) THEN
  215.                   DO JIDIM=IIDIM+1,IDIM
  216.                      JMET=MET.COEF(IIDIM,JIDIM)
  217.                      ME.XMAT(IIDIM,JIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM)
  218.                      ME.XMAT(JIDIM,IIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM)
  219.                   ENDDO
  220.                ENDIF
  221.             ENDDO
  222. *
  223. * Copie du jacobien
  224. *
  225.             CALL MAMA(JMAJAC.WELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF,
  226.      $           'COPIE   ',
  227.      $           JAC.XMAT,IREL,IREF,
  228.      $           IMPR,IRET)
  229.             IF (IRET.NE.0) GOTO 9999
  230. *            SEGPRT,JAC
  231. *
  232. * Calcul de la métrique G
  233. *
  234. *     Calcul de Jt
  235.             CALL MAMA(JAC.XMAT,IREL,IREF,
  236.      $           'TRANSPOS',JT.XMAT,IREF,IREL,
  237.      $           IMPR,IRET)
  238.             IF (IRET.NE.0) GOTO 9999
  239. *     Calcul de MJ
  240.             CALL MAMAMA(ME.XMAT,IREL,IREL,JAC.XMAT,IREL,IREF,
  241.      $           'FOIS    ',MJ.XMAT,IREL,IREF,
  242.      $           IMPR,IRET)
  243.             IF (IRET.NE.0) GOTO 9999
  244. *     Calcul de G=JtMJ
  245.             CALL MAMAMA(JT.XMAT,IREF,IREL,MJ.XMAT,IREL,IREF,
  246.      $           'FOIS    ',G.XMAT,IREF,IREF,
  247.      $           IMPR,IRET)
  248.             IF (IRET.NE.0) GOTO 9999
  249. *     Calcul de gh-1, de sa trace et de son déterminant, qui peut être
  250. *     nul. La trace ne l'est pas, sauf si tous les points de l'element
  251. *     se confondent
  252.  
  253.             CALL MAMAMA(G.XMAT,IREF,IREF,IH.XMAT,IREF,IREF,
  254.      $           'FOIS    ',GIH.XMAT,IREF,IREF,IMPR,IRET)
  255.             IF (IRET.NE.0) GOTO 9999
  256.             CALL GEOLI2(IREF,1,1,GIH.XMAT,HIG.XMAT,DETGIH,
  257.      $           IMPR,IRET)
  258.             IF (IRET.GT.0) GOTO 9999
  259.             XM=SQRT(MAX(DETGIH(1,1),XZERO))
  260.             CALL MARE(GIH.XMAT,IREF,IREF,'TRACE   ',
  261.      $           XL,IMPR,IRET)
  262.             IF (IRET.NE.0) GOTO 9999
  263.             IF (ABS(XL).LE.XPET) THEN
  264.                WRITE(IOIMP,*) 'Nil jacobian matrix, check your mesh'
  265.                GOTO 9999
  266.             ENDIF
  267. *
  268. * Calcul des qualités de maillage
  269. *
  270.             IF (NOMLOI.EQ.'QEQU    ') THEN
  271.                CONTRI=XM
  272.             ELSEIF (NOMLOI.EQ.'QALI') THEN
  273.                IF (IREF.EQ.1) THEN
  274.                   CONTRI=1.D0
  275.                ELSE
  276.                   XIREF=DBLE(IREF)
  277.                   XNUM=XL
  278.                   XDEN=XIREF*(XM**(2.D0/XIREF))
  279.                   XEXP=XIREF/(2.D0*(XIREF-1.D0))
  280.                   CONTRI=(XDEN/XNUM)**XEXP
  281. *                  CONTRI=(XNUM/XDEN)**XEXP
  282. C                  SEGPRT,GIH
  283. C                  WRITE(IOIMP,*) 'TRGIH =',XL
  284. C                  WRITE(IOIMP,*) 'DETGIH=',(1.D0/DETHIG)
  285. C                  WRITE(IOIMP,*) 'XNUM=',XNUM
  286. C                  WRITE(IOIMP,*) 'XDEN=',XDEN
  287.                ENDIF
  288.             ELSE
  289.                WRITE(IOIMP,*) 'Erreur grave'
  290.                GOTO 9999
  291.             ENDIF
  292.             FC.WELCHE(1,1,1,1,IPFC,ILFC)=
  293.      $           FC.WELCHE(1,1,1,1,IPFC,ILFC)+
  294.      $           CONTRI
  295.          ENDDO
  296.       ENDDO
  297.       SEGSUP,ME
  298.       SEGSUP,GIH
  299.       SEGSUP,HIG
  300.       SEGSUP,H
  301.       SEGSUP,IH
  302.       SEGSUP,G
  303.       SEGSUP,JT
  304.       SEGSUP,MJ
  305.       SEGSUP,JAC
  306.       SEGSUP,MET
  307. *
  308. * Normal termination
  309. *
  310.          IRET=0
  311.          RETURN
  312. *
  313. * Format handling
  314. *
  315. *
  316. * Error handling
  317. *
  318.  9999    CONTINUE
  319.          IRET=1
  320.          WRITE(IOIMP,*) 'An error was detected in subroutine ccgqme'
  321.          RETURN
  322. *
  323. * End of subroutine CCGQME
  324. *
  325.          END
  326.  
  327.  

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