Sources Esope | Cast3M

idmat5
C IDMAT5    SOURCE    BP208322  16/12/14    21:15:15     9258           
      SUBROUTINE IDMAT5 (NUDIR1,NUMP1,NUDIR2,NUMP3,ANG,MELEME,MINTE,
     1                   MINTE1,IPEPAI,IPEXCE,NPG2,ITHER,IPVAL)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
*
************************************************************************
*
*                             I D M A T 5
*                             -----------
*
* FONCTION:
* ---------
*
*     CALCUL DES COSINUS-DIRECTEURS DU 1ER AXE D'ORTHOTROPIE /REPERES
*
*     LOCAUX D'UNE COQUE EPAISSE
*
*
* MODULES UTILISES:
* -----------------
*
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC CCREEL
-INC SMELEME
-INC SMINTE
-INC SMCHAML
*
* PARAMETRES:   (E)=ENTREE   (S)=SORTIE   (+ = CONTENU DANS UN COMMUN)
* -----------
*
*     NUDIR1  (E)  NUMERO DE DIRECTIVE UTILISEE DANS LA LISTE:
*                  "DIRECTION", "RADIAL",
*                  POUR LA DEFINITION DES DIRECTIONS D'ORTHOTROPIE.
*     NUMP1   (E)  NUMERO DU POINT P ASSOCIE A NUDIR1
*     NUDIR2  (E)  NUMERO DE DIRECTIVE UTILISEE DANS LA LISTE:
*                  "PARALLELE", "PERPENDIC.", "INCLINE",
*                  POUR LA DEFINITION DES DIRECTIONS D'ORTHOTROPIE.
*     NUMP3   (E)  NUMERO DU POINT P ASSOCIE A NUDIR2
*     ANG     (E)  ANGLE UTILISE DANS LA DEFINITION DES DIRECTIONS
*                  D ORTHOTROPIE (INCLINE)
*     MELEME  (E)  POINTEUR DE "MAILLAGE" A 1 SEUL TYPE D'ELEMENT.
*     MINTE   (E)  SEGMENT CONTENANT LES FONCTIONS DE FORME ET LEURS
*                  DERIVEES AU POINTS DE GAUSS
*     MINTE1  (E)  SEGMENT CONTENANT LES FONCTIONS DE FORME ET LEURS
*                  DERIVEES AU NOEUDS
*     IPEPAI  (E)  POINTEUR DU SEGMENT CONTENANT LES EPAISSEURS
*     IPEXCE  (E)  POINTEUR DU SEGMENT CONTENANT LES EXCENTREMENTS
*     XVAL    (S)  DIRECTIONS D'ORTHOTROPIE PAR POINT DE GAUSS
*                  ON FOURNIT LE COSINUS ET LE SINUS DE L'ANGLE QUE FAIT
*                  LA DIRECTION 1 D'ORTHOTROPIE AVEC L'AXE 1 DU REPRE
*                  LIE  AU POINT DE GAUSS
*     NPG2    (E)  NOMBRE DE POINTS DE GAUSS
*
*
* VARIABLES:
* ----------
*
*     COSA     = COSINUS DE L'ANGLE DE LA DIRECTION 1 D'ORTHOTROPIE AVEC
*               L'AXE 1 DU REPERE LOCAL SITUE AU POINT DE GAUSS CONSIDER
*     SINA     = SINUS DE L'ANGLE DE LA DIRECTION 1 D'ORTHOTROPIE AVEC L
*               'AXE 1 DU REPERE LOCAL SITUE AU POINT DE GAUSS CONSIDERE
*     XI,YI,ZI = COSINUS-DIRECTEURS DE L'AXE 1 DU REPERE LOCAL SITUE AU
*                POINT DE GAUSS CONSIDERE
*     XN,YN,ZN = COSINUS-DIRECTEURS DE L'AXE 3 DU REPERE LOCAL SITUE AU
*                 POINT DE GAUSS CONSIDERE
*     XDIR1,.. = COORDONNEES DE LA DIRECTION FOURNIE POUR DETERMINER LES
*                DIRECTIONS D'ORTHOTROPIE.
*
* REMARQUES:
* ----------
*
*  LA PREMIERE DIRECTION D'ORTHOTROPIE  EST  DEFINIE A CHAQUE POINT DE
*  GAUSS PAR RAPPORT A UN REPERE LOCAL DONT LES DEUX PREMIERS AXES SONT
*  SITUES DANS LE PLAN TANGENT
*
* AUTEUR, DATE DE CREATION:
* -------------------------
*
*     P DOWLATYARI       JANVIER 1991
*
* LANGAGE:
* --------
*
*     ESOPE77 FORTRAN77 + EXTENSION
*
*
************************************************************************
*
*
      REAL*8 XJ(3,3),TT(9)
*
      SEGMENT GEOM
         REAL*8 EP(NBNN),TXR(3,3,NBNN),EXC(NBNN),XE(3,NBNN)
      ENDSEGMENT
