Télécharger la plaquette Cast3M


fr en

Télécharger tens1.eso

Retour à la liste

Numérotation des lignes : 

tens1

  1. C TENS1     SOURCE    GOUNAND   24/09/27    21:15:21     12020          
  2.       SUBROUTINE TENS1(ICHA,TYCHA,MLMOTS,IOTENS,ICHA1)
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM         : TENS1
  7. C DESCRIPTION : Opérations sur des tenseurs (unaires pour l'instant)
  8. C
  9. C
  10. C
  11. C LANGAGE     : ESOPE
  12. C AUTEUR      : Stephane GOUNAND (CEA/DES/ISAS/DM2S/SEMT/LTA)
  13. C               mel : gounand@semt2.smts.cea.fr
  14. C***********************************************************************
  15. C APPELES          :
  16. C APPELES (E/S)    :
  17. C APPELES (BLAS)   :
  18. C APPELES (CALCUL) :
  19. C APPELE PAR       :
  20. C***********************************************************************
  21. C SYNTAXE GIBIANE    :
  22. C ENTREES            :
  23. C ENTREES/SORTIES    :
  24. C SORTIES            :
  25. C***********************************************************************
  26. C VERSION    : v1, 28/08/2024, version initiale
  27. C HISTORIQUE : v1, 28/08/2024, creation
  28. C HISTORIQUE :
  29. C HISTORIQUE :
  30. C***********************************************************************
  31. -INC PPARAM
  32. -INC CCOPTIO
  33. -INC CCREEL
  34. -INC SMCHPOI
  35. -INC SMCHAML
  36. -INC SMLMOTS
  37. -INC TMTRAV
  38.       SEGMENT MTRAV1.MTRAV
  39.       SEGMENT MTRCML
  40.       REAL*8 CC(NNIN,N1PTEL,N1EL)
  41.       ENDSEGMENT
  42.  
  43.       CHARACTER*8 TYCHA,CCCOMP
  44. *
  45. C   ENTREEES
  46. C        A(3,3) = MATRICE A DIAGONALISER
  47. C   SORTIES
  48. C        D(3)   =  VALEURS PROPRES D(1) > D(2) ET D(3)=A(3,3)
  49. C        R(3,3) =  RESULTAT DES OPERATIONS TENSORIELLES
  50.       DIMENSION A(3,3),D(3),R(3,3)
  51.       DIMENSION IDXSYM(3,3,3)
  52.       DIMENSION IDXGEN(3,3,3)
  53.       DIMENSION INVSYM(2,6,3)
  54.       DIMENSION INVGEN(2,9,3)
  55.       DATA ((A(I,J),I=1,3),J=1,3) /9*0.D0/
  56.       DATA (D(I),I=1,3) /3*0.D0/
  57.       DATA ((R(I,J),I=1,3),J=1,3) /9*0.D0/
  58.       DATA (((IDXSYM(I,J,K),I=1,1),J=1,1),K=1,1) /1/
  59.       DATA (((IDXSYM(I,J,K),I=1,2),J=1,2),K=2,2) /1,2,2,3/
  60.       DATA (((IDXSYM(I,J,K),I=1,3),J=1,3),K=3,3) /1,2,4,2,3,5,4,5,6/
  61.       DATA (((IDXGEN(I,J,K),I=1,1),J=1,1),K=1,1) /1/
  62.       DATA (((IDXGEN(I,J,K),I=1,2),J=1,2),K=2,2) /1,2,3,4/
  63.       DATA (((IDXGEN(I,J,K),I=1,3),J=1,3),K=3,3) /1,2,3,4,5,6,7,8,9/
  64. *
  65.       DATA (((INVSYM(I,J,K),I=1,2),J=1,1),K=1,1) /1,1/
  66.       DATA (((INVSYM(I,J,K),I=1,2),J=1,3),K=2,2) /1,1,2,1,2,2/
  67.       DATA (((INVSYM(I,J,K),I=1,2),J=1,6),K=3,3)
  68.      $     /1,1,2,1,2,2,3,1,3,2,3,3/
  69.       DATA (((INVGEN(I,J,K),I=1,2),J=1,1),K=1,1) /1,1/
  70.       DATA (((INVGEN(I,J,K),I=1,2),J=1,4),K=2,2) /1,1,1,2,2,1,2,2/
  71.       DATA (((INVGEN(I,J,K),I=1,2),J=1,9),K=3,3)
  72.      $     /1,1,1,2,1,3,2,1,2,2,2,3,3,1,3,2,3,3/
  73. *
  74. *     DATA MOTENS/'NORM2','NORMINF','DET','TRACE','INVE','IDEN','LOG'
  75. *     ,'EXP','INVS','ABS','SPECTRAL'/
  76.  
  77. *
  78. * Executable statements
  79. *
  80.       SEGACT MLMOTS
  81.       NCOMP=MOTS(/2)
  82. * Cas entrée scalaire ou vecteur, symétrique, ordre 2, général
  83.       IF (NCOMP.EQ.1) THEN
  84.          IKAS=1
  85.       ELSEIF (NCOMP.EQ.IDIM) THEN
  86.          IKAS=2
  87.       ELSEIF (NCOMP*2.EQ.IDIM*(IDIM+1)) THEN
  88.          IKAS=3
  89.       ELSEIF (NCOMP.EQ.IDIM*IDIM) THEN
  90.          IKAS=4
  91.       ELSE
  92.          GOTO 9997
  93.       ENDIF
  94. * Nombre de composantes du résultat
  95.       IF (IOTENS.GE.1.AND.IOTENS.LE.4) THEN
  96.          NCOM1=1
  97.       ELSEIF (IOTENS.EQ.11) THEN
  98.          NCOM1=IDIM*(IDIM+1)
  99.       ELSE
  100.          NCOM1=NCOMP
  101.       ENDIF
  102.       IF (TYCHA.EQ.'CHPOINT ') THEN
  103.          MCHPOI=ICHA
  104. * Transformation du CHPOINT entré en objet de travail
  105.          CALL CP2TR2(MLMOTS,0,MCHPOI,MTRAV)
  106.          IF (IERR.NE.0) RETURN
  107.          SEGACT MTRAV*MOD
  108.          NNNOE=IGEO(/1)
  109.          IF (NCOM1.NE.NCOMP) THEN
  110.             NNIN=NCOM1
  111.             SEGINI MTRAV1
  112.             IF (NCOM1.EQ.1) THEN
  113.                MTRAV1.INCO(1)='SCAL'
  114.                MTRAV1.NHAR(1)=NHAR(1)
  115.             ELSEIF (IOTENS.EQ.11) THEN
  116.                ICMP=0
  117.                CCCOMP='SI'
  118.                DO i=1,IDIM
  119.                   WRITE(CCCOMP(3:3),FMT='(I1)') I
  120.                   WRITE(CCCOMP(4:4),FMT='(I1)') I
  121.                   ICMP=ICMP+1
  122.                   MTRAV1.INCO(ICMP)=CCCOMP
  123.                   MTRAV1.NHAR(ICMP)=NHAR(1)
  124.                ENDDO
  125.                CCCOMP='CO'
  126.                DO i=1,IDIM
  127.                   WRITE(CCCOMP(4:4),FMT='(I1)') I
  128.                   DO j=1,IDIM
  129.                      WRITE(CCCOMP(3:3),FMT='(I1)') J
  130.                      ICMP=ICMP+1
  131.                      MTRAV1.INCO(ICMP)=CCCOMP
  132.                      MTRAV1.NHAR(ICMP)=NHAR(1)
  133.                   ENDDO
  134.                ENDDO
  135.             ELSE
  136.                write(ioimp,*) 'CHPO COMP ncom1=',ncom1,' ????'
  137.                goto 9999
  138.             ENDIF
  139.             DO INNOE=1,NNNOE
  140.                MTRAV1.IGEO(INNOE)=IGEO(INNOE)
  141.             ENDDO
  142.          ENDIF
  143.       ELSEIF (TYCHA.EQ.'MCHAML') THEN
  144.          MCHELM=ICHA
  145.          CALL ACTOBJ('MCHAML  ',MCHELM,1)
  146.          NOID=1
  147. * Extrayons les composantes intéressantes
  148.          CALL EXCOC2(MCHELM,MLMOTS,MCHEL1,MLMOTS,NOID)
  149.          IF (IERR.NE.0) RETURN
  150.          SEGACT,MCHEL1*MOD
  151.       ELSE
  152.          GOTO 9998
  153.       ENDIF
  154. *      segprt,mtrav
  155. *
  156. *     write(ioimp,*) 'IOTENS=',IOTENS
  157. *
  158. * Un peu de gestion d'erreur avant la boucle...
  159. *
  160. * Seuls NORM2 et NORMINF disponibles pour les vecteurs
  161.       IF (IKAS.EQ.2.AND.IOTENS.GT.2) THEN
  162.          CALL ERREUR(803)
  163.          RETURN
  164.       ENDIF
  165.       MOTERR(1:8)=TYCHA
  166. *
  167. * La boucle sur les noeuds ou les éléments : on pourra paralléliser ici
  168. *
  169.       IF (TYCHA.EQ.'CHPOINT ') THEN
  170.          DO INNOE=1,NNNOE
  171. * Remplissage des tableaux de travail
  172.             IF (IKAS.EQ.1.OR.IKAS.EQ.2) THEN
  173.                DO I=1,NCOMP
  174.                   A(I,1)=BB(I,INNOE)
  175.                ENDDO
  176.             ELSEIF (IKAS.EQ.3) THEN
  177.                DO J=1,IDIM
  178.                   DO I=1,IDIM
  179.                      A(I,J)=BB(IDXSYM(I,J,IDIM),INNOE)
  180.                   ENDDO
  181.                ENDDO
  182.             ELSEIF (IKAS.EQ.4) THEN
  183.                DO J=1,IDIM
  184.                   DO I=1,IDIM
  185.                      A(I,J)=BB(IDXGEN(I,J,IDIM),INNOE)
  186.                   ENDDO
  187.                ENDDO
  188.             ELSE
  189.                write(ioimp,*) 'ikas=',ikas,' ????'
  190.                MOTERR(1:8)='TENS1'
  191.                call erreur(1039)
  192.                return
  193.             endif
  194. *     On prépare le message d'erreur 49 en cas de déterminant nul pour
  195. *     l'inversion d'une matrice...
  196.             IF (IOTENS.EQ.5) INTERR(1)=IGEO(INNOE)
  197. *
  198. * Faire les opérations tensorielles
  199. *
  200.             CALL TENS2(IOTENS,IKAS,A,D,R)
  201.             IF (IERR.NE.0) RETURN
  202. *
  203.             IF (NCOM1.NE.NCOMP) THEN
  204.                IF(NCOM1.EQ.1) THEN
  205.                   MTRAV1.BB(1,INNOE)=R(1,1)
  206.                   MTRAV1.IBIN(1,INNOE)=1
  207.                ELSEIF (IOTENS.EQ.11) THEN
  208.                   ICMP=0
  209.                   DO i=1,IDIM
  210.                      ICMP=ICMP+1
  211.                      MTRAV1.BB(ICMP,INNOE)=D(i)
  212.                      MTRAV1.IBIN(ICMP,INNOE)=1
  213.                   ENDDO
  214.                   DO i=1,IDIM
  215.                      DO j=1,IDIM
  216.                         ICMP=ICMP+1
  217.                         MTRAV1.BB(ICMP,INNOE)=R(j,i)
  218.                         MTRAV1.IBIN(ICMP,INNOE)=1
  219.                      ENDDO
  220.                   ENDDO
  221.                ELSE
  222.                   write(ioimp,*) 'CHPO VAL ncom1=',ncom1,' ????'
  223.                   goto 9999
  224.                ENDIF
  225.             ELSE
  226.                IF (IKAS.EQ.1.OR.IKAS.EQ.2) THEN
  227.                   DO I=1,NCOMP
  228.                      BB(I,INNOE)=R(I,1)
  229.                      IBIN(I,INNOE)=1
  230.                   ENDDO
  231.                ELSEIF (IKAS.EQ.3) THEN
  232.                   DO ICOMP=1,NCOMP
  233.                      I=INVSYM(1,ICOMP,IDIM)
  234.                      J=INVSYM(2,ICOMP,IDIM)
  235.                      BB(ICOMP,INNOE)=R(I,J)
  236.                      IBIN(ICOMP,INNOE)=1
  237.                   ENDDO
  238.                ELSEIF (IKAS.EQ.4) THEN
  239.                   DO ICOMP=1,NCOMP
  240.                      I=INVGEN(1,ICOMP,IDIM)
  241.                      J=INVGEN(2,ICOMP,IDIM)
  242.                      BB(ICOMP,INNOE)=R(I,J)
  243.                      IBIN(ICOMP,INNOE)=1
  244.                   ENDDO
  245.                ELSE
  246.                   write(ioimp,*) 'ikas=',ikas,' ????'
  247.                   goto 9999
  248.                endif
  249.             ENDIF
  250. * Fin de la boucle sur les noeuds
  251.          ENDDO
  252. *     segprt,mtrav
  253.       ELSEIF (TYCHA.EQ.'MCHAML  ') THEN
  254.          N1=ICHAML(/1)
  255.          DO I1=1,N1
  256.             MCHAML=ICHAML(I1)
  257.             SEGACT MCHAML
  258.             N2=IELVAL(/1)
  259.             IF (N2.NE.NCOMP) THEN
  260. * Composante inexistante
  261.                CALL ERREUR(280)
  262.                RETURN
  263.             ENDIF
  264.             N1EL=0
  265.             N1PTEL=0
  266.             N2EL=0
  267.             N2PTEL=0
  268.             DO I2=1,NCOMP
  269.                IF (TYPCHE(I1).NE.'REAL*8') THEN
  270.                   MOTERR(1:8)='MCHAML'
  271.                   MOTERR(9:16)=TYPCHE(I1)
  272. * 131 2               On n'attend pas un objet de type %m1:8 de sous-type %m9:16
  273.                   CALL ERREUR(131)
  274.                   RETURN
  275.                ENDIF
  276.                MELVAL=IELVAL(I2)
  277.                SEGACT MELVAL
  278.                N1EL=MAX(N1EL,VELCHE(/2))
  279.                N1PTEL=MAX(N1PTEL,VELCHE(/1))
  280.             ENDDO
  281.             IF (NCOM1.EQ.NCOMP) THEN
  282.                SEGINI,MCHAM1=MCHAML
  283.             ELSE
  284.                N2=NCOM1
  285.                SEGINI,MCHAM1
  286.                IF(NCOM1.EQ.1) THEN
  287.                   MCHAM1.NOMCHE(1)='SCAL'
  288.                   MCHAM1.TYPCHE(1)='REAL*8'
  289.                ELSEIF (IOTENS.EQ.11) THEN
  290.                   ICMP=0
  291.                   CCCOMP='SI'
  292.                   DO i=1,IDIM
  293.                      ICMP=ICMP+1
  294.                      WRITE(CCCOMP(3:3),FMT='(I1)') I
  295.                      WRITE(CCCOMP(4:4),FMT='(I1)') I
  296.                      ICMP=ICMP+1
  297.                      MCHAM1.NOMCHE(ICMP)=CCCOMP
  298.                      MCHAM1.TYPCHE(ICMP)='REAL*8'
  299.                   ENDDO
  300.                   CCCOMP='CO'
  301.                   DO i=1,IDIM
  302.                      WRITE(CCCOMP(4:4),FMT='(I1)') I
  303.                      DO j=1,IDIM
  304.                         WRITE(CCCOMP(3:3),FMT='(I1)') J
  305.                         ICMP=ICMP+1
  306.                         MCHAM1.NOMCHE(ICMP)=CCCOMP
  307.                         MCHAM1.TYPCHE(ICMP)='REAL*8'
  308.                      ENDDO
  309.                   ENDDO
  310.                ELSE
  311.                   write(ioimp,*) 'CHAM COMP ncom1=',ncom1,' ????'
  312.                   goto 9999
  313.                ENDIF
  314.             ENDIF
  315.             NNIN=NCOMP
  316.             SEGINI MTRCML
  317.             DO I2=1,NCOMP
  318.                MELVAL=IELVAL(I2)
  319.                J1PTEL=VELCHE(/1)
  320.                J1EL=VELCHE(/2)
  321.                DO I1EL=1,N1EL
  322.                   DO I1PTEL=1,N1PTEL
  323.                      CC(I2,I1PTEL,I1EL)=VELCHE(MIN(I1PTEL,J1PTEL)
  324.      $                    ,MIN(I1EL,J1EL))
  325.                   ENDDO
  326.                ENDDO
  327.             ENDDO
  328. *
  329.             DO I1EL=1,N1EL
  330.                DO I1PTEL=1,N1PTEL
  331. * Remplissage des tableaux de travail
  332.                   IF (IKAS.EQ.1.OR.IKAS.EQ.2) THEN
  333.                      DO I=1,NCOMP
  334.                         A(I,1)=CC(I,I1PTEL,I1EL)
  335.                      ENDDO
  336.                   ELSEIF (IKAS.EQ.3) THEN
  337.                      DO J=1,IDIM
  338.                         DO I=1,IDIM
  339.                            A(I,J)=CC(IDXSYM(I,J,IDIM),I1PTEL,I1EL)
  340.                         ENDDO
  341.                      ENDDO
  342.                   ELSEIF (IKAS.EQ.4) THEN
  343.                      DO J=1,IDIM
  344.                         DO I=1,IDIM
  345.                            A(I,J)=CC(IDXGEN(I,J,IDIM),I1PTEL,I1EL)
  346.                         ENDDO
  347.                      ENDDO
  348.                   ELSE
  349.                      write(ioimp,*) 'ikas=',ikas,' ????'
  350.                      GOTO 9999
  351.                   endif
  352. *
  353. * Faire les opérations tensorielles
  354. *
  355.                   CALL TENS2(IOTENS,IKAS,A,D,R)
  356.                   IF (IERR.NE.0) RETURN
  357. *
  358.                   IF (NCOM1.NE.NCOMP) THEN
  359.                      IF(NCOM1.EQ.1) THEN
  360.                         CC(1,I1PTEL,I1EL)=R(1,1)
  361.                      ELSEIF (IOTENS.EQ.11) THEN
  362.                         ICMP=0
  363.                         DO i=1,IDIM
  364.                            ICMP=ICMP+1
  365.                            CC(ICMP,I1PTEL,I1EL)=D(i)
  366.                         ENDDO
  367.                         DO i=1,IDIM
  368.                            DO j=1,IDIM
  369.                               ICMP=ICMP+1
  370.                               CC(ICMP,I1PTEL,I1EL)=R(j,i)
  371.                            ENDDO
  372.                         ENDDO
  373.                      ELSE
  374.                         write(ioimp,*) 'CHAM VAL ncom1=',ncom1,' ????'
  375.                         goto 9999
  376.                      ENDIF
  377.  
  378.                   ELSE
  379.                      IF (IKAS.EQ.1.OR.IKAS.EQ.2) THEN
  380.                         DO I=1,NCOMP
  381.                            CC(I,I1PTEL,I1EL)=R(I,1)
  382.                         ENDDO
  383.                      ELSEIF (IKAS.EQ.3) THEN
  384.                         DO ICOMP=1,NCOMP
  385.                            I=INVSYM(1,ICOMP,IDIM)
  386.                            J=INVSYM(2,ICOMP,IDIM)
  387.                            CC(ICOMP,I1PTEL,I1EL)=R(I,J)
  388.                         ENDDO
  389.                      ELSEIF (IKAS.EQ.4) THEN
  390.                         DO ICOMP=1,NCOMP
  391.                            I=INVGEN(1,ICOMP,IDIM)
  392.                            J=INVGEN(2,ICOMP,IDIM)
  393.                            CC(ICOMP,I1PTEL,I1EL)=R(I,J)
  394.                         ENDDO
  395.                      ELSE
  396.                         write(ioimp,*) 'ikas=',ikas,' ????'
  397.                         goto 9999
  398.                      endif
  399.                   ENDIF
  400. * Fin de la boucle sur les noeuds et éléments
  401.                ENDDO
  402. * Fin de la boucle sur les éléments
  403.             ENDDO
  404.             DO I2=1,NCOM1
  405.                SEGINI MELVA1
  406.                DO I1EL=1,N1EL
  407.                   DO I1PTEL=1,N1PTEL
  408.                      MELVA1.VELCHE(I1PTEL,I1EL)=CC(I2,I1PTEL,I1EL)
  409.                   ENDDO
  410.                ENDDO
  411.                MCHAM1.IELVAL(I2)=MELVA1
  412.             ENDDO
  413.             SEGSUP MTRCML
  414.             MCHEL1.ICHAML(I1)=MCHAM1
  415. *     Fin de la boucle sur les sous-zones
  416.          ENDDO
  417.       ELSE
  418.          GOTO 9998
  419.       ENDIF
  420. *
  421. * Transformation des objets de travail en champs
  422. *
  423.       IF (TYCHA.EQ.'CHPOINT ') THEN
  424.          IF (NCOM1.NE.NCOMP) THEN
  425. * Le résultat est dans MTRAV1
  426.             CALL CRECHP(MTRAV1,MCHPO1)
  427.             SEGSUP MTRAV1
  428.          ELSE
  429.             IF (IOTENS.GE.1.AND.IOTENS.LE.4) INCO(1)='SCAL'
  430.             CALL CRECHP(MTRAV,MCHPO1)
  431.          ENDIF
  432.          SEGSUP MTRAV
  433.          ICHA1=MCHPO1
  434.       ELSEIF (TYCHA.EQ.'MCHAML  ') THEN
  435.          ICHA1=MCHEL1
  436.       ELSE
  437.          GOTO 9998
  438.       ENDIF
  439. *
  440. * Normal termination
  441. *
  442.       RETURN
  443. *
  444. * Format handling
  445. *
  446. *
  447. * Error handling
  448. *
  449.  9997 CONTINUE
  450. * 980 2
  451. * L'objet %m1:8 n'a pas le bon nombre de composantes
  452.       MOTERR(1:8)='LISTMOTS'
  453.       CALL ERREUR(980)
  454.       RETURN
  455.  9998 CONTINUE
  456. * On ne veut pas d'objet de type %m1:8
  457.       MOTERR(1:8)=TYCHA
  458.       CALL ERREUR(39)
  459.       RETURN
  460.  9999 CONTINUE
  461.       MOTERR(1:8)='TENS1'
  462.       call erreur(1039)
  463.       return
  464. *
  465.  
  466. * End of subroutine TENS1
  467. *
  468.       END
  469.  
  470.  

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