Télécharger la plaquette Cast3M


fr en

Télécharger ccgtsu.eso

Retour à la liste

Numérotation des lignes : 

ccgtsu

  1. C CCGTSU    SOURCE    GOUNAND   26/01/09    21:15:13     12441          
  2.       SUBROUTINE CCGTSU(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         : CCGTSU
  9. C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss :
  10. C               Tension de surface
  11. C
  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, 17/01/07, version initiale
  26. C HISTORIQUE : v1, 17/01/07, création
  27. C HISTORIQUE :
  28. C HISTORIQUE :
  29. C***********************************************************************
  30. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  31. C en cas de modification de ce sous-programme afin de faciliter
  32. C la maintenance !
  33. C***********************************************************************
  34.  
  35. -INC PPARAM
  36. -INC CCOPTIO
  37. -INC CCREEL
  38. -INC TNLIN
  39. *-INC SMCHAEL
  40.       INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
  41.       POINTEUR FC.MCHEVA
  42.       POINTEUR LCOF.LCHEVA
  43.       POINTEUR JMAJAC.MCHEVA
  44.       POINTEUR JMIJAC.MCHEVA
  45.       POINTEUR JDTJAC.MCHEVA
  46.       POINTEUR JMAREG.MCHEVA
  47.       POINTEUR JPC.MCHEVA
  48.       POINTEUR JTS.MCHEVA
  49.       CHARACTER*8 NOMLOI
  50.       INTEGER ICOF
  51. *
  52. -INC TMXMAT
  53.       POINTEUR JAC.MXMAT,JT.MXMAT
  54.       POINTEUR G.MXMAT
  55.       POINTEUR IG.MXMAT
  56.       POINTEUR JPT.MXMAT
  57.       POINTEUR JPTJT.MXMAT
  58. *
  59.       INTEGER LAXSP
  60.       REAL*8 DEUPI,QUATPI,XR
  61.       REAL*8 DETG(1)
  62.       PARAMETER (DEUPI=2.D0*XPI,QUATPI=4.D0*XPI)
  63. *
  64.       INTEGER IMPR,IRET
  65. *
  66. * Executable statements
  67. *
  68.       IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgtsu'
  69. C      IF (.NOT.(IDIM.EQ.1)) THEN
  70. C         WRITE(IOIMP,*) 'IDIM=',IDIM,' ?'
  71. C         GOTO 9999
  72. C      ENDIF
  73.       NLFC=FC.WELCHE(/6)
  74.       NPFC=FC.WELCHE(/5)
  75.       ICOF=0
  76. *
  77. *     Récupération du coefficent de tension de surface
  78. *
  79.       ICOF=ICOF+1
  80.       JTS =LCOF.LISCHE(ICOF)
  81.       NLTS=JTS.WELCHE(/6)
  82.       NPTS=JTS.WELCHE(/5)
  83. *
  84.       ICOF=ICOF+1
  85.       JMAJAC=LCOF.LISCHE(ICOF)
  86.       NLJA=JMAJAC.WELCHE(/6)
  87.       NPJA=JMAJAC.WELCHE(/5)
  88.       IREF=JMAJAC.WELCHE(/4)
  89.       IREL=JMAJAC.WELCHE(/3)
  90. *
  91.       IF (IREL.NE.IDIM) THEN
  92.          WRITE(IOIMP,*) 'Erreur dims JMAJAC'
  93.          GOTO 9999
  94.       ENDIF
  95. *
  96.       ICOF=ICOF+1
  97.       ICOF=ICOF+1
  98.       ICOF=ICOF+1
  99. *
  100. *!      WRITE(IOIMP,*) 'IFOMOD=',IFOMOD
  101.       IF (IFOMOD.EQ.0.OR.IFOMOD.EQ.4) THEN
  102. *!         WRITE(IOIMP,*) 'MODE AXISYMETRIQUE !!!!!!'
  103.          LAXSP=1
  104.       ELSEIF (IFOMOD.EQ.5) THEN
  105. *!         WRITE(IOIMP,*) 'MODE SPHERIQUE !!!!!!'
  106.          LAXSP=2
  107.       ELSE
  108. *!         WRITE(IOIMP,*) 'MODE PLAN !!!!!!'
  109.          LAXSP=0
  110.       ENDIF
  111. *!      WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  112. *!      WRITE(IOIMP,*) 'LAXSP=',LAXSP
  113. *
  114.       IF (LAXSP.GT.0) THEN
  115.          ICOF=ICOF+1
  116.          ICOF=ICOF+1
  117.          JPC=LCOF.LISCHE(ICOF)
  118.          NLPC=JPC.WELCHE(/6)
  119.          NPPC=JPC.WELCHE(/5)
  120.          IF (((NLPC.NE.1).AND.(NLPC.NE.NLFC)).OR.
  121.      $        ((NPPC.NE.1).AND.(NPPC.NE.NPFC))) THEN
  122.             WRITE(IOIMP,*) 'Erreur dims JPC'
  123.             GOTO 9999
  124.          ENDIF
  125.       ELSE
  126.          JPC=0
  127.          NLPC=0
  128.          NPPC=0
  129.       ENDIF
  130. *
  131. * Objets temporaires
  132. *
  133.       LDIM1=IREL
  134.       LDIM2=IREF
  135.       SEGINI,JAC
  136.       SEGINI,JPT
  137.       LDIM1=IREF
  138.       LDIM2=IREL
  139.       SEGINI,JT
  140.       LDIM1=IREF
  141.       LDIM2=IREF
  142.       SEGINI,G
  143.       SEGINI,IG
  144.       LDIM1=IREL
  145.       LDIM2=IREL
  146.       SEGINI,JPTJT
  147. *
  148.       DO ILFC=1,NLFC
  149.          IF (NLTS.EQ.1) THEN
  150.             ILTS=1
  151.          ELSE
  152.             ILTS=ILFC
  153.          ENDIF
  154.          IF (NLJA.EQ.1) THEN
  155.             ILJA=1
  156.          ELSE
  157.             ILJA=ILFC
  158.          ENDIF
  159.          IF (NLPC.EQ.1) THEN
  160.             ILPC=1
  161.          ELSE
  162.             ILPC=ILFC
  163.          ENDIF
  164.  
  165.  
  166.          DO IPFC=1,NPFC
  167.             IF (NPTS.EQ.1) THEN
  168.                IPTS=1
  169.             ELSE
  170.                IPTS=IPFC
  171.             ENDIF
  172.             IF (NPJA.EQ.1) THEN
  173.                IPJA=1
  174.             ELSE
  175.                IPJA=IPFC
  176.             ENDIF
  177.             IF (NPPC.EQ.1) THEN
  178.                IPPC=1
  179.             ELSE
  180.                IPPC=IPFC
  181.             ENDIF
  182. *
  183. * Valeur du coefficient
  184. *
  185.             IF (JTS.NE.0) THEN
  186.                XTS=JTS.WELCHE(1,1,1,1,IPTS,ILTS)
  187.             ENDIF
  188. *
  189. * Première coordonnée
  190. *
  191.             IF (JPC.NE.0) THEN
  192.                XR=JPC.WELCHE(1,1,1,1,IPPC,ILPC)
  193.             ENDIF
  194. *
  195. * Copie du jacobien
  196. *
  197.             CALL MAMA(JMAJAC.WELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF,
  198.      $           'COPIE   ',
  199.      $           JAC.XMAT,IREL,IREF,
  200.      $           IMPR,IRET)
  201.             IF (IRET.NE.0) GOTO 9999
  202. C            SEGPRT,JAC
  203. *
  204. * Calcul de la métrique G
  205. *
  206. *     Calcul de Jt
  207.             CALL MAMA(JAC.XMAT,IREL,IREF,
  208.      $           'TRANSPOS',JT.XMAT,IREF,IREL,
  209.      $           IMPR,IRET)
  210.             IF (IRET.NE.0) GOTO 9999
  211. *     Calcul de G=JtJ
  212.             CALL MAMAMA(JT.XMAT,IREF,IREL,JAC.XMAT,IREL,IREF,
  213.      $           'FOIS    ',G.XMAT,IREF,IREF,
  214.      $           IMPR,IRET)
  215.             IF (IRET.NE.0) GOTO 9999
  216. *
  217. * Calcul de l'inverse, du déterminant et trace de l'inverse
  218. *
  219.             CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,IMPR,IRET)
  220.             IF (IRET.NE.0) THEN
  221.               WRITE(IOIMP,*)
  222.      $             'JtJ non inversible'
  223.                GOTO 9999
  224.             ENDIF
  225.             SDETG=SQRT(ABS(DETG(1)))
  226. *     Calcul de JPT=J (JtJ)-1
  227.             CALL MAMAMA(JAC.XMAT,IREL,IREF,IG.XMAT,IREF,IREF,
  228.      $           'FOIS    ',JPT.XMAT,IREL,IREF,
  229.      $           IMPR,IRET)
  230.             IF (IRET.NE.0) GOTO 9999
  231. *     Calcul de JPT=J (JtJ)-1
  232.             CALL MAMAMA(JPT.XMAT,IREL,IREF,JT.XMAT,IREF,IREL,
  233.      $           'FOIS    ',JPTJT.XMAT,IREL,IREL,
  234.      $           IMPR,IRET)
  235.             IF (IRET.NE.0) GOTO 9999
  236. *
  237. * Calcul de la fonctionnelle
  238. *
  239.             IF (NOMLOI.EQ.'TSUF    ') THEN
  240.                CONTRI=XTS*SDETG
  241.                IF (LAXSP.EQ.1) THEN
  242.                   CONTRI=CONTRI*2.D0*XPI*XR
  243.                ELSEIF (LAXSP.EQ.2) THEN
  244.                   CONTRI=CONTRI*4.D0*XPI*XR*XR
  245.                ENDIF
  246. *
  247. * Calcul du résidu
  248. *
  249.             ELSEIF (NOMLOI(1:4).EQ.'TSUR') THEN
  250.                CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
  251.                IF (IRET.NE.0) GOTO 9999
  252.                CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
  253.                IF (IRET.NE.0) GOTO 9999
  254.                IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN
  255.                   IF (LAXSP.EQ.1) THEN
  256.                      CONTRI=XTS*2.D0*XPI*SDETG
  257.                   ELSEIF (LAXSP.EQ.2) THEN
  258.                      CONTRI=XTS*8.D0*XPI*XR*SDETG
  259.                   ELSE
  260.                      WRITE(IOIMP,*) 'LAXSP=',LAXSP
  261.                      GOTO 9999
  262.                   ENDIF
  263.                ELSE
  264.                   CONTRI=XTS*SDETG*JPT.XMAT(IDIM1,IDIM2)
  265.                   IF (LAXSP.EQ.1) THEN
  266.                      CONTRI=CONTRI*DEUPI*XR
  267.                   ELSEIF (LAXSP.EQ.2) THEN
  268.                      CONTRI=CONTRI*QUATPI*XR*XR
  269.                   ENDIF
  270.                ENDIF
  271. *
  272. * Calcul du jacobien
  273. *
  274.             ELSEIF (NOMLOI(1:4).EQ.'TSUJ') THEN
  275.                CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
  276.                IF (IRET.NE.0) GOTO 9999
  277.                CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
  278.                IF (IRET.NE.0) GOTO 9999
  279.                CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET)
  280.                IF (IRET.NE.0) GOTO 9999
  281.                CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET)
  282.                IF (IRET.NE.0) GOTO 9999
  283.                IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)
  284.      $              .AND.(IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN
  285.                   IF (LAXSP.EQ.2) THEN
  286.                      CONTRI=XTS*8.D0*XPI*SDETG
  287.                   ELSE
  288.                      WRITE(IOIMP,*) 'LAXSP=',LAXSP
  289.                      WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  290.                      GOTO 9999
  291.                   ENDIF
  292.                ELSEIF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN
  293.                   IF (LAXSP.GE.1) THEN
  294.                      CONTRI=XTS*DEUPI*SDETG*JPT.XMAT(IDIM3,IDIM4)
  295.                      IF (LAXSP.EQ.2) THEN
  296.                         CONTRI=CONTRI*4.D0*XR
  297.                      ENDIF
  298.                   ELSE
  299.                      WRITE(IOIMP,*) 'LAXSP=',LAXSP
  300.                      WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  301.                      GOTO 9999
  302.                   ENDIF
  303.                ELSEIF ((IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN
  304.                   IF (LAXSP.GE.1) THEN
  305.                      CONTRI=XTS*DEUPI*SDETG*JPT.XMAT(IDIM1,IDIM2)
  306.                      IF (LAXSP.EQ.2) THEN
  307.                         CONTRI=CONTRI*4.D0*XR
  308.                      ENDIF
  309.                   ELSE
  310.                      WRITE(IOIMP,*) 'LAXSP=',LAXSP
  311.                      WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  312.                      GOTO 9999
  313.                   ENDIF
  314.                ELSE
  315.                   CONTRI=JPT.XMAT(IDIM1,IDIM2)*JPT.XMAT(IDIM3,IDIM4)
  316.                   IF (IDIM1.EQ.IDIM3) THEN
  317.                      CONTRI=CONTRI+IG.XMAT(IDIM2,IDIM4)
  318.                   ENDIF
  319.                   CONTRI=CONTRI-
  320.      $                 (JPTJT.XMAT(IDIM1,IDIM3)*IG.XMAT(IDIM2,IDIM4))
  321.                   CONTRI=CONTRI-
  322.      $                 (JPT.XMAT(IDIM3,IDIM2)*JPT.XMAT(IDIM1,IDIM4))
  323.                   CONTRI=CONTRI*XTS*SDETG
  324.                   IF (LAXSP.EQ.1) THEN
  325.                      CONTRI=CONTRI*DEUPI*XR
  326.                   ELSEIF (LAXSP.EQ.2) THEN
  327.                      CONTRI=CONTRI*QUATPI*XR*XR
  328.                   ENDIF
  329.                ENDIF
  330. *
  331. * Calcul d'une partie du jacobien
  332. *
  333.             ELSEIF (NOMLOI(1:3).EQ.'TSU') THEN
  334.                CALL CH2INT(NOMLOI(4:4),IDIM0,IMPR,IRET)
  335.                IF (IRET.NE.0) GOTO 9999
  336.                CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
  337.                IF (IRET.NE.0) GOTO 9999
  338.                CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
  339.                IF (IRET.NE.0) GOTO 9999
  340.                CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET)
  341.                IF (IRET.NE.0) GOTO 9999
  342.                CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET)
  343.                IF (IRET.NE.0) GOTO 9999
  344.                IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)
  345.      $              .AND.(IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN
  346.                   WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  347.                   GOTO 9999
  348.                ELSEIF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN
  349.                   WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  350.                   GOTO 9999
  351.                ELSEIF ((IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN
  352.                   WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  353.                   GOTO 9999
  354.                ELSE
  355.                   IF (IDIM0.EQ.1) THEN
  356.                      CONTRI=JPT.XMAT(IDIM1,IDIM2)*JPT.XMAT(IDIM3,IDIM4)
  357.                   ELSEIF (IDIM0.EQ.2) THEN
  358.                      IF (IDIM1.EQ.IDIM3) THEN
  359.                         CONTRI=IG.XMAT(IDIM2,IDIM4)
  360.                      ELSE
  361.                         CONTRI=0.D0
  362.                      ENDIF
  363.                   ELSEIF (IDIM0.EQ.3) THEN
  364.                      CONTRI=
  365.      $                    -(JPTJT.XMAT(IDIM1,IDIM3)*IG.XMAT(IDIM2,IDIM4)
  366.      $                    )
  367.                   ELSEIF (IDIM0.EQ.4) THEN
  368.                      CONTRI=-(JPT.XMAT(IDIM3,IDIM2)*JPT.XMAT(IDIM1,IDIM4
  369.      $                    ))
  370.                   ELSE
  371.                      WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  372.                      GOTO 9999
  373.                   ENDIF
  374.                   CONTRI=CONTRI*XTS*SDETG
  375.                   IF (LAXSP.EQ.1) THEN
  376.                      CONTRI=CONTRI*DEUPI*XR
  377.                   ELSEIF (LAXSP.EQ.2) THEN
  378.                      CONTRI=CONTRI*QUATPI*XR*XR
  379.                   ENDIF
  380.                ENDIF
  381. *!               WRITE(IOIMP,*) 'CONTRI=',CONTRI
  382.             ELSE
  383.                WRITE(IOIMP,*) 'Erreur grave'
  384.                GOTO 9999
  385.             ENDIF
  386.             FC.WELCHE(1,1,1,1,IPFC,ILFC)=
  387.      $           FC.WELCHE(1,1,1,1,IPFC,ILFC)+
  388.      $           CONTRI
  389.          ENDDO
  390.       ENDDO
  391.       SEGSUP,JPTJT
  392.       SEGSUP,JPT
  393.       SEGSUP,IG
  394.       SEGSUP,G
  395.       SEGSUP,JT
  396.       SEGSUP,JAC
  397. *
  398. * Normal termination
  399. *
  400.       IRET=0
  401.       RETURN
  402. *
  403. * Format handling
  404. *
  405. *
  406. * Error handling
  407. *
  408.  9999 CONTINUE
  409.       IRET=1
  410.       WRITE(IOIMP,*) 'An error was detected in subroutine ccgtsu'
  411.       RETURN
  412. *
  413. * End of subroutine CCGTSU
  414. *
  415.       END
  416.  
  417.  

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