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   25/07/28    21:15:11     12338          
  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.       CHARACTER*8 TYCHA,CCCOMP
  43.       segment idxcom(ncomp)
  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.          NOID=1
  146. * Extrayons les composantes intéressantes
  147. *         write(ioimp,*) 'Avant excoc2' ;
  148. *         CALL ECROBJ('MCHAML  ',MCHELM)
  149. *         call prlist
  150. *         segprt,MLMOTS
  151.          CALL EXCOC2(MCHELM,MLMOTS,MCHEL1,MLMOTS,NOID)
  152.          IF (IERR.NE.0) RETURN
  153. *         write(ioimp,*) 'Apres excoc2' ;
  154. *         CALL ECROBJ('MCHAML  ',MCHEL1)
  155. *         call prlist
  156.          CALL ACTOBJ('MCHAML  ',MCHEL1,1)
  157.          SEGACT,MCHEL1*MOD
  158.          SEGACT MLMOTS
  159.       ELSE
  160.          GOTO 9998
  161.       ENDIF
  162. *      segprt,mtrav
  163. *
  164. *     write(ioimp,*) 'IOTENS=',IOTENS
  165. *
  166. * Un peu de gestion d'erreur avant la boucle...
  167. *
  168. * Seuls NORM2 et NORMINF disponibles pour les vecteurs
  169.       IF (IKAS.EQ.2.AND.IOTENS.GT.2) THEN
  170.          CALL ERREUR(803)
  171.          RETURN
  172.       ENDIF
  173.       MOTERR(1:8)=TYCHA
  174. *
  175. * La boucle sur les noeuds ou les éléments : on pourra paralléliser ici
  176. *
  177.       IF (TYCHA.EQ.'CHPOINT ') THEN
  178.          DO INNOE=1,NNNOE
  179. * Remplissage des tableaux de travail
  180.             IF (IKAS.EQ.1.OR.IKAS.EQ.2) THEN
  181.                DO I=1,NCOMP
  182.                   A(I,1)=BB(I,INNOE)
  183.                ENDDO
  184.             ELSEIF (IKAS.EQ.3) THEN
  185.                DO J=1,IDIM
  186.                   DO I=1,IDIM
  187.                      A(I,J)=BB(IDXSYM(I,J,IDIM),INNOE)
  188.                   ENDDO
  189.                ENDDO
  190.             ELSEIF (IKAS.EQ.4) THEN
  191.                DO J=1,IDIM
  192.                   DO I=1,IDIM
  193.                      A(I,J)=BB(IDXGEN(I,J,IDIM),INNOE)
  194.                   ENDDO
  195.                ENDDO
  196.             ELSE
  197.                write(ioimp,*) 'ikas=',ikas,' ????'
  198.                MOTERR(1:8)='TENS1'
  199.                call erreur(1039)
  200.                return
  201.             endif
  202. *     On prépare le message d'erreur 49 en cas de déterminant nul pour
  203. *     l'inversion d'une matrice...
  204.             IF (IOTENS.EQ.5) INTERR(1)=IGEO(INNOE)
  205. *
  206. * Faire les opérations tensorielles
  207. *
  208.             CALL TENS2(IOTENS,IKAS,A,D,R)
  209.             IF (IERR.NE.0) RETURN
  210. *
  211.             IF (NCOM1.NE.NCOMP) THEN
  212.                IF(NCOM1.EQ.1) THEN
  213.                   MTRAV1.BB(1,INNOE)=R(1,1)
  214.                   MTRAV1.IBIN(1,INNOE)=1
  215.                ELSEIF (IOTENS.EQ.11) THEN
  216.                   ICMP=0
  217.                   DO i=1,IDIM
  218.                      ICMP=ICMP+1
  219.                      MTRAV1.BB(ICMP,INNOE)=D(i)
  220.                      MTRAV1.IBIN(ICMP,INNOE)=1
  221.                   ENDDO
  222.                   DO i=1,IDIM
  223.                      DO j=1,IDIM
  224.                         ICMP=ICMP+1
  225.                         MTRAV1.BB(ICMP,INNOE)=R(j,i)
  226.                         MTRAV1.IBIN(ICMP,INNOE)=1
  227.                      ENDDO
  228.                   ENDDO
  229.                ELSE
  230.                   write(ioimp,*) 'CHPO VAL ncom1=',ncom1,' ????'
  231.                   goto 9999
  232.                ENDIF
  233.             ELSE
  234.                IF (IKAS.EQ.1.OR.IKAS.EQ.2) THEN
  235.                   DO I=1,NCOMP
  236.                      BB(I,INNOE)=R(I,1)
  237.                      IBIN(I,INNOE)=1
  238.                   ENDDO
  239.                ELSEIF (IKAS.EQ.3) THEN
  240.                   DO ICOMP=1,NCOMP
  241.                      I=INVSYM(1,ICOMP,IDIM)
  242.                      J=INVSYM(2,ICOMP,IDIM)
  243.                      BB(ICOMP,INNOE)=R(I,J)
  244.                      IBIN(ICOMP,INNOE)=1
  245.                   ENDDO
  246.                ELSEIF (IKAS.EQ.4) THEN
  247.                   DO ICOMP=1,NCOMP
  248.                      I=INVGEN(1,ICOMP,IDIM)
  249.                      J=INVGEN(2,ICOMP,IDIM)
  250.                      BB(ICOMP,INNOE)=R(I,J)
  251.                      IBIN(ICOMP,INNOE)=1
  252.                   ENDDO
  253.                ELSE
  254.                   write(ioimp,*) 'ikas=',ikas,' ????'
  255.                   goto 9999
  256.                endif
  257.             ENDIF
  258. * Fin de la boucle sur les noeuds
  259.          ENDDO
  260. *     segprt,mtrav
  261.       ELSEIF (TYCHA.EQ.'MCHAML  ') THEN
  262.          N1=MCHEL1.ICHAML(/1)
  263.          DO I1=1,N1
  264.             MCHAML=MCHEL1.ICHAML(I1)
  265.             SEGACT MCHAML
  266.             N2=IELVAL(/1)
  267.             IF (N2.NE.NCOMP) THEN
  268.                write(ioimp,*) ' TENS1 1 MCHEL1,MCHAML,N2,NCOMP=',MCHEL1
  269.      $              ,MCHAML,N2,NCOMP
  270.                write(ioimp,'(10(1X,A))') 'NOMCHE=',(NOMCHE(ii),ii=1,N2)
  271.                write(ioimp,'(10(1X,A))') 'MLMOTS=',(MOTS(ii),ii=1,NCOMP)
  272. * Composante inexistante
  273.                CALL ERREUR(280)
  274.                RETURN
  275.             ENDIF
  276. * Les composantes du chaml ne sont pas forcément dans l'ordre de MLMOTS...
  277.             segini idxcom
  278.             DO ICOMP=1,NCOMP
  279.                call place(nomche(1),ncomp,i2,mots(icomp))
  280.                if (i2.le.0) then
  281.                   write(ioimp,*) ' TENS1 2 MCHEL1,MCHAML,NCOMP,I2='
  282.      $                 ,MCHEL1,MCHAML,NCOMP,I2
  283.                   write(ioimp,'(10(1X,A))') 'NOMCHE=',(NOMCHE(ii),ii=1
  284.      $                 ,NCOMP)
  285.                   write(ioimp,'(10(1X,A))') 'MLMOTS=',(MOTS(ii),ii=1
  286.      $                 ,NCOMP)
  287.                   goto 9999
  288.                endif
  289.                idxcom(ICOMP)=i2
  290.             ENDDO
  291.             N1EL=0
  292.             N1PTEL=0
  293.             N2EL=0
  294.             N2PTEL=0
  295.             DO ICOMP=1,NCOMP
  296.                i2=idxcom(icomp)
  297.                IF (TYPCHE(I2).NE.'REAL*8') THEN
  298.                   MOTERR(1:8)='MCHAML'
  299.                   MOTERR(9:16)=TYPCHE(I2)
  300. * 131 2               On n'attend pas un objet de type %m1:8 de sous-type %m9:16
  301.                   CALL ERREUR(131)
  302.                   RETURN
  303.                ENDIF
  304.                MELVAL=IELVAL(I2)
  305.                SEGACT MELVAL
  306.                N1EL=MAX(N1EL,VELCHE(/2))
  307.                N1PTEL=MAX(N1PTEL,VELCHE(/1))
  308.             ENDDO
  309.             IF (NCOM1.EQ.NCOMP) THEN
  310.                SEGINI,MCHAM1=MCHAML
  311.             ELSE
  312.                N2=NCOM1
  313.                SEGINI,MCHAM1
  314.                IF(NCOM1.EQ.1) THEN
  315.                   MCHAM1.NOMCHE(1)='SCAL'
  316.                   MCHAM1.TYPCHE(1)='REAL*8'
  317.                ELSEIF (IOTENS.EQ.11) THEN
  318.                   ICMP=0
  319.                   CCCOMP='SI'
  320.                   DO i=1,IDIM
  321.                      ICMP=ICMP+1
  322.                      WRITE(CCCOMP(3:3),FMT='(I1)') I
  323.                      WRITE(CCCOMP(4:4),FMT='(I1)') I
  324.                      MCHAM1.NOMCHE(ICMP)=CCCOMP
  325.                      MCHAM1.TYPCHE(ICMP)='REAL*8'
  326.                   ENDDO
  327.                   CCCOMP='CO'
  328.                   DO i=1,IDIM
  329.                      WRITE(CCCOMP(4:4),FMT='(I1)') I
  330.                      DO j=1,IDIM
  331.                         WRITE(CCCOMP(3:3),FMT='(I1)') J
  332.                         ICMP=ICMP+1
  333.                         MCHAM1.NOMCHE(ICMP)=CCCOMP
  334.                         MCHAM1.TYPCHE(ICMP)='REAL*8'
  335.                      ENDDO
  336.                   ENDDO
  337.                ELSE
  338.                   write(ioimp,*) 'CHAM COMP ncom1=',ncom1,' ????'
  339.                   goto 9999
  340.                ENDIF
  341.             ENDIF
  342.             NNIN=max(NCOM1,NCOMP)
  343.             SEGINI MTRCML
  344.             DO ICOMP=1,NCOMP
  345.                i2=idxcom(icomp)
  346.                MELVAL=IELVAL(I2)
  347.                J1PTEL=VELCHE(/1)
  348.                J1EL=VELCHE(/2)
  349.                DO I1EL=1,N1EL
  350.                   DO I1PTEL=1,N1PTEL
  351.                      CC(ICOMP,I1PTEL,I1EL)=VELCHE(MIN(I1PTEL,J1PTEL)
  352.      $                    ,MIN(I1EL,J1EL))
  353.                   ENDDO
  354.                ENDDO
  355.             ENDDO
  356. *
  357.             DO I1EL=1,N1EL
  358.                DO I1PTEL=1,N1PTEL
  359. * Remplissage des tableaux de travail
  360.                   IF (IKAS.EQ.1.OR.IKAS.EQ.2) THEN
  361.                      DO I=1,NCOMP
  362.                         A(I,1)=CC(I,I1PTEL,I1EL)
  363.                      ENDDO
  364.                   ELSEIF (IKAS.EQ.3) THEN
  365.                      DO J=1,IDIM
  366.                         DO I=1,IDIM
  367.                            A(I,J)=CC(IDXSYM(I,J,IDIM),I1PTEL,I1EL)
  368.                         ENDDO
  369.                      ENDDO
  370.                   ELSEIF (IKAS.EQ.4) THEN
  371.                      DO J=1,IDIM
  372.                         DO I=1,IDIM
  373.                            A(I,J)=CC(IDXGEN(I,J,IDIM),I1PTEL,I1EL)
  374.                         ENDDO
  375.                      ENDDO
  376.                   ELSE
  377.                      write(ioimp,*) 'ikas=',ikas,' ????'
  378.                      GOTO 9999
  379.                   endif
  380. *
  381. * Faire les opérations tensorielles
  382. *
  383.                   CALL TENS2(IOTENS,IKAS,A,D,R)
  384.                   IF (IERR.NE.0) RETURN
  385. *
  386.                   IF (NCOM1.NE.NCOMP) THEN
  387.                      IF(NCOM1.EQ.1) THEN
  388.                         CC(1,I1PTEL,I1EL)=R(1,1)
  389.                      ELSEIF (IOTENS.EQ.11) THEN
  390.                         ICMP=0
  391.                         DO i=1,IDIM
  392.                            ICMP=ICMP+1
  393.                            CC(ICMP,I1PTEL,I1EL)=D(i)
  394.                         ENDDO
  395.                         DO i=1,IDIM
  396.                            DO j=1,IDIM
  397.                               ICMP=ICMP+1
  398.                               CC(ICMP,I1PTEL,I1EL)=R(j,i)
  399.                            ENDDO
  400.                         ENDDO
  401.                      ELSE
  402.                         write(ioimp,*) 'CHAM VAL ncom1=',ncom1,' ????'
  403.                         goto 9999
  404.                      ENDIF
  405.  
  406.                   ELSE
  407.                      IF (IKAS.EQ.1.OR.IKAS.EQ.2) THEN
  408.                         DO I=1,NCOMP
  409.                            CC(I,I1PTEL,I1EL)=R(I,1)
  410.                         ENDDO
  411.                      ELSEIF (IKAS.EQ.3) THEN
  412.                         DO ICOMP=1,NCOMP
  413.                            I=INVSYM(1,ICOMP,IDIM)
  414.                            J=INVSYM(2,ICOMP,IDIM)
  415.                            CC(ICOMP,I1PTEL,I1EL)=R(I,J)
  416.                         ENDDO
  417.                      ELSEIF (IKAS.EQ.4) THEN
  418.                         DO ICOMP=1,NCOMP
  419.                            I=INVGEN(1,ICOMP,IDIM)
  420.                            J=INVGEN(2,ICOMP,IDIM)
  421.                            CC(ICOMP,I1PTEL,I1EL)=R(I,J)
  422.                         ENDDO
  423.                      ELSE
  424.                         write(ioimp,*) 'ikas=',ikas,' ????'
  425.                         goto 9999
  426.                      endif
  427.                   ENDIF
  428. * Fin de la boucle sur les noeuds et éléments
  429.                ENDDO
  430. * Fin de la boucle sur les éléments
  431.             ENDDO
  432.             DO ICOM1=1,NCOM1
  433.                SEGINI MELVA1
  434.                DO I1EL=1,N1EL
  435.                   DO I1PTEL=1,N1PTEL
  436.                      MELVA1.VELCHE(I1PTEL,I1EL)=CC(ICOM1,I1PTEL,I1EL)
  437.                   ENDDO
  438.                ENDDO
  439.                IF (NCOMP.EQ.NCOM1) THEN
  440.                   I2=IDXCOM(ICOM1)
  441.                   MCHAM1.IELVAL(I2)=MELVA1
  442.                ELSE
  443.                   MCHAM1.IELVAL(ICOM1)=MELVA1
  444.                ENDIF
  445.             ENDDO
  446.             segsup idxcom
  447.             SEGSUP MTRCML
  448.             MCHEL1.ICHAML(I1)=MCHAM1
  449. *     Fin de la boucle sur les sous-zones
  450.          ENDDO
  451.       ELSE
  452.          GOTO 9998
  453.       ENDIF
  454. *
  455. * Transformation des objets de travail en champs
  456. *
  457.       IF (TYCHA.EQ.'CHPOINT ') THEN
  458.          IF (NCOM1.NE.NCOMP) THEN
  459. * Le résultat est dans MTRAV1
  460.             CALL CRECHP(MTRAV1,MCHPO1)
  461.             SEGSUP MTRAV1
  462.          ELSE
  463.             IF (IOTENS.GE.1.AND.IOTENS.LE.4) INCO(1)='SCAL'
  464.             CALL CRECHP(MTRAV,MCHPO1)
  465.          ENDIF
  466.          SEGSUP MTRAV
  467.          ICHA1=MCHPO1
  468.       ELSEIF (TYCHA.EQ.'MCHAML  ') THEN
  469.          ICHA1=MCHEL1
  470.       ELSE
  471.          GOTO 9998
  472.       ENDIF
  473. *
  474. * Normal termination
  475. *
  476.       RETURN
  477. *
  478. * Format handling
  479. *
  480. *
  481. * Error handling
  482. *
  483.  9997 CONTINUE
  484. * 980 2
  485. * L'objet %m1:8 n'a pas le bon nombre de composantes
  486.       MOTERR(1:8)='LISTMOTS'
  487.       CALL ERREUR(980)
  488.       RETURN
  489.  9998 CONTINUE
  490. * On ne veut pas d'objet de type %m1:8
  491.       MOTERR(1:8)=TYCHA
  492.       CALL ERREUR(39)
  493.       RETURN
  494.  9999 CONTINUE
  495.       MOTERR(1:8)='TENS1'
  496.       call erreur(1039)
  497.       return
  498. *
  499.  
  500. * End of subroutine TENS1
  501. *
  502.       END
  503.  
  504.  

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