Télécharger la plaquette Cast3M


fr en

Télécharger testpg.eso

Retour à la liste

Numérotation des lignes : 

testpg

  1. C TESTPG    SOURCE    GOUNAND   23/07/31    21:15:06     11713          
  2.       SUBROUTINE TESTPG()
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM         : TESTPG
  7. C PROJET      : Noyau linéaire NLIN
  8. C DESCRIPTION : On vérifie les méthodes d'intégration numériques. On
  9. C               vérifie l'exactitude de l'intégration des monomes
  10. C               d'ordre inférieur ou égal à celui de la méthode.
  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 APPELES (E/S)    :
  18. C APPELES (BLAS)   :
  19. C APPELES (CALCUL) :
  20. C APPELE PAR       : PILOT
  21. C***********************************************************************
  22. C SYNTAXE GIBIANE    :
  23. C ENTREES            :
  24. C ENTREES/SORTIES    :
  25. C SORTIES            :
  26. C***********************************************************************
  27. C VERSION    : v1, 22/07/99, version initiale
  28. C HISTORIQUE : v1, 22/07/99, création
  29. C HISTORIQUE :
  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 SPOGAU
  42.       POINTEUR MYPGS.POGAUS
  43.       POINTEUR PGCOUR.POGAU
  44. *
  45. * Exposant des monomes
  46. *
  47.       SEGMENT XPMONO
  48.       INTEGER EXPOSA(NDESP)
  49.       ENDSEGMENT
  50.       INTEGER NDESP
  51.       POINTEUR MYMONO.XPMONO
  52. *
  53.       INTEGER IMPR,IRET
  54. *
  55.       INTEGER NBMETH,I,J,K
  56.       INTEGER NORDM,NESREF,NPOGAU
  57.       INTEGER IMETH,IORD,IEXP,IEXP1,IEXP2,IPG,ICOOR
  58.       REAL*8 XFACT
  59.       REAL*8 XPRMON,SOLEX,SOLAPP,SERR
  60. *
  61. * Executable statements
  62. *
  63. *
  64. * Initialisation du segment contenant les informations sur les
  65. * méthodes d'intégration (type Gauss).
  66. *
  67.       IMPR=0
  68.       CALL INPGS(MYPGS,IMPR,IRET)
  69.       IF (IRET.NE.0) GOTO 9999
  70. *
  71. * Puis on teste chacune de méthodes
  72. *
  73.       IMPR=0
  74.       SEGACT MYPGS
  75.       NBMETH=MYPGS.LISPG(/1)
  76.       DO 1 IMETH=1,NBMETH
  77.          PGCOUR=MYPGS.LISPG(IMETH)
  78.          SEGACT PGCOUR
  79.          WRITE(IOIMP,*) 'Methode d''integration ',IMETH,' : '
  80.      $        ,PGCOUR.NOMPG
  81.          IF (IMPR.GT.1) THEN
  82.             CALL PRPG(PGCOUR,IMPR,IRET)
  83.             IF (IRET.NE.0) GOTO 9999
  84.          ENDIF
  85.          NORDM=PGCOUR.NORDPG
  86.          NESREF=PGCOUR.XCOPG(/1)
  87.          NPOGAU=PGCOUR.XCOPG(/2)
  88.          XYZMAX=XZERO
  89.          DO IPG=1,NPOGAU
  90.             DO ICOOR=1,NESREF
  91.                XYZMAX=MAX(abs(PGCOUR.XCOPG(ICOOR,IPG)),xyzmax)
  92.             ENDDO
  93.          ENDDO
  94. * Le monôme 1 intègre à 2 ** d sur un element de refernce
  95.          xtprec=max(((XYZMAX*2)**NESREF)*xzprec*3,sqrt(xpetit))
  96.          write(ioimp,*) 'xtprec=',xtprec
  97. *
  98.          DO 2 IORD=0,NORDM+1
  99.             IF (PGCOUR.FORLPG.EQ.'SEGMENT') THEN
  100.                NDESP=NESREF
  101.                SEGINI MYMONO
  102.                MYMONO.EXPOSA(1)=IORD
  103.                I=MYMONO.EXPOSA(1)
  104.                IF (PGCOUR.TYPMPG.EQ.'GAUSS'.OR.
  105.      $              PGCOUR.TYPMPG.EQ.'NEWTON-COTES') THEN
  106.                   IF (MOD(I,2).EQ.0) THEN
  107.                      SOLEX=2.D0/(DBLE(I)+1.D0)
  108.                   ELSE
  109.                      SOLEX=0.D0
  110.                   ENDIF
  111.                ELSEIF (PGCOUR.TYPMPG.EQ.'GAUSS-JACOBI10') THEN
  112.                   SOLEX=(1.D0/(DBLE(I+1)))
  113.      $                 -(1.D0/(DBLE(I+2)))
  114.                ELSEIF (PGCOUR.TYPMPG.EQ.'GAUSS-JACOBI20') THEN
  115.                   SOLEX=(1.D0/(DBLE(I+1)))
  116.      $                 -(2.D0/(DBLE(I+2)))
  117.      $                 +(1.D0/(DBLE(I+3)))
  118.                ELSE
  119.                   WRITE(IOIMP,*) 'Type de méthode non reconnue'
  120.                   GOTO 9999
  121.                ENDIF
  122.                SOLAPP=0.D0
  123.                DO 22 IPG=1,NPOGAU
  124.                   XPRMON=1.D0
  125.                   DO 222 ICOOR=1,NDESP
  126.                      XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG)
  127.      $                    **MYMONO.EXPOSA(ICOOR))
  128.  222              CONTINUE
  129.                   SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON)
  130.  22            CONTINUE
  131.                SEGSUP MYMONO
  132.                SERR=ABS(SOLEX-SOLAPP)
  133.                IF (IMPR.GT.3) THEN
  134.                   WRITE(IOIMP,*) 'Monome : ksi**',I,' :'
  135.                   WRITE(IOIMP,*) '  Solution exacte : ',SOLEX
  136.                   WRITE(IOIMP,*) '  Solution appr.  : ',SOLAPP
  137.                   WRITE(IOIMP,*) '  Erreur relative : ',SERR
  138.                ENDIF
  139.                IF (SERR.GT.XTPREC.AND.IORD.LE.NORDM) THEN
  140.                   GOTO 9998
  141.                ENDIF
  142.             ELSEIF (PGCOUR.FORLPG.EQ.'TRIANGLE') THEN
  143.                DO 26 IEXP=IORD,0,-1
  144.                   NDESP=NESREF
  145.                   SEGINI MYMONO
  146.                   MYMONO.EXPOSA(1)=IEXP
  147.                   MYMONO.EXPOSA(2)=IORD-IEXP
  148.                   I=MYMONO.EXPOSA(1)
  149.                   J=MYMONO.EXPOSA(2)
  150.                   SOLEX=(DFACT(I)*DFACT(J))/DFACT(I+J+2)
  151.                   SOLAPP=0.D0
  152.                   DO 262 IPG=1,NPOGAU
  153.                      XPRMON=1.D0
  154.                      DO 2622 ICOOR=1,NDESP
  155.                         XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG)
  156.      $                       **MYMONO.EXPOSA(ICOOR))
  157.  2622                CONTINUE
  158.                      SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON)
  159.  262              CONTINUE
  160.                   SEGSUP MYMONO
  161.                   SERR=ABS(SOLEX-SOLAPP)
  162.                   IF (IMPR.GT.3) THEN
  163.                      WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J,' :'
  164.                      WRITE(IOIMP,*) '  Solution exacte : ',SOLEX
  165.                      WRITE(IOIMP,*) '  Solution appr.  : ',SOLAPP
  166.                      WRITE(IOIMP,*) '  Erreur relative : ',SERR
  167.                   ENDIF
  168.                   IF (SERR.GT.XTPREC.AND.IORD.LE.NORDM) THEN
  169.                      GOTO 9998
  170.                   ENDIF
  171.  26            CONTINUE
  172.             ELSEIF (PGCOUR.FORLPG.EQ.'CARRE') THEN
  173.                DO 24 IEXP=IORD,0,-1
  174.                   NDESP=NESREF
  175.                   SEGINI MYMONO
  176.                   MYMONO.EXPOSA(1)=IEXP
  177.                   MYMONO.EXPOSA(2)=IORD-IEXP
  178.                   I=MYMONO.EXPOSA(1)
  179.                   J=MYMONO.EXPOSA(2)
  180.                   IF (MOD(I,2).EQ.0.AND.MOD(J,2).EQ.0) THEN
  181.                      SOLEX=4.D0/((DBLE(I)+1.D0)*(DBLE(J)+1.D0))
  182.                   ELSE
  183.                      SOLEX=0.D0
  184.                   ENDIF
  185.                   SOLAPP=0.D0
  186.                   DO 242 IPG=1,NPOGAU
  187.                      XPRMON=1.D0
  188.                      DO 2422 ICOOR=1,NDESP
  189.                         XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG)
  190.      $                       **MYMONO.EXPOSA(ICOOR))
  191.  2422                CONTINUE
  192.                      SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON)
  193.  242              CONTINUE
  194.                   SEGSUP MYMONO
  195.                   SERR=ABS(SOLEX-SOLAPP)
  196.                   IF (IMPR.GT.3) THEN
  197.                      WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J,' :'
  198.                      WRITE(IOIMP,*) '  Solution exacte : ',SOLEX
  199.                      WRITE(IOIMP,*) '  Solution appr.  : ',SOLAPP
  200.                      WRITE(IOIMP,*) '  Erreur relative : ',SERR
  201.                   ENDIF
  202.                   IF (SERR.GT.XTPREC.AND.IORD.LE.NORDM) THEN
  203.                      GOTO 9998
  204.                   ENDIF
  205.  24            CONTINUE
  206.             ELSEIF (PGCOUR.FORLPG.EQ.'TETRAEDRE') THEN
  207.                DO 28 IEXP1=IORD,0,-1
  208.                   DO 282 IEXP2=IORD-IEXP1,0,-1
  209.                      NDESP=NESREF
  210.                      SEGINI MYMONO
  211.                      MYMONO.EXPOSA(1)=IEXP1
  212.                      MYMONO.EXPOSA(2)=IEXP2
  213.                      MYMONO.EXPOSA(3)=IORD-IEXP1-IEXP2
  214.                      I=MYMONO.EXPOSA(1)
  215.                      J=MYMONO.EXPOSA(2)
  216.                      K=MYMONO.EXPOSA(3)
  217.                      SOLEX=(DFACT(I)*DFACT(J)*DFACT(K))/DFACT(I+J+K+3)
  218.                      SOLAPP=0.D0
  219.                      DO 2822 IPG=1,NPOGAU
  220.                         XPRMON=1.D0
  221.                         DO 28222 ICOOR=1,NDESP
  222.                            XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG)
  223.      $                          **MYMONO.EXPOSA(ICOOR))
  224. 28222                   CONTINUE
  225.                         SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON)
  226.  2822                CONTINUE
  227.                      SEGSUP MYMONO
  228.                      SERR=ABS(SOLEX-SOLAPP)
  229.                      IF (IMPR.GT.3) THEN
  230.                         WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J
  231.      $                       ,' :',' dzeta**',K,' :'
  232.                         WRITE(IOIMP,*) '  Solution exacte : ',SOLEX
  233.                         WRITE(IOIMP,*) '  Solution appr.  : ',SOLAPP
  234.                         WRITE(IOIMP,*) '  Erreur relative : ',SERR
  235.                      ENDIF
  236.                      IF (SERR.GT.XTPREC.AND.IORD.LE.NORDM) THEN
  237.                         GOTO 9998
  238.                      ENDIF
  239.  282              CONTINUE
  240.  28            CONTINUE
  241.             ELSEIF (PGCOUR.FORLPG.EQ.'PYRAMIDE') THEN
  242.                NDESP=NESREF
  243.                SEGINI MYMONO
  244.                MYMONO.EXPOSA(1)=0
  245.                MYMONO.EXPOSA(2)=0
  246.                MYMONO.EXPOSA(3)=IORD
  247.                I=MYMONO.EXPOSA(1)
  248.                J=MYMONO.EXPOSA(2)
  249.                K=MYMONO.EXPOSA(3)
  250.                SOLEX=2.D0*((1.D0/(DBLE(IORD+3)))
  251.      $              +(-2.D0/(DBLE(IORD+2)))
  252.      $              +(1.D0/(DBLE(IORD+1))))
  253.                SOLAPP=0.D0
  254.                DO 2922 IPG=1,NPOGAU
  255.                   XPRMON=1.D0
  256.                   DO 29222 ICOOR=1,NDESP
  257.                      XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG)
  258.      $                    **MYMONO.EXPOSA(ICOOR))
  259. 29222             CONTINUE
  260.                   SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON)
  261.  2922          CONTINUE
  262.                SEGSUP MYMONO
  263.                SERR=ABS(SOLEX-SOLAPP)
  264.                IF (IMPR.GT.3) THEN
  265.                   WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J
  266.      $                 ,' :',' dzeta**',K,' :'
  267.                   WRITE(IOIMP,*) '  Solution exacte : ',SOLEX
  268.                   WRITE(IOIMP,*) '  Solution appr.  : ',SOLAPP
  269.                   WRITE(IOIMP,*) '  Erreur relative : ',SERR
  270.                ENDIF
  271. * Intégration non exacte
  272.                IF (SERR.GT.SQRT(XTPREC).AND.IORD.LE.NORDM) THEN
  273.                   GOTO 9998
  274.                ENDIF
  275.             ELSEIF (PGCOUR.FORLPG.EQ.'PRISME') THEN
  276.                DO 30 IEXP1=IORD,0,-1
  277.                   DO 302 IEXP2=IORD-IEXP1,0,-1
  278.                      NDESP=NESREF
  279.                      SEGINI MYMONO
  280.                      MYMONO.EXPOSA(1)=IEXP1
  281.                      MYMONO.EXPOSA(2)=IEXP2
  282.                      MYMONO.EXPOSA(3)=IORD-IEXP1-IEXP2
  283.                      I=MYMONO.EXPOSA(1)
  284.                      J=MYMONO.EXPOSA(2)
  285.                      K=MYMONO.EXPOSA(3)
  286.                      IF (MOD(K,2).EQ.0) THEN
  287.                         SOLEX=((DFACT(I)*DFACT(J))/DFACT(I+J+2))*
  288.      $                       (2.D0/(DBLE(K)+1.D0))
  289.                      ELSE
  290.                         SOLEX=0.D0
  291.                      ENDIF
  292.                      SOLAPP=0.D0
  293.                      DO 3022 IPG=1,NPOGAU
  294.                         XPRMON=1.D0
  295.                         DO 30222 ICOOR=1,NDESP
  296.                            XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG)
  297.      $                          **MYMONO.EXPOSA(ICOOR))
  298. 30222                   CONTINUE
  299.                         SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON)
  300.  3022                CONTINUE
  301.                      SEGSUP MYMONO
  302.                      SERR=ABS(SOLEX-SOLAPP)
  303.                      IF (IMPR.GT.3) THEN
  304.                         WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J
  305.      $                       ,' :',' dzeta**',K,' :'
  306.                         WRITE(IOIMP,*) '  Solution exacte : ',SOLEX
  307.                         WRITE(IOIMP,*) '  Solution appr.  : ',SOLAPP
  308.                         WRITE(IOIMP,*) '  Erreur relative : ',SERR
  309.                      ENDIF
  310.                      IF (SERR.GT.XTPREC.AND.IORD.LE.NORDM) THEN
  311.                         GOTO 9998
  312.                      ENDIF
  313.  302              CONTINUE
  314.  30            CONTINUE
  315.             ELSEIF (PGCOUR.FORLPG.EQ.'CUBE') THEN
  316.                DO 32 IEXP1=IORD,0,-1
  317.                   DO 322 IEXP2=IORD-IEXP1,0,-1
  318.                      NDESP=NESREF
  319.                      SEGINI MYMONO
  320.                      MYMONO.EXPOSA(1)=IEXP1
  321.                      MYMONO.EXPOSA(2)=IEXP2
  322.                      MYMONO.EXPOSA(3)=IORD-IEXP1-IEXP2
  323.                      I=MYMONO.EXPOSA(1)
  324.                      J=MYMONO.EXPOSA(2)
  325.                      K=MYMONO.EXPOSA(3)
  326.                      IF (MOD(I,2).EQ.0.AND.MOD(J,2).EQ.0
  327.      $                    .AND.MOD(K,2).EQ.0) THEN
  328.                         SOLEX=8.D0/((DBLE(I)+1.D0)*(DBLE(J)+1.D0)
  329.      $                       *(DBLE(K)+1.D0))
  330.                      ELSE
  331.                         SOLEX=0.D0
  332.                      ENDIF
  333.                      SOLAPP=0.D0
  334.                      DO 3222 IPG=1,NPOGAU
  335.                         XPRMON=1.D0
  336.                         DO 32222 ICOOR=1,NDESP
  337.                            XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG)
  338.      $                          **MYMONO.EXPOSA(ICOOR))
  339. 32222                   CONTINUE
  340.                         SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON)
  341.  3222                CONTINUE
  342.                      SEGSUP MYMONO
  343.                      SERR=ABS(SOLEX-SOLAPP)
  344.                      IF (IMPR.GT.3) THEN
  345.                         WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J
  346.      $                       ,' :',' dzeta**',K,' :'
  347.                         WRITE(IOIMP,*) '  Solution exacte : ',SOLEX
  348.                         WRITE(IOIMP,*) '  Solution appr.  : ',SOLAPP
  349.                         WRITE(IOIMP,*) '  Erreur relative : ',SERR
  350.                      ENDIF
  351.                      IF (SERR.GT.XTPREC.AND.IORD.LE.NORDM) THEN
  352.                         GOTO 9998
  353.                      ENDIF
  354.  322              CONTINUE
  355.  32            CONTINUE
  356.             ELSE
  357.                WRITE(IOIMP,*) PGCOUR.FORLPG,' n''est pas une forme'
  358.                WRITE(IOIMP,*) 'd''élément reconnue.'
  359.                GOTO 9999
  360.             ENDIF
  361.  2       CONTINUE
  362.          SEGDES PGCOUR
  363.          WRITE(IOIMP,*) 'Ok'
  364.  1    CONTINUE
  365.       SEGDES MYPGS
  366. *
  367. * Normal termination
  368. *
  369.       IRET=0
  370.       CALL ECRLOG(.TRUE.)
  371.       RETURN
  372. *
  373. * Format handling
  374. *
  375. *
  376. * Error handling
  377. *
  378.  9998 CONTINUE
  379.       WRITE(IOIMP,*) 'SERR=',SERR,' => Pb. ds la méthode !!'
  380.  9999 CONTINUE
  381.       IRET=1
  382.       WRITE(IOIMP,*) 'An error was detected in subroutine testpg'
  383.       CALL ECRLOG(.FALSE.)
  384.       RETURN
  385. *
  386. * End of subroutine testpg
  387. *
  388.       END
  389.  
  390.  

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