*
      SEGMENT XVAL
         REAL*8 CVAL(NPG2,NBELEM),SVAL(NPG2,NBELEM)
      ENDSEGMENT
*
*
* CONSTANTES:
* -----------
*
*     COS1D  = VALEUR APPROCHEE DE  COSINUS(1 DEGRE).
*     SIN1D  = VALEUR APPROCHEE DE  SINUS(1 DEGRE).
*
      PARAMETER (COS1D = 0.999 85 D0)
      PARAMETER (SIN1D = 0.017 45 D0)
      PARAMETER (UNDEMI= 0.5D0  )
*
         SEGACT MELEME
         SEGACT MINTE
         SEGACT MINTE1
         NBELEM=NUM(/2)
         NBNN  = NUM(/1)
*
*        DIRECTION FOURNIE, AVANT PROJECTION SUR LA COQUE.
*
*        ON CHERCHE LE POINT NUMP1 RATACHER A LA DIRECTIVE DIRECTION
*
         CALL EXCOO1 (NUMP1,XDIR1,YDIR1,ZDIR1,REEL1)
         IF (NUDIR1 .EQ. 1) THEN
            XNORM = SQRT(XDIR1**2 + YDIR1**2 + ZDIR1**2)
            XDIR1 = XDIR1 / XNORM
            YDIR1 = YDIR1 / XNORM
            ZDIR1 = ZDIR1 / XNORM
         END IF
         XDIR2 = XDIR1
         YDIR2 = YDIR1
         ZDIR2 = ZDIR1
*
         SEGINI XVAL
         IPVAL=XVAL
         SEGINI GEOM
*
*
*
         DO 10 IEL=1,NBELEM
*
*        ON CHERCHE LES COORDONNEES DES NOEUDS DE L'ELEMENT IEL,
*        DANS LE REPERE GLOBAL
*
         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
*
*
*       CALCUL DES AXES LOCAUX A TOUS LES NOEUDS DE L'ELEMENT
*
*
         CALL CQ8LOC (XE,NBNN,MINTE1.SHPTOT,TXR,IRR)
         IF(IRR.EQ.0)THEN
*
*        ECHEC DANS LE CALCUL DES AXES LOCAUX
*
        CALL ERREUR (515)
         SEGSUP ,GEOM,XVAL
         SEGDES,MELEME,MINTE,MINTE1
         RETURN
         ENDIF
*
*        ON CHERCHE LES EPAISSEURS ET LES EXCENTREMENTS
           DO 20  IG=1,NBNN
             MELVAL=IPEPAI
             IGMN=MIN(IG,VELCHE(/1))
             IBMN=MIN(IEL,VELCHE(/2))
             EP(IG)=VELCHE(IGMN,IBMN)
             IF(IPEXCE.NE.0)THEN
               MELVAL=IPEXCE
               IGMN=MIN(IG,VELCHE(/1))
               IBMN=MIN(IEL,VELCHE(/2))
               EXC(IG)=VELCHE(IGMN,IBMN)
           ELSE
              EXC(IG)=XZERO
             ENDIF
 20        CONTINUE
         NPG1=NPG2/2
*
*          BOUCLE SUR LES POINTS DE GAUSS
*
         DO  50  IGAU = 1,NPG2
*
*
*        CALCUL DU JACOBIEN ET DE SON DETERMINENT EN CE POINT DE GAUSS
*
         E3=DZEGAU(IGAU)
*
         CALL CQ8JCE(IGAU,NBNN,E3,XE,EP,EXC,TXR,SHPTOT,XJ,DJAC,IRR)
         IF (IRR.LT.0)THEN
*        JACOBIEN NUL DANS L'ELEMENT IEL
         INTERR(1)=IEL
         CALL ERREUR (405)
          SEGSUP ,GEOM,XVAL
          SEGDES,MELEME,MINTE,MINTE1
          RETURN
         ENDIF
*
*    DETERMINATION DES COSINUS DIRECTEURS DES AXES LOCAUX EN CE POINT
*
*
         DO 30 I=1,3
         DO 30 J=1,2
          K=3*(J-1)+I
          TT(K)    = XJ(J,I)
   30    CONTINUE
*
*    PRODUITS VECTORIELS ET NORMALISATIONS
*
         CALL CROSS2(TT(1),TT(4),TT(7),IRR)
         CALL CROSS2(TT(7),TT(1),TT(4),IRR)
         CALL CROSS2(TT(4),TT(7),TT(1),IRR)
*
        IF (NUDIR1 .EQ. 1) THEN
             XDIR1 = XDIR2
             YDIR1 = YDIR2
             ZDIR1 = ZDIR2
           ELSE
