Télécharger la plaquette Cast3M


fr en

Télécharger ccgnor.eso

Retour à la liste

Numérotation des lignes : 

ccgnor

  1. C CCGNOR    SOURCE    GOUNAND   26/01/09    21:15:09     12441          
  2.       SUBROUTINE CCGNOR(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         : CCGNOR
  9. C DESCRIPTION : Calcul des composantes d'un vecteur normal
  10. C
  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, 09/03/07, version initiale
  25. C HISTORIQUE : v1, 09/03/07, création
  26. C HISTORIQUE :
  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.  
  34. -INC PPARAM
  35. -INC CCOPTIO
  36. -INC CCREEL
  37. -INC TNLIN
  38. *-INC SMCHAEL
  39.       INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
  40.       POINTEUR FC.MCHEVA
  41.       POINTEUR LCOF.LCHEVA
  42.       POINTEUR JMAJAC.MCHEVA
  43.       POINTEUR JMIJAC.MCHEVA
  44.       POINTEUR JDTJAC.MCHEVA
  45.       CHARACTER*8 NOMLOI
  46.       INTEGER ICOF
  47. *
  48. -INC TMXMAT
  49. * Objets temporaires
  50.       POINTEUR JAC.MXMAT,JT.MXMAT,JP.MXMAT
  51.       POINTEUR G.MXMAT,IG.MXMAT
  52. *
  53.       SEGMENT MIMAT2
  54.       INTEGER IMAT2(2,2)
  55.       ENDSEGMENT
  56.       POINTEUR EIJ.MIMAT2
  57.       SEGMENT MIMAT3
  58.       INTEGER IMAT3(3,3,3)
  59.       ENDSEGMENT
  60.       POINTEUR EIJK.MIMAT3
  61. *
  62.       LOGICAL LBID
  63.       INTEGER LAXSP
  64. *
  65.       INTEGER IMPR,IRET
  66.       REAL*8 DETG(1)
  67. *
  68. * Executable statements
  69. *
  70. *      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgnor'
  71. *      WRITE(IOIMP,*) 'Entrée dans ccgnor'
  72. *      WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  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. *
  81.       ICOF=ICOF+1
  82.       JMAJAC=LCOF.LISCHE(ICOF)
  83.       NLJA=JMAJAC.WELCHE(/6)
  84.       NPJA=JMAJAC.WELCHE(/5)
  85.       IREF=JMAJAC.WELCHE(/4)
  86.       IREL=JMAJAC.WELCHE(/3)
  87. *      SEGPRT,JMAJAC
  88. *
  89. *      WRITE(IOIMP,*) 'IREL=',IREL
  90. *      WRITE(IOIMP,*) 'IREF=',IREF
  91.       IF (IREL.NE.IDIM) THEN
  92.          WRITE(IOIMP,*) 'Erreur dims JMAJAC'
  93.          GOTO 9999
  94.       ENDIF
  95.       IF (IREL.NE.IREF+1) THEN
  96.          WRITE(IOIMP,*) 'Le maillage donne nest pas une surface'
  97.          GOTO 9999
  98.       ENDIF
  99.       IF ((IREL.NE.2).AND.(IREL.NE.3)) THEN
  100.          WRITE(IOIMP,*) 'Ne marche quen dimension despace 2 ou 3'
  101.          GOTO 9999
  102.       ENDIF
  103. *
  104. * Objets temporaires et à préconditionner
  105. *
  106.       LDIM1=IREL
  107.       LDIM2=IREF
  108.       SEGINI,JAC
  109.       LDIM1=IREF
  110.       LDIM2=IREL
  111.       SEGINI,JT
  112. *      SEGINI,JP
  113.       LDIM1=IREF
  114.       LDIM2=IREF
  115.       SEGINI,G
  116.       SEGINI,IG
  117. *
  118. * Initialisation des tenseurs de permutation
  119. *
  120.       SEGINI,EIJ
  121.       IMULT=1
  122.       ICPT=0
  123.       DO I=1,2
  124.          DO J=1,2
  125.             IF (I.NE.J) THEN
  126.                ICPT=ICPT+1
  127.                IF (ICPT.EQ.2) THEN
  128.                   ICPT=0
  129.                   IMULT=IMULT*(-1)
  130.                ENDIF
  131.                EIJ.IMAT2(I,J)=IMULT
  132.             ENDIF
  133.          ENDDO
  134.       ENDDO
  135.       SEGINI,EIJK
  136.       IMULT=1
  137.       ICPT=0
  138.       DO I=1,3
  139.          DO J=1,3
  140.             IF (I.NE.J) THEN
  141.                DO K=1,3
  142.                   IF ((K.NE.I).AND.(K.NE.J)) THEN
  143.                      ICPT=ICPT+1
  144.                      IF (ICPT.EQ.2) THEN
  145.                         ICPT=0
  146.                         IMULT=IMULT*(-1)
  147.                      ENDIF
  148.                      EIJK.IMAT3(I,J,K)=IMULT
  149.                   ENDIF
  150.                ENDDO
  151.             ENDIF
  152.          ENDDO
  153.       ENDDO
  154. *      SEGPRT,EIJ
  155. *      SEGPRT,EIJK
  156. *
  157.       DO ILFC=1,NLFC
  158.          IF (NLJA.EQ.1) THEN
  159.             ILJA=1
  160.          ELSE
  161.             ILJA=ILFC
  162.          ENDIF
  163.          DO IPFC=1,NPFC
  164.             IF (NPJA.EQ.1) THEN
  165.                IPJA=1
  166.             ELSE
  167.                IPJA=IPFC
  168.             ENDIF
  169. *
  170. * Copie du jacobien
  171. *
  172.             CALL MAMA(JMAJAC.WELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF,
  173.      $           'COPIE   ',
  174.      $           JAC.XMAT,IREL,IREF,
  175.      $           IMPR,IRET)
  176.             IF (IRET.NE.0) GOTO 9999
  177. *     Calcul de Jt
  178.             CALL MAMA(JAC.XMAT,IREL,IREF,
  179.      $           'TRANSPOS',JT.XMAT,IREF,IREL,
  180.      $           IMPR,IRET)
  181.             IF (IRET.NE.0) GOTO 9999
  182. *     Calcul de G=JtJ
  183.             CALL MAMAMA(JT.XMAT,IREF,IREL,JAC.XMAT,IREL,IREF,
  184.      $           'FOIS    ',G.XMAT,IREF,IREF,
  185.      $           IMPR,IRET)
  186.             IF (IRET.NE.0) GOTO 9999
  187. * Calcul de l'inverse, du déterminant et trace de l'inverse de g
  188.             CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,IMPR,IRET)
  189.             IF (IRET.NE.0) THEN
  190.               WRITE(IOIMP,*)
  191.      $             'JtJ non inversible'
  192.                GOTO 9999
  193.             ENDIF
  194.             XM=SQRT(ABS(DETG(1)))
  195. **
  196. ** Calcul de la pseudo-inverse J+ = g-1 Jt
  197. **
  198. *               CALL MAMAMA(IG.XMAT,IREF,IREF,JT.XMAT,IREF,IREL,
  199. *     $              'FOIS    ',JP.XMAT,IREF,IREL,
  200. *     $              IMPR,IRET)
  201. *               IF (IRET.NE.0) GOTO 9999
  202.  
  203.             IF (NOMLOI(1:4).EQ.'VNOR') THEN
  204.                CALL CH2INT(NOMLOI(5:5),I,IMPR,IRET)
  205.                IF (IRET.NE.0) GOTO 9999
  206. **
  207. ** Calcul de la pseudo-inverse J+ = g-1 Jt
  208. **
  209. *               CALL MAMAMA(IG.XMAT,IREF,IREF,JT.XMAT,IREF,IREL,
  210. *     $              'FOIS    ',JP.XMAT,IREF,IREL,
  211. *     $              IMPR,IRET)
  212. *               IF (IRET.NE.0) GOTO 9999
  213.  
  214. *
  215.                IF (IREL.EQ.2) THEN
  216.                   CONTRI=0.D0
  217.                   DO J=1,2
  218.                      CONTRI=CONTRI-(EIJ.IMAT2(I,J)*
  219.      $                    JMAJAC.WELCHE(1,1,J,1,IPJA,ILJA))
  220.                   ENDDO
  221.                   CONTRI=CONTRI/XM
  222.                ELSEIF (IREL.EQ.3) THEN
  223. C                  XNUM=XM
  224. C                  XDENO=0.D0
  225. C* C'est louche parce que II ne varie pas de 1 à IREL
  226. C                  DO II=1,IREF
  227. C                     DO IJ=1,IREF
  228. C*                        DO IN=1,IREF
  229. C*                           DO IO=1,IREF
  230. C                              XDENO=XDENO+
  231. C     $                             (EIJ.IMAT2(II,IJ)*
  232. C     $                       JAC.XMAT(II,1)*JAC.XMAT(IJ,2))**2
  233. CC     $                             (EIJ.IMAT2(II,IJ)*EIJ.IMAT2(IN,IO)*
  234. CC     $                             JAC.XMAT(II,1)*JAC.XMAT(IJ,2)*
  235. CC     $                             JAC.XMAT(IN,1)*JAC.XMAT(IO,2))
  236. CC                           ENDDO
  237. CC                        ENDDO
  238. C                     ENDDO
  239. C                  ENDDO
  240. C                  XAL=XNUM/XDENO
  241. C                  CONTRI=0.D0
  242. C                  DO J=1,IREL
  243. C                     DO K=1,IREL
  244. C                        CONTRI=CONTRI+
  245. C     $                       (EIJK.IMAT3(I,J,K)*
  246. C     $                       JAC.XMAT(J,1)*JAC.XMAT(K,2))
  247. C                     ENDDO
  248. C                  ENDDO
  249. C                  CONTRI=CONTRI*XAL
  250.                   CONTRI=0.D0
  251.                   DO IA=1,IREF
  252.                      DO IB=1,IREF
  253.                         DO J=1,IREL
  254.                            DO K=1,IREL
  255.                               CONTRI=CONTRI+(EIJ.IMAT2(IA,IB)*
  256.      $                             EIJK.IMAT3(I,J,K)*
  257.      $                             JMAJAC.WELCHE(1,1,J,IA,IPJA,ILJA)*
  258.      $                             JMAJAC.WELCHE(1,1,K,IB,IPJA,ILJA))
  259. C     $                             JP.XMAT(IA,J)*JP.XMAT(IB,K))
  260.                            ENDDO
  261.                         ENDDO
  262.                      ENDDO
  263.                   ENDDO
  264. *                  SEGPRT,JAC
  265. *                  SEGPRT,JP
  266. *                  WRITE(IOIMP,*) 'XM=',XM
  267.                   CONTRI=CONTRI/(2.D0*XM)
  268. *                  WRITE(IOIMP,*) 'CONTRI=',CONTRI
  269.                ELSE
  270.                   WRITE(IOIMP,*) 'Erreur grave IREL=',IREL
  271.                   GOTO 9999
  272.                ENDIF
  273.             ELSEIF (NOMLOI(1:4).EQ.'VNOJ') THEN
  274.                CALL CH2INT(NOMLOI(5:5),I,IMPR,IRET)
  275.                IF (IRET.NE.0) GOTO 9999
  276.                CALL CH2INT(NOMLOI(6:6),J,IMPR,IRET)
  277.                IF (IRET.NE.0) GOTO 9999
  278.                CALL CH2INT(NOMLOI(7:7),L,IMPR,IRET)
  279.                IF (IRET.NE.0) GOTO 9999
  280.                IF (IREL.EQ.2) THEN
  281.                   CONTRI=0.D0
  282.                   CONTRI=CONTRI-(EIJ.IMAT2(I,J)*
  283.      $                 JMAJAC.WELCHE(1,1,L,1,IPJA,ILJA))
  284.                   CONTRI=CONTRI/XM
  285.                ELSEIF (IREL.EQ.3) THEN
  286.                   CONTRI=0.D0
  287.                   DO IA=1,IREF
  288.                      DO IB=1,IREF
  289. *                        DO J=1,IREL
  290.                            DO K=1,IREL
  291.                               CONTRI=CONTRI+(EIJ.IMAT2(IA,IB)*
  292.      $                             EIJK.IMAT3(I,J,K)*
  293.      $                             JMAJAC.WELCHE(1,1,L,IA,IPJA,ILJA)*
  294.      $                             JMAJAC.WELCHE(1,1,K,IB,IPJA,ILJA))
  295.                               CONTRI=CONTRI+(EIJ.IMAT2(IA,IB)*
  296.      $                             EIJK.IMAT3(I,K,J)*
  297.      $                             JMAJAC.WELCHE(1,1,K,IA,IPJA,ILJA)*
  298.      $                             JMAJAC.WELCHE(1,1,L,IB,IPJA,ILJA))
  299.                            ENDDO
  300. *                        ENDDO
  301.                      ENDDO
  302.                   ENDDO
  303. *                  SEGPRT,JAC
  304. *                  SEGPRT,JP
  305. *                  WRITE(IOIMP,*) 'XM=',XM
  306.                   CONTRI=CONTRI/(2.D0*XM)
  307. *                  WRITE(IOIMP,*) 'CONTRI=',CONTRI
  308.                ELSE
  309.                   WRITE(IOIMP,*) 'Erreur grave IREL=',IREL
  310.                   GOTO 9999
  311.                ENDIF
  312.             ELSE
  313.                WRITE(IOIMP,*) 'Erreur grave NOMLOI=',NOMLOI
  314.                GOTO 9999
  315.             ENDIF
  316.             FC.WELCHE(1,1,1,1,IPFC,ILFC)=
  317.      $           FC.WELCHE(1,1,1,1,IPFC,ILFC)+
  318.      $           CONTRI
  319.          ENDDO
  320.       ENDDO
  321.       SEGSUP,EIJK
  322.       SEGSUP,EIJ
  323. *      SEGPRT,FC
  324.       SEGSUP,JAC
  325. *      SEGSUP,JP
  326.       SEGSUP,JT
  327.       SEGSUP,G
  328.       SEGSUP,IG
  329. *
  330. * Normal termination
  331. *
  332.       IRET=0
  333.       RETURN
  334. *
  335. * Format handling
  336. *
  337. *
  338. * Error handling
  339. *
  340.  9999 CONTINUE
  341.       IRET=1
  342.       WRITE(IOIMP,*) 'An error was detected in subroutine ccgnor'
  343.       RETURN
  344. *
  345. * End of subroutine CCGNOR
  346. *
  347.       END
  348.  
  349.  

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