Télécharger la plaquette Cast3M


fr en

Télécharger ccgmus.eso

Retour à la liste

Numérotation des lignes : 

ccgmus

  1. C CCGMUS    SOURCE    GOUNAND   21/06/02    21:15:15     11022          
  2.       SUBROUTINE CCGMUS(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         : CCGMUS
  9. C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss :
  10. C               Taille d'un élément suivant une direction donnée
  11. C
  12. C LANGAGE     : ESOPE
  13. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  14. C               mél : gounand@semt2.smts.cea.fr
  15. C***********************************************************************
  16. C APPELES          :
  17. C APPELE PAR       :
  18. C***********************************************************************
  19. C ENTREES            :
  20. C ENTREES/SORTIES    :
  21. C SORTIES            : -
  22. C TRAVAIL            :
  23. C***********************************************************************
  24. C VERSION    : v1, 30/08/06, version initiale
  25. C HISTORIQUE : v1, 30/08/06, création
  26. C HISTORIQUE : 09/06/2015 : SG ajout de tests de positivite pour eviter
  27. C              des pbs dans les SQRTs
  28. C              Ajout du cas ou l'espace de reference est de dimension
  29. C              strictement inferieure a l'espace reel
  30. C HISTORIQUE :
  31. C***********************************************************************
  32. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  33. C en cas de modification de ce sous-programme afin de faciliter
  34. C la maintenance !
  35. C***********************************************************************
  36.  
  37. -INC PPARAM
  38. -INC CCOPTIO
  39. -INC CCREEL
  40. -INC TNLIN
  41. *-INC SMCHAEL
  42.       INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
  43.       POINTEUR FC.MCHEVA
  44.       POINTEUR LCOF.LCHEVA
  45.       POINTEUR JMAJAC.MCHEVA
  46.       POINTEUR JMIJAC.MCHEVA
  47.       POINTEUR JDTJAC.MCHEVA
  48.       POINTEUR JMAREG.MCHEVA
  49.       POINTEUR JDIAMA.MCHEVA
  50.       POINTEUR JPC.MCHEVA
  51.       POINTEUR JVIT.MCHEVA
  52.       POINTEUR JRHO.MCHEVA
  53.       POINTEUR JMU.MCHEVA
  54.       POINTEUR JPEC.MCHEVA
  55.       CHARACTER*8 NOMLOI
  56.       INTEGER ICOF
  57. *
  58. -INC TMXMAT
  59.       POINTEUR A.MXMAT
  60.       POINTEUR AP.MXMAT
  61.       POINTEUR JMA.MXMAT
  62.       POINTEUR KJMA.MXMAT
  63.       POINTEUR JM1.MXMAT,J.MXMAT
  64.       POINTEUR K.MXMAT
  65. *
  66.       SEGMENT MVIT
  67.       POINTEUR MVCOMP(IDIM).MCHEVA
  68.       ENDSEGMENT
  69.       LOGICAL LCRIT,LFIRST
  70.       LOGICAL LRRHO,LRMU,LRPEC
  71. *
  72.       INTEGER IMPR,IRET
  73. *
  74. * Executable statements
  75. *
  76.       IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgmus'
  77.       XPETI=SQRT(XPETIT)
  78.       CALL CH2INT(NOMLOI(7:7),IMETH,IMPR,IRET)
  79.       IF (IRET.NE.0) GOTO 9999
  80.       CALL CH2INT(NOMLOI(8:8),IDIMI,IMPR,IRET)
  81.       IF (IRET.NE.0) GOTO 9999
  82.       NLFC=FC.WELCHE(/6)
  83.       NPFC=FC.WELCHE(/5)
  84.       ICOF=0
  85. *
  86. * Récupération de RHO et MU
  87. *
  88.       ICOF=ICOF+1
  89.       JRHO=LCOF.LISCHE(ICOF)
  90.       NLRO=JRHO.WELCHE(/6)
  91.       NPRO=JRHO.WELCHE(/5)
  92.       ICOF=ICOF+1
  93.       JMU=LCOF.LISCHE(ICOF)
  94.       NLMU=JMU.WELCHE(/6)
  95.       NPMU=JMU.WELCHE(/5)
  96. *
  97. *     Récupération des coefficients de la metrique
  98. *
  99.       ICOF1=ICOF+1
  100.       SEGINI MVIT
  101.       DO IIDIM=1,IDIM
  102.          ICOF=ICOF+1
  103.          JVIT=LCOF.LISCHE(ICOF)
  104.          IF (ICOF.EQ.ICOF1) THEN
  105.             NLJV=JVIT.WELCHE(/6)
  106.             NPJV=JVIT.WELCHE(/5)
  107.          ELSE
  108.             NLJV2=JVIT.WELCHE(/6)
  109.             NPJV2=JVIT.WELCHE(/5)
  110.             IF (NLJV2.NE.NLJV.OR.NPJV2.NE.NPJV) THEN
  111.                WRITE(IOIMP,*) 'Erreur grave dims JVIT'
  112.                GOTO 9999
  113.             ENDIF
  114.          ENDIF
  115.          MVIT.MVCOMP(IIDIM)=JVIT
  116.       ENDDO
  117. *
  118. * Récupération du Péclet critique
  119. *
  120.       ICOF=ICOF+1
  121.       JPEC=LCOF.LISCHE(ICOF)
  122.       NLPEC=JPEC.WELCHE(/6)
  123.       NPPEC=JPEC.WELCHE(/5)
  124. *
  125. * sg : 09/06/2015
  126. * Verification que rho, mu et pec sont positifs.
  127. * Meme si les champs sont positifs, il peut y avoir un souci
  128. * avec les quadratiques aux points de Gauss
  129. * 220 0
  130. *Attention: une des composantes du champ est négative ou nulle
  131. * 220 0
  132. *sa puissance est mise à zéro
  133.       LRRHO=.FALSE.
  134.       DO ILRO=1,NLRO
  135.          DO IPRO=1,NPRO
  136.             XRHO=JRHO.WELCHE(1,1,1,1,IPRO,ILRO)
  137.             IF (XRHO.LT.XZERO) THEN
  138.                IF (.NOT.LRRHO) THEN
  139.                   LRRHO=.TRUE.
  140.                   WRITE(IOIMP,*) 'ccgmus.eso : RHO'
  141.                   CALL ERREUR(220)
  142.                ENDIF
  143.             ENDIF
  144.          ENDDO
  145.       ENDDO
  146.       LRMU=.FALSE.
  147.       DO ILMU=1,NLMU
  148.          DO IPMU=1,NPMU
  149.             XMU=JMU.WELCHE(1,1,1,1,IPMU,ILMU)
  150.             IF (XMU.LT.XZERO) THEN
  151.                IF (.NOT.LRMU) THEN
  152.                   LRMU=.TRUE.
  153.                   WRITE(IOIMP,*) 'ccgmus.eso : MU'
  154.                   CALL ERREUR(220)
  155.                ENDIF
  156.             ENDIF
  157.          ENDDO
  158.       ENDDO
  159.       LRPEC=.FALSE.
  160.       DO ILPEC=1,NLPEC
  161.          DO IPPEC=1,NPPEC
  162.             XPEC=JPEC.WELCHE(1,1,1,1,IPPEC,ILPEC)
  163.             IF (XPEC.LT.XZERO) THEN
  164.                IF (.NOT.LRPEC) THEN
  165.                   LRPEC=.TRUE.
  166.                   WRITE(IOIMP,*) 'ccgmus.eso : PEC'
  167.                   CALL ERREUR(220)
  168.                ENDIF
  169.             ENDIF
  170.          ENDDO
  171.       ENDDO
  172. *dbg                  IF (CONTRI.LT.0.D0) THEN
  173. *dbg                     WRITE(IOIMP,*) 'XRHO=',XRHO,' XMU=',XMU
  174. *dbg     $                    ,' XV=',XV,' XH=',XH,' XPEC=',XPEC
  175. *dbg                     WRITE(IOIMP,*) 'NPRO=',NPRO,'ILRO=',ILRO
  176. *dbg                     do ii=1,npro
  177. *dbg                        XRHO=JRHO.WELCHE(1,1,1,1,Ii,ILRO)
  178. *dbg                        WRITE(IOIMP,*) 'ii XRHO=',ii,XRHO
  179. *dbg                     enddo
  180. *dbg                     GOTO 9999
  181. *dbg                  ENDIF
  182. *
  183.       ICOF=ICOF+1
  184.       JMAJAC=LCOF.LISCHE(ICOF)
  185. C      NLJA=JMAJAC.WELCHE(/6)
  186. C      NPJA=JMAJAC.WELCHE(/5)
  187.       IREF=JMAJAC.WELCHE(/4)
  188.       IREL=JMAJAC.WELCHE(/3)
  189. *dbg      WRITE(IOIMP,*) 'IREF=',IREF,' IREL=',IREL
  190. *dbg      IF (IREL.NE.IDIM.OR.IREF.NE.IDIM) THEN
  191.       IF (IREL.NE.IDIM) THEN
  192.          WRITE(IOIMP,*) 'Erreur dims JMAJAC'
  193.          GOTO 9999
  194.       ENDIF
  195. *
  196.       ICOF=ICOF+1
  197.       JMIJAC=LCOF.LISCHE(ICOF)
  198.       IF (JMIJAC.EQ.0) THEN
  199.          WRITE(IOIMP,*) 'Erreur JMIJAC=0'
  200.          GOTO 9999
  201.       ENDIF
  202.       NLJI=JMIJAC.WELCHE(/6)
  203.       NPJI=JMIJAC.WELCHE(/5)
  204.       IREL2=JMIJAC.WELCHE(/4)
  205.       IREF2=JMIJAC.WELCHE(/3)
  206. *
  207.       IF (IREL2.NE.IREL.OR.IREF2.NE.IREF) THEN
  208.          WRITE(IOIMP,*) 'Erreur dims JMIJAC'
  209.          GOTO 9999
  210.       ENDIF
  211. *
  212.       ICOF=ICOF+1
  213.       ICOF=ICOF+1
  214.       JMAREG=LCOF.LISCHE(ICOF)
  215.       NLJR=JMAREG.WELCHE(/6)
  216.       NPJR=JMAREG.WELCHE(/5)
  217.       I1  =JMAREG.WELCHE(/4)
  218.       I2  =JMAREG.WELCHE(/3)
  219.       IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF))
  220.      $     THEN
  221.          WRITE(IOIMP,*) 'Erreur dims JMAREG'
  222.          GOTO 9999
  223.       ENDIF
  224.       ICOF=ICOF+1
  225.       JDIAMA=LCOF.LISCHE(ICOF)
  226.       NLJD=JDIAMA.WELCHE(/6)
  227.       NPJD=JDIAMA.WELCHE(/5)
  228.       I1  =JDIAMA.WELCHE(/4)
  229.       I2  =JDIAMA.WELCHE(/3)
  230.       IF ((NLJD.NE.1).OR.(NPJD.NE.1).OR.(I1.NE.1).OR.(I2.NE.1))
  231.      $     THEN
  232.          WRITE(IOIMP,*) 'Erreur dims JDIAMA'
  233.          GOTO 9999
  234.       ENDIF
  235.       XDIAMA=JDIAMA.WELCHE(1,1,1,1,1,1)
  236. *dbg
  237. *dbg      WRITE(IOIMP,*) 'XDIAMA=',XDIAMA
  238. *
  239. * Objets temporaires
  240. *
  241.       LDIM1=IREL
  242.       LDIM2=1
  243.       SEGINI,A
  244.       SEGINI,AP
  245.       LDIM1=IREF
  246.       LDIM2=1
  247.       SEGINI,JMA
  248.       LDIM1=IREF
  249.       LDIM2=1
  250.       SEGINI,KJMA
  251.       LDIM1=IREL
  252.       LDIM2=IREF
  253.       SEGINI,J
  254.       LDIM1=IREF
  255.       LDIM2=IREL
  256.       SEGINI,JM1
  257.       LDIM1=IREF
  258.       LDIM2=IREF
  259.       SEGINI,K
  260. *
  261. * Calcul de la métrique des éléments réguliers
  262. *
  263.       CALL MAMA(JMAREG.WELCHE,IREF,IREF,
  264.      $     'COPIE   ',K.XMAT,IREF,IREF,
  265.      $     IMPR,IRET)
  266.       IF (IRET.NE.0) GOTO 9999
  267. *      SEGPRT,H
  268. *
  269.       DO ILFC=1,NLFC
  270.          IF (NLRO.EQ.1) THEN
  271.             ILRO=1
  272.          ELSE
  273.             ILRO=ILFC
  274.          ENDIF
  275.          IF (NLMU.EQ.1) THEN
  276.             ILMU=1
  277.          ELSE
  278.             ILMU=ILFC
  279.          ENDIF
  280.          IF (NLJV.EQ.1) THEN
  281.             ILJV=1
  282.          ELSE
  283.             ILJV=ILFC
  284.          ENDIF
  285.          IF (NLPEC.EQ.1) THEN
  286.             ILPEC=1
  287.          ELSE
  288.             ILPEC=ILFC
  289.          ENDIF
  290. C         IF (NLJA.EQ.1) THEN
  291. C            ILJA=1
  292. C         ELSE
  293. C            ILJA=ILFC
  294. C         ENDIF
  295.          IF (NLJI.EQ.1) THEN
  296.             ILJI=1
  297.          ELSE
  298.             ILJI=ILFC
  299.          ENDIF
  300. *
  301.          DO IPFC=1,NPFC
  302.             IF (NPRO.EQ.1) THEN
  303.                IPRO=1
  304.             ELSE
  305.                IPRO=IPFC
  306.             ENDIF
  307.             IF (NPMU.EQ.1) THEN
  308.                IPMU=1
  309.             ELSE
  310.                IPMU=IPFC
  311.             ENDIF
  312.             IF (NPJV.EQ.1) THEN
  313.                IPJV=1
  314.             ELSE
  315.                IPJV=IPFC
  316.             ENDIF
  317.             IF (NPPEC.EQ.1) THEN
  318.                IPPEC=1
  319.             ELSE
  320.                IPPEC=IPFC
  321.             ENDIF
  322. C            IF (NPJA.EQ.1) THEN
  323. C               IPJA=1
  324. C            ELSE
  325. C               IPJA=IPFC
  326. C            ENDIF
  327.             IF (NPJI.EQ.1) THEN
  328.                IPJI=1
  329.             ELSE
  330.                IPJI=IPFC
  331.             ENDIF
  332.             XRHO=JRHO.WELCHE(1,1,1,1,IPRO,ILRO)
  333.             XMU=JMU.WELCHE(1,1,1,1,IPMU,ILMU)
  334.             XPEC=JPEC.WELCHE(1,1,1,1,IPPEC,ILPEC)
  335.             IF (LRRHO) XRHO=MAX(XRHO,XZERO)
  336.             IF (LRMU)  XMU=MAX(XMU,XZERO)
  337.             IF (LRPEC) XPEC=MAX(XPEC,XZERO)
  338.  
  339. *dbg            LFIRST=(ILFC.EQ.1).AND.(IPFC.EQ.1)
  340.             LFIRST=.FALSE.
  341.             IF (LFIRST) WRITE(IOIMP,*) 'XRHO=',XRHO,' XMU=',XMU,' XPEC='
  342.      $           ,XPEC
  343. *
  344. * Copie de la vitesse
  345. *
  346.             DO IIDIM=1,IDIM
  347.                JVIT=MVIT.MVCOMP(IIDIM)
  348.                A.XMAT(IIDIM,1)=JVIT.WELCHE(1,1,1,1,IPJV,ILJV)
  349.             ENDDO
  350.             IF (LFIRST) WRITE(IOIMP,*) 'A ',(A.XMAT(ii,1),ii=1,IDIM)
  351. *            SEGPRT,A
  352. *
  353. * Copie de l'inverse du jacobien ref->reel
  354. *
  355.             CALL MAMA(JMIJAC.WELCHE(1,1,1,1,IPJI,ILJI),IREF,IREL,
  356.      $           'COPIE   ',
  357.      $           JM1.XMAT,IREF,IREL,
  358.      $           IMPR,IRET)
  359.             IF (IRET.NE.0) GOTO 9999
  360. C            SEGPRT,JM1
  361. *
  362. * Copie du jacobien ref->reel
  363. *
  364.             CALL MAMA(JMAJAC.WELCHE(1,1,1,1,IPJI,ILJI),IREL,IREF,
  365.      $           'COPIE   ',
  366.      $           J.XMAT,IREL,IREF,
  367.      $           IMPR,IRET)
  368.             IF (IRET.NE.0) GOTO 9999
  369. C            SEGPRT,J
  370. *
  371. * Vecteur A dans le repère de référence
  372. *
  373.             CALL MAMAMA(JM1.XMAT,IREF,IREL,A.XMAT,IREL,1,
  374.      $           'FOIS    ',JMA.XMAT,IREF,1,
  375.      $           IMPR,IRET)
  376.             IF (IRET.NE.0) GOTO 9999
  377. *
  378. * Vecteur A' dans le repère réel (projection sur le domaine considéré)
  379. *
  380.             CALL MAMAMA(J.XMAT,IREL,IREF,JMA.XMAT,IREF,1,
  381.      $           'FOIS    ',AP.XMAT,IREL,1,
  382.      $           IMPR,IRET)
  383.             IF (IRET.NE.0) GOTO 9999
  384.             IF (LFIRST) WRITE(IOIMP,*) 'AP ',(AP.XMAT(ii,1),ii=1,IDIM)
  385. *
  386. * Vecteur A dans le repère de l'élément régulier
  387. *
  388.             CALL MAMAMA(K.XMAT,IREF,IREF,JMA.XMAT,IREF,1,
  389.      $           'FOIS    ',KJMA.XMAT,IREF,1,
  390.      $           IMPR,IRET)
  391.             IF (IRET.NE.0) GOTO 9999
  392. *
  393. * Normes de A
  394. *
  395.             XAREG2=XZERO
  396.             DO IIREF=1,IREF
  397.                XAREG2=XAREG2+(KJMA.XMAT(IIREF,1)**2)
  398.             ENDDO
  399.             XA2=XZERO
  400.             DO IIDIM=1,IDIM
  401.                XA2=XA2+(A.XMAT(IIDIM,1)**2)
  402.             ENDDO
  403.             XAP2=XZERO
  404.             DO IIDIM=1,IDIM
  405.                XAP2=XAP2+(AP.XMAT(IIDIM,1)**2)
  406.             ENDDO
  407.             IF (LFIRST) WRITE(IOIMP,*) 'XAREG2=',XAREG2,' XA2=',XA2,
  408.      $           ' XAP2=',XAP2
  409. *
  410. *old            IF (XA2.LT.XPETI) THEN
  411.             IF (XAP2.LT.XPETI) THEN
  412.                CONTRI=XZERO
  413.             ELSE
  414. *old               XV=SQRT(XA2)
  415. *old* Pondération par un diamètre externe !
  416. *old*               XH=SQRT(XA2/XAREG2)
  417. *old               XH=SQRT(XA2/XAREG2)*XDIAMA
  418. * Pour que le décentrement marche sur une surface
  419.                XV=SQRT(XAP2)
  420.                XH=SQRT(XAP2/XAREG2)*XDIAMA
  421. *
  422.                IF (IMETH.EQ.1) THEN
  423. * Diffusion artificielle UPWIND
  424.                   CONTRI=((XRHO*XV*XH)/XPEC)
  425.                   CONTRI=SQRT(CONTRI)
  426.                   IF (IDIMI.GT.0) THEN
  427.                      CONTRI=CONTRI*
  428. *old     $                    (A.XMAT(IDIMI,1)/XV)
  429.      $                    (AP.XMAT(IDIMI,1)/XV)
  430.                   ENDIF
  431.                   IF (LFIRST) WRITE(IOIMP,*) 'CONTRI=',CONTRI
  432.                ELSEIF (IMETH.EQ.2) THEN
  433. * Diffusion artificielle comme Alberto (equiv critical approx)
  434. * Doit-on décentrer ?
  435.                   LCRIT=((XRHO*XV*XH).GT.(XPEC*XMU))
  436.                   IF (LCRIT) THEN
  437.                      CONTRI=((XRHO*XV*XH)/XPEC)-XMU
  438.                      CONTRI=SQRT(CONTRI)
  439.                      IF (IDIMI.GT.0) THEN
  440.                         CONTRI=CONTRI*
  441. *old     $                       (A.XMAT(IDIMI,1)/XV)
  442.      $                       (AP.XMAT(IDIMI,1)/XV)
  443.                      ENDIF
  444.                   ELSE
  445.                      CONTRI=XZERO
  446.                   ENDIF
  447.                ELSEIF (IMETH.EQ.3) THEN
  448. * Diffusion artificielle (equiv doubly asymptotic)
  449.                   LCRIT=((XRHO*XV*XH).GT.(XPEC*XMU*3.D0))
  450.                   IF (LCRIT) THEN
  451.                      CONTRI=((XRHO*XV*XH)/XPEC)
  452.                   ELSE
  453.                      CONTRI=(((XRHO*XV*XH)/XPEC)**2)/(3.D0*XMU)
  454.                   ENDIF
  455.                   CONTRI=SQRT(CONTRI)
  456.                   IF (IDIMI.GT.0) THEN
  457.                      CONTRI=CONTRI*
  458. *old     $                    (A.XMAT(IDIMI,1)/XV)
  459.      $                    (AP.XMAT(IDIMI,1)/XV)
  460.                   ENDIF
  461.                ELSE
  462.                   WRITE(IOIMP,*) 'IMETH=',IMETH,' unknown'
  463.                   GOTO 9999
  464.                ENDIF
  465.             ENDIF
  466.             FC.WELCHE(1,1,1,1,IPFC,ILFC)=
  467.      $           FC.WELCHE(1,1,1,1,IPFC,ILFC)+
  468.      $           CONTRI
  469.          ENDDO
  470.       ENDDO
  471.       SEGSUP,K
  472.       SEGSUP,J
  473.       SEGSUP,JM1
  474.       SEGSUP,KJMA
  475.       SEGSUP,JMA
  476.       SEGSUP,AP
  477.       SEGSUP,A
  478.       SEGSUP,MVIT
  479. *
  480. * Normal termination
  481. *
  482.       IRET=0
  483.       RETURN
  484. *
  485. * Format handling
  486. *
  487. *
  488. * Error handling
  489. *
  490.  9999 CONTINUE
  491.       IRET=1
  492.       WRITE(IOIMP,*) 'An error was detected in subroutine ccgmus'
  493.       RETURN
  494. *
  495. * End of subroutine CCGMUS
  496. *
  497.       END
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  
  504.  

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