*           OPTION "RADIAL": L'ORIENTATION DE L'ORTHOTROPIE EST
*           DETERMINEE A PARTIR DU VECTEUR DE VALEUR (PTGAUSS - P1
             XDIR1 = - XDIR2
             YDIR1 = - YDIR2
             ZDIR1 = - ZDIR2
             DO 40 INOE=1,NBNN
        XDIR1 = XDIR1 + SHPTOT(1,INOE,IGAU)*(XE(1,INOE)+
     1        UNDEMI*EP(INOE)*E3*TXR(1,3,INOE)+EXC(INOE)*TXR(1,3,INOE))
        YDIR1 = YDIR1 + SHPTOT(1,INOE,IGAU)*(XE(2,INOE)+
     1        UNDEMI*EP(INOE)*E3*TXR(2,3,INOE)+EXC(INOE)*TXR(2,3,INOE))
        ZDIR1 = ZDIR1 + SHPTOT(1,INOE,IGAU)*(XE(3,INOE)+
     1        UNDEMI*EP(INOE)*E3*TXR(3,3,INOE)+EXC(INOE)*TXR(3,3,INOE))
  40        CONTINUE
            XNORM = SQRT(XDIR1**2 + YDIR1**2 + ZDIR1**2)
            XDIR1 = XDIR1 / XNORM
            YDIR1 = YDIR1 / XNORM
            ZDIR1 = ZDIR1 / XNORM
        END IF
        XI=TT(1)
        YI=TT(2)
        ZI=TT(3)
        XN=TT(7)
        YN=TT(8)
        ZN=TT(9)
*
*     PROJECTION SUR LE PLAN TANGENT LE VECTEUR D'ORIENTATION DE
*     L'ORTROTHOPIE.
*
            P1N = XDIR1*XN + YDIR1*YN + ZDIR1*ZN
            IF (ABS(P1N) .GT. COS1D) THEN
*               DIRECTION FOURNIE ORTHOGONALE A LA COQUE.
                  CALL ERREUR (344)
                  SEGSUP ,GEOM,XVAL
                  SEGDES,MELEME,MINTE,MINTE1
                  RETURN
             END IF
             XDIR1 = XDIR1 - P1N*XN
             YDIR1 = YDIR1 - P1N*YN
             ZDIR1 = ZDIR1 - P1N*ZN
             XNORM = SQRT(XDIR1**2 + YDIR1**2 + ZDIR1**2)
             XDIR1 = XDIR1 / XNORM
             YDIR1 = YDIR1 / XNORM
             ZDIR1 = ZDIR1 / XNORM
*
*              ANGLE DU VECTEUR PROJETE AVEC LA DIRECTION 1 DE LA COQUE:
             COSA = XI * XDIR1 + YI * YDIR1 + ZI * ZDIR1
             SINA = XN * (YI*ZDIR1 - ZI*YDIR1)
     &              + YN * (ZI*XDIR1 - XI*ZDIR1)
     &              + ZN * (XI*YDIR1 - YI*XDIR1)
*
*
             IF (NUDIR2 .EQ. 2) THEN
*                 'PERPENDICULAIRE' -> ON AJOUTE PI/2 A L'ANGLE:
                  REEL1 = COSA
                  COSA = - SINA
                  SINA = REEL1
             ELSE IF (NUDIR2 .EQ. 3) THEN
*                 'INCLINE' -> ON TOURNE de +/-ANG
                  COSB = COS(ANG)
                  SINB = SIN(ANG)
c                 si pas de point P3, on oriente selon la normale a la coque
                  if (NUMP3.eq.0) then
                    REEL1=1.D0
                  else
                     CALL EXCOO1 (NUMP3,XDIR3,YDIR3,ZDIR3,REEL1)
                     XNORM = SQRT(XDIR3**2 + YDIR3**2 + ZDIR3**2)
                     IF(XNORM .LT. XPETIT) THEN
C                       POINT P3 est trop proche de l'origine
                        CALL ERREUR (239)
                        RETURN
                     ENDIF
                    XDIR3 = XDIR3 / XNORM
                    YDIR3 = YDIR3 / XNORM
                    ZDIR3 = ZDIR3 / XNORM         
                    REEL1 = XN * XDIR3 + YN * YDIR3 + ZN * ZDIR3
                  endif
                  IF (REEL1 .GT. XPETIT) THEN
*                    prod scal >0 -> ON AJOUTE "ANGLE":
                     REEL1 = SINA*COSB + COSA*SINB
                     COSA = COSA*COSB - SINA*SINB
                     SINA = REEL1
                  ELSE IF (REEL1 .LT. (-XPETIT) ) THEN
*                    prod scal >0 -> ON RETRANCHE "ANGLE":
                     REEL1 = SINA*COSB - COSA*SINB
                     COSA = COSA*COSB + SINA*SINB
                     SINA = REEL1
                  ELSE
*                    POINT P3 DONNE POUR ORIENTATION INAPROPRIE.
                     CALL ERREUR (278)
                     SEGSUP ,GEOM,XVAL
                     SEGDES,MELEME,MINTE,MINTE1
                     RETURN
                  END IF
            END IF
*
            CVAL(IGAU,IEL) = COSA
            SVAL(IGAU,IEL) = SINA
*
  50        CONTINUE
*
  10   CONTINUE
*
*
        SEGSUP ,GEOM
        SEGDES,XVAL
        SEGDES,MELEME,MINTE,MINTE1
*
      RETURN
      END