Sources Esope | Cast3M

tailca
C TAILCA    SOURCE    FANDEUR   10/08/30    21:15:47     6729
C                     responsable : Mr MILLARD
 
C=======================================================================
C
       SUBROUTINE TAILCA(IPTRA,POI,IDIM,IFOUR,NBNO,MELE,IELE,MFR, IRET)
C
C=======================================================================
C
C     CALCULE LE GRADIENT DE LA GEOMETRIE DEFORMEE
C     CALCULE LA TRANSFORMEE DES VECTEURS UNITAIRES DE LA
C     GEOMETRIE DE REFERENCE DANS LE REPERE GLOBAL
C   ENTREES
C     IPTRA=POINTEUR SUR UN SEGMENT CONTENANT LES INFOS :
C   -XEL(3,NBNO)=COORDONNEES LOCALES DE L ELEMENT
C   -SHP(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE
C   -SHPQSI(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE de (1 0 0)
C   _SHPETA(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE de (0 1 0)
C   -SHPDZE(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE de (0 0 1)
C     POI = POI POUR L INTEGRATION DU GRADIANT AU PT DE GAUSS
C     IDIM=DIMENSION
C     NBNO=NOMBRE DE NOEUDS
C
C   SORTIES
C    VCOMP = COMPOSANTES DU PARAMETRE DE TAILLE DE L ELEMENT
C            ET DU POINT DE GAUSS COURANT
C=================================================================
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (XZER=0.D0,UNDEMI=.5D0,UN=1.D0,DEUX=2.D0,
     &           x1sur3=0.33333333333333333333333333333333333333333,
     &           x1sra2=0.70710678118654752440084436210484904)
C Note : x1sur3 = 1./3. et x1sra2 = 1./SQRT(2.)
C
      SEGMENT MTRA1
        REAL*8 XEL(3,NBNN)
        REAL*8 VCOMP(NCOMP)
        REAL*8 SHP(6,NBNN),SHPZER(6,NBNN)
        REAL*8 SHPQSI(6,NBNN),SHPETA(6,NBNN),SHPDZE(6,NBNN)
C*      REAL*8 SHPGAU(6,NBNN)
      ENDSEGMENT
 
      NBNN=NBNO
      MTRA1=IPTRA
 
C=============================
C     CAS 2 DIMENSIONS
C  OU CAS 3D COQUES MINCES
C=============================
      IF (IDIM.EQ.3.AND.(MFR.EQ.1.OR.MFR.EQ.31)) GOTO 30
C=============================
 
      POIG=ABS(POI)
 
      DXDQSI = XZER
      DXDETA = XZER
      DYDQSI = XZER
      DYDETA = XZER
      XORI = XZER
      YORI = XZER
      XQSI = XZER
      YQSI = XZER
      XETA = XZER
      YETA = XZER
 
      DO 100 I=1,NBNN
        DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
        DXDETA=DXDETA+SHP(3,I)*XEL(1,I)
        DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I)
        DYDETA=DYDETA+SHP(3,I)*XEL(2,I)
        XORI=XORI+SHPZER(1,I)*XEL(1,I)
        YORI=YORI+SHPZER(1,I)*XEL(2,I)
        XQSI=XQSI+SHPQSI(1,I)*XEL(1,I)
        YQSI=YQSI+SHPQSI(1,I)*XEL(2,I)
        XETA=XETA+SHPETA(1,I)*XEL(1,I)
        YETA=YETA+SHPETA(1,I)*XEL(2,I)
 100  CONTINUE
 
      UX=XQSI-XORI
      UY=YQSI-YORI
      VX=XETA-XORI
      VY=YETA-YORI
      WX=XQSI-XETA
      WY=YQSI-YETA
 
      RNU=SQRT(Ux*Ux+Uy*Uy)
      RNV=SQRT(Vx*Vx+Vy*Vy)
      RNW=SQRT(Wx*Wx+Wy*Wy)
C En cas de vecteur(s) nul(s), on met sa(leur) norme(s) a un.
C Le determinant "det" calcule plus loin sera nul d'où une erreur !
      IF (RNU.EQ.XZER) RNU=UN
      IF (RNV.EQ.XZER) RNV=UN
      IF (RNW.EQ.XZER) RNW=UN
C
C------Calcule de la norme du gradient dans les directions associees-----
C------------------------------a QSI ETA DZE-----------------------------
C------------------et integration par la methode de Gauss----------------
C
      XLQSI=SQRT(POIG*(DXDQSI*DXDQSI+DYDQSI*DYDQSI))
      YLETA=SQRT(POIG*(DXDETA*DXDETA+DYDETA*DYDETA))
      ZLDZE=UN
C
C========================================================================
C--------------CALCUL DU TENSEUR DE PARAMETRE TAILLE---------------------
C========================================================================
C
C----------CALCUL DE LA MATRICE DE PASSAGE ET DE SON INVERSE-------------
C
C==========================================================
C---------------lorsqu il s agit de triangle,--------------
C---------- les deux + petit cotes sont pris en compte-----
C==========================================================
      IF (IELE.EQ.4.OR.IELE.EQ.6) THEN
         IF (IELE.EQ.4) THEN
           rconst = x1sra2
         ELSE
C        ELSE IF (IELE.EQ.6) THEN
           rconst = UNDEMI
         ENDIF
         IF (RNW.LT.RNU.OR.RNW.LT.RNV) THEN
            IF (RNU.LT.RNV) THEN
               YLETA=SQRT(POIG)*RNW
               UUX=UX
               UUY=UY
               VVX=WX
               VVY=WY
            ELSE
               XLQSI=SQRT(POIG)*RNW
               UUX=WX
               UUY=WY
               VVX=VX
               VVY=VY
            ENDIF
            RNUU=SQRT(UUX*UUX+UUY*UUY)
            RNVV=SQRT(VVX*VVX+VVY*VVY)
            a=UUX/RNUU
            b=UUY/RNUU
            c=VVX/RNVV
            d=VVY/RNVV
         ELSE
            a=UX/RNU
            b=UY/RNU
            c=VX/RNV
            d=VY/RNV
         ENDIF
         XLQSI = XLQSI*rconst
         YLETA = YLETA*rconst
C------------------------------------------
C pas de modification pour les quadrangles
C------------------------------------------
      ELSE
          a=UX/RNU
          b=UY/RNU
          c=VX/RNV
          d=VY/RNV
      ENDIF
C
C-------------------------------------------
C      calcul du produit scalaire
C-------------------------------------------
C
      XY=UNDEMI*(a*c+b*d)*SQRT(XLQSI*YLETA)
C
C-----------------------------------------------------------------
C           CALCUL DES COEFFICIENTS DANS LE REPERE GLOBAL
C-----------------------------------------------------------------
C      write (*,*) 'XLQSI : ',XLQSI,' YLETA : ',YLETA
C      write (*,*)
C      write (*,*) 'RNU : ',RNU,' RNV : ',RNV,' RNW : ',RNW
C      write (*,*) 'a : ',a,'   c : ',c
C      write (*,*) 'b : ',b,'   d : ',d
C      write (*,*) '================================================='
C
      det=a*d-b*c
      IF (det.EQ.XZER) GOTO 999
C
      XLXX=(d*d*XLQSI+b*b*YLETA-DEUX*b*d*XY)/(det*det)
      XLXY=(-c*d*XLQSI-b*a*YLETA+(b*c+a*d)*XY)/(det*det)
C--------------------------------------------
C               XLYX=XLXY
C--------------------------------------------
      YLYY=(c*c*XLQSI+a*a*YLETA-DEUX*a*c*XY)/(det*det)
C------------------------------------------
C   dans la direction normal au plan
C   il n y a pas de probleme d objectivite
C   a l utilisateur de rentrer epsi rupture
C------------------------------------------
      ZLZZ=ZLDZE
      PXX=(d*d+b*b)/(det*det)
      PXY=-(c*d+b*a)/(det*det)
C--------------------------------------------
C               PYX=PXY
C--------------------------------------------
      PYY=(a*a+c*c)/(det*det)
 
      VCOMP(1)=XLXX
      VCOMP(2)=YLYY
      VCOMP(3)=ZLZZ
      VCOMP(4)=XLXY
      VCOMP(5)=PXX
      VCOMP(6)=PYY
      VCOMP(7)=PXY
 
      GOTO 666
 
C================================
C     CAS TRIDIMENSIONEL
C================================
 30   CONTINUE
      IF (IDIM.NE.3) GOTO 666
      IF (MFR.NE.1.AND.MFR.NE.31) GOTO 666
 
      POIG=ABS(POI)**x1sur3
 
      DXDQSI=XZER
      DXDETA=XZER
      DXDDZE=XZER
      DYDQSI=XZER
      DYDETA=XZER
      DYDDZE=XZER
      DZDQSI=XZER
      DZDETA=XZER
      DZDDZE=XZER
 
      XORI=XZER
      YORI=XZER
      ZORI=XZER
      XQSI=XZER
      YQSI=XZER
      ZQSI=XZER
      XETA=XZER
      YETA=XZER
      ZETA=XZER
      XDZE=XZER
      YDZE=XZER
      ZDZE=XZER
 
      DO 300 I=1,NBNO
C
         DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
         DXDETA=DXDETA+SHP(3,I)*XEL(1,I)
         DXDDZE=DXDDZE+SHP(4,I)*XEL(1,I)
         DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I)
         DYDETA=DYDETA+SHP(3,I)*XEL(2,I)
         DYDDZE=DYDDZE+SHP(4,I)*XEL(2,I)
         DZDQSI=DZDQSI+SHP(2,I)*XEL(3,I)
         DZDETA=DZDETA+SHP(3,I)*XEL(3,I)
         DZDDZE=DZDDZE+SHP(4,I)*XEL(3,I)
 
         XORI=XORI+SHPZER(1,I)*XEL(1,I)
         YORI=YORI+SHPZER(1,I)*XEL(2,I)
         ZORI=ZORI+SHPZER(1,I)*XEL(3,I)
         XQSI=XQSI+SHPQSI(1,I)*XEL(1,I)
         YQSI=YQSI+SHPQSI(1,I)*XEL(2,I)
         ZQSI=ZQSI+SHPQSI(1,I)*XEL(3,I)
         XETA=XETA+SHPETA(1,I)*XEL(1,I)
         YETA=YETA+SHPETA(1,I)*XEL(2,I)
         ZETA=ZETA+SHPETA(1,I)*XEL(3,I)
         XDZE=XDZE+SHPDZE(1,I)*XEL(1,I)
         YDZE=YDZE+SHPDZE(1,I)*XEL(2,I)
         ZDZE=ZDZE+SHPDZE(1,I)*XEL(3,I)
 
 300  CONTINUE
C
C------Calcule de la norme du gradient dans les directions associees-----
C------------------------------a QSI ETA DZE-----------------------------
C------------------et integration par la methode de Gauss----------------
C
      XLQSI = POIG * SQRT(DXDQSI*DXDQSI+DYDQSI*DYDQSI+DZDQSI*DZDQSI)
      YLETA = POIG * SQRT(DXDETA*DXDETA+DYDETA*DYDETA+DZDETA*DZDETA)
      ZLDZE = POIG * SQRT(DXDDZE*DXDDZE+DYDDZE*DYDDZE+DZDDZE*DZDDZE)
C
C========================================================================
C--------------CALCUL DU TENSEUR DE PARAMETRE TAILLE---------------------
C========================================================================
      UX=XQSI-XORI
      UY=YQSI-YORI
      UZ=ZQSI-ZORI
      VX=XETA-XORI
      VY=YETA-YORI
      VZ=ZETA-ZORI
      WX=XDZE-XORI
      WY=YDZE-YORI
      WZ=ZDZE-ZORI
C
      VQSI=SQRT(UX*UX+UY*UY+UZ*UZ)
      VETA=SQRT(VX*VX+VY*VY+VZ*VZ)
      VDZE=SQRT(WX*WX+WY*WY+WZ*WZ)
C En cas de vecteur(s) nul(s), on met sa(leur) norme(s) a un.
C Le determinant "det" calcule plus loin sera nul d'où une erreur !
      IF (VQSI.EQ.XZER) VQSI=UN
      IF (VETA.EQ.XZER) VETA=UN
      IF (VDZE.EQ.XZER) VDZE=UN
 
C----------CALCUL DE LA MATRICE DE PASSAGE ET DE SON INVERSE-------------
      a=UX/VQSI
      b=UY/VQSI
      c=UZ/VQSI
      d=VX/VETA
      e=VY/VETA
      f=VZ/VETA
      o=WX/VDZE
      p=WY/VDZE
      q=WZ/VDZE
C
C-----------------------------------------
C   calcul des demi produits scalaires
C-----------------------------------------
C
      XY=UNDEMI*(a*d+b*e+c*f)*SQRT(XLQSI*YLETA)
      XZ=UNDEMI*(a*o+b*p+c*q)*SQRT(XLQSI*ZLDZE)
      YZ=UNDEMI*(d*o+e*p+f*q)*SQRT(YLETA*ZLDZE)
C
C---------------------------------------------------------------
C       CALCUL DES COEFFICIENTS DANS LE REPERE GLOBAL
C---------------------------------------------------------------
      det=a*e*q+b*f*o+c*d*p-o*e*c-d*b*q-a*p*f
C
      if (det.EQ.XZER) GOTO 999
C
C----calcul de l inverse de la matrice de passage----
C
      r=(e*q-f*p)/det
      s=(c*p-b*q)/det
      t=(b*f-c*e)/det
      u=(f*o-d*q)/det
      v=(a*q-c*o)/det
      w=(c*d-a*f)/det
      x=(d*p-o*e)/det
      y=(b*o-a*p)/det
      z=(a*e-b*d)/det
C
      XLXX= XLQSI*r*r+YLETA*s*s+ZLDZE*t*t
     &      +DEUX*(r*s*XY+s*t*YZ+r*t*XZ)
C
      XLXY=XLQSI*r*u+YLETA*s*v+ZLDZE*t*w
     &         +(s*u+r*v)*XY+(t*u+r*w)*XZ+(t*v+s*w)*YZ
C
      XLXZ=XLQSI*r*x+YLETA*s*y+ZLDZE*t*z
     &         +(s*x+r*y)*XY+(t*x+r*z)*XZ+(t*y+s*z)*YZ
C
      YLYY=XLQSI*u*u+YLETA*v*v+ZLDZE*w*w
     &         +DEUX*(u*v*XY+u*w*XZ+v*w*YZ)
C
      YLYZ=XLQSI*x*u+YLETA*y*v+ZLDZE*z*w
     &         +(x*v+y*u)*XY+(x*w+z*u)*XZ+(y*w+z*v)*YZ
C
      ZLZZ=XLQSI*x*x+YLETA*y*y+ZLDZE*z*z
     &         +DEUX*(x*y*XY+x*z*XZ+y*z*YZ)
C
      PXX=r*r+s*s+t*t
C
      PXY=r*u+s*v+t*w
C
      PXZ=r*x+s*y+t*z
C
      PYY=u*u+v*v+w*w
C
      PYZ=x*u+y*v+z*w
C
      PZZ=x*x+y*y+z*z
C
      VCOMP( 1)=XLXX
      VCOMP( 2)=YLYY
      VCOMP( 3)=ZLZZ
      VCOMP( 4)=XLXY
      VCOMP( 5)=XLXZ
      VCOMP( 6)=YLYZ
      VCOMP( 7)=PXX
      VCOMP( 8)=PYY
      VCOMP( 9)=PZZ
      VCOMP(10)=PXY
      VCOMP(11)=PXZ
      VCOMP(12)=PYZ
 
 666  CONTINUE
      IRET=1
      RETURN
 
C-------------------------
C   DETERMINANT NUL
C-------------------------
 999  CONTINUE
      CALL ERREUR(664)
      IRET=0
 
      RETURN
      END