Sources Esope | Cast3M

cqtgr1
C CQTGR1    SOURCE    CHAT      05/01/12    22:27:34     5004
      SUBROUTINE CQTGR1(IGAU,MELE,NBNN,LRE,IFOU,NGRA,XE,SHPTOT,
     1                   XDDL,SHP,BGR,DJAC,GRADI)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
*****************************************************************
*      SOUS-PROGRAMME DE L'OPERATEUR GRADIENT
*
*      CALCUL DES GRADIENTS DE TEMPERATURE POUR LES ELEMENTS COQ2,COQ3
*
*
* ENTREES :
*
*    IGAU   NUMERO DU POINT DE GAUSS
*    MELE   NUMERO DE L'ELEMENT DANS NOMTP
*    NBNN   NOMBRE DE NOEUDS
*    LRE    NOMBRE DE COLONNES DE LA MATRICE BGR
*    IFOU   IFOUR DE CCOPTIO
*    NGRA   NOMBRE DE COMPOSANTES DE GRADIENTS
*    XE     COORDONNEES DE L'ELEMENT
*  SHPTOT   FONCTIONS DE FORMES ET DERIVEES
*   SHP     TABLEAU DE TRAVAIL
*   BGR     TABLEAU DE TRAVAIL
*
* SORTIES  :
*
*   GRADI   TABLEAU CONTENANT LES COMPOSANTES DU GRADIENT
*   DJAC    JACOBIEN
*
*  AUTEUR   : P.DOWLATYARI   27/05/91
************************************************************************
 
-INC PPARAM
-INC CCOPTIO
       DIMENSION XE(3,*),BGR(NGRA,*),SHP(6,*),SHPTOT(6,NBNN,*)
       DIMENSION XDDL(*),GRADI(*)
       DIMENSION COSD1(3),COSD2(3),COSD3(3),XE1(3,3)
*
       CALL ZERO(BGR,NGRA,LRE)
*
       IF(MELE.EQ.44)THEN
*
*      COQ2
*
         IF(IFOU.EQ.0)THEN
*
*     LA LONGUEUR DE L'ELEMENT
*
           D=SQRT((XE(1,2)-XE(1,1))**2.D0 +(XE(2,2)-XE(2,1))**2.D0)
           DJAC=D
           IF (D.EQ.0.D0) RETURN
           BGR(1,1)=-1.D0/D
           BGR(1,2)= 1.D0/D
          ELSE
*
*         ELEMENT AXISYM. DE FOURIER OU AUTRE
*         OPTION NON DISPONIBLE ACTUELLEMENT
*
            CALL ERREUR(19)
             RETURN
           ENDIF
         ELSE
*
*            COQ3
*
*
*       CALCUL DES COORDONNEES DES NOEUDS DANS LE REPERE LOCAL DE L'
*       ELEMENT COQUE
*
           CALL ZERO(XE1,3,3)
           DO 10 I=1,3
            COSD1(I)=XE(I,2)-XE(I,1)
            COSD2(I)=XE(I,3)-XE(I,1)
  10       CONTINUE
*
           COSD3(1)=COSD1(2)*COSD2(3)-COSD1(3)*COSD2(2)
           COSD3(2)=COSD1(3)*COSD2(1)-COSD1(1)*COSD2(3)
           COSD3(3)=COSD1(1)*COSD2(2)-COSD1(2)*COSD2(1)
*
           COSD1L=SQRT(COSD1(1)*COSD1(1)+COSD1(2)*COSD1(2)+
     .               COSD1(3)*COSD1(3))
           COSD3L=SQRT(COSD3(1)*COSD3(1)+COSD3(2)*COSD3(2)+
     .              COSD3(3)*COSD3(3))
*
           DO 20 I=1,3
             COSD1(I)=COSD1(I)/COSD1L
             COSD3(I)=COSD3(I)/COSD3L
   20      CONTINUE
*
           COSD2(1)=COSD3(2)*COSD1(3)-COSD3(3)*COSD1(2)
           COSD2(2)=COSD3(3)*COSD1(1)-COSD3(1)*COSD1(3)
           COSD2(3)=COSD3(1)*COSD1(2)-COSD3(2)*COSD1(1)
*
           DO 30 NOE=1,NBNN
           DO 30  I=1,3
            XE1(1,NOE)=XE1(1,NOE)+XE(I,NOE)*COSD1(I)
            XE1(2,NOE)=XE1(2,NOE)+XE(I,NOE)*COSD2(I)
   30      CONTINUE
*
           DO 40 NP=1,NBNN
             SHP(1,NP)=SHPTOT(1,NP,IGAU)
             SHP(2,NP)=SHPTOT(2,NP,IGAU)
             SHP(3,NP)=SHPTOT(3,NP,IGAU)
  40        CONTINUE
*
           CALL JACOBI(XE1,SHP,2,NBNN,DJAC)
           K=1
           DO 50 NP=1,NBNN
             BGR(1,K)=SHP(2,NP)
             BGR(2,K)=SHP(3,NP)
             K=K+1
  50       CONTINUE
*
          ENDIF
*
*      CALCUL DES COMPOSANTES DU GRADIENT
*
          LRE1=LRE/3
          DO 60  IA=1,NGRA
             CC=0.D0
              DO 70 IB=1,LRE1
                IBB=(IB-1)*3+2
                CC=CC+XDDL(IBB)*BGR(IA,IB)
  70          CONTINUE
             GRADI(IA)=CC
  60      CONTINUE
*
          RETURN
          END