Télécharger la plaquette Cast3M


fr en

Télécharger idmat4.eso

Retour à la liste

Numérotation des lignes : 

idmat4

  1. C IDMAT4    SOURCE    BP208322  17/03/01    21:17:40     9325           
  2.       SUBROUTINE IDMAT4 (NUMP1,NUMP2,NUDIR1,NUDIR2,ANG,ANG2,
  3.      .                   MELEME,MINTE,IPVAL,NPG2,MINTE2)
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8(A-H,O-Z)
  6. *
  7. ************************************************************************
  8. *
  9. *                             I D M A T 4
  10. *                             -----------
  11. *
  12. * FONCTION:
  13. * ---------
  14. *
  15. *     CALCUL DES COS-DIR. DES AXES ORTH. (OU ANISO.)PAR RAPPORT AU
  16. *
  17. *     REPERE LOCAL D'ELEMENT , A PARTIR D'UNE DEFINITION GLOBALE
  18. *
  19. *
  20. * MODULES UTILISES:
  21. * -----------------
  22. *
  23. -INC SMCOORD
  24.  
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. -INC CCREEL
  28. -INC SMELEME
  29. -INC SMINTE
  30. *
  31. * PARAMETRES:   (E)=ENTREE   (S)=SORTIE   (+ = CONTENU DANS UN COMMUN)
  32. * -----------
  33. *
  34. *     NUMP1   (E)  NUMERO DU POINT P1 ASSOCIE A NUDIR1
  35. *     NUMP2   (E)  NUMERO DU POINT P2 ASSOCIE A NUDIR1
  36. *     NUDIR1  (E)  NUMERO DE DIRECTIVE UTILISEE DANS LA LISTE:
  37. *                  "DIRECTION", "RADIAL"
  38. *     NUDIR2  (E)  NUMERO DE DIRECTIVE UTILISEE DANS LA LISTE:
  39. *                  "PARALLELE", "PERPENDIC.", "INCLINE",
  40. *                  POUR LA DEFINITION DES DIRECTIONS D'ORTHOTROPIE.
  41. *     ANG ANG2 (E) ANGLEs UTILISEs DANS LA DEFINITION DES DIRECTIONS
  42. *                  D ORTHOTROPIE (INCLINE)
  43. *     MELEME  (E)  POINTEUR DE "MAILLAGE" A 1 SEUL TYPE D'ELEMENT.
  44. *     MINTE   (E)  SEGMENT CONTENANT LES FONCTIONS DE FORME,ET LEURS
  45. *                  DERIVEES
  46. *     XVAL    (S)  DIRECTIONS D'ORTHOTROPIE  PAR ELEMENT,
  47. *                  ON FOURNIT LES COSINUS DIRECTEURS  DE L'AXE 1
  48. *                  EN BIDIM ,ET LES COSINUS DIRECTEURS DES AXES 1 ET
  49. *                  2 EN TRIDIM,PAR RAPPORT AU REPERE LOCAL DE L'ELEMENT
  50. *
  51. *
  52. * VARIABLES:
  53. *-----------
  54. *  VGLOB1,VGLOB2,VGLOB3  = COS.DIRECTEURS DES AXES 1, 2 ET 3 D'ORTH. PAR
  55. *   RAPPOT AU REPERE GLOBAL
  56. *  TXR   = TABLEAU CONTENANT LES COS.DIRECTEURS DES AXES LOCAUX PAR
  57. *           RAPPORT AU REPERE  GLOBAL
  58. *
  59. *
  60. * REMARQUES:
  61. * ----------
  62. *
  63. *  LES DIRECTION P1 ET P2 DEFINISSENT LE PLAN QUI CONTIENT LES AXES 1,2
  64. *  NUDIR2 =1 SIGNIFIE QUE L'AXE 1 EST PARALLELE A LA DIRECTION P1,
  65. *  NUDIR2 =2 SIGNIFIE QUE L'AXE 1 EST PERPENDICULAIRE A LA DIRECTION P1,
  66. *  NUDIR2 =3 SIGNIFIE QUE L'AXE 1 FAIT AN ANGLE DE (ANG) AVEC LA
  67. *  DIRECTION P1
  68. *
  69. *
  70. * AUTEUR, DATE DE CREATION:
  71. * -------------------------
  72. *
  73. *     P. DOWLATYARI OCT. 1990
  74. *
  75. * LANGAGE:
  76. * --------
  77. *       FORTRAN 77 + ESOPE
  78. *
  79. *
  80. ************************************************************************
  81. *
  82. *
  83.       DIMENSION VGLOB1(3),VGLOB2(3),VGLOB3(3),VAUX1(3),VAUX2(3)
  84.       DIMENSION VPT1(3),VPT2(3)
  85. *
  86.       SEGMENT YVAL
  87.          REAL*8 VLOC1(IDIM2,NPG2,NBELEM),VLOC2(IDIM2,NPG2,NBELEM)
  88.       ENDSEGMENT
  89. *
  90.       SEGMENT WORK
  91.          REAL*8 XE(3,NBNN),TXR(IDIM,IDIM)
  92.       ENDSEGMENT
  93. *
  94.       SEGACT MELEME
  95.       NBELEM=NUM(/2)
  96.       NBNN=NUM(/1)
  97.       SEGACT MINTE
  98.       IF(NUDIR1.EQ.2) SEGACT MINTE2
  99.  
  100.  
  101. ************************************************************************
  102. *     CALCUL DE VGLOB1   
  103. ************************************************************************
  104.  
  105.       VGLOB1(3)=0.D0
  106.       VGLOB2(3)=0.D0
  107.  
  108. **********      CAS 2D      *************
  109.  
  110.       IF (IDIM.EQ.2)THEN
  111. *
  112. *       COORDONNEES DU POINT P1
  113.         CALL EXCOO1 (NUMP1,VPT1(1),VPT1(2),VPT1(3),REEL1)
  114. *
  115. *       OPTION 'RADIAL' : VPT1 = CENTRE -> ON STOCKE PUIS GOTO 15
  116.         IF(NUDIR1.EQ.2) THEN
  117.           VAUX1(1)=VPT1(1)
  118.           VAUX1(2)=VPT1(2)
  119.           GO TO 15
  120.         ENDIF
  121. *
  122. *       OPTION 'DIRECTION' : VPT1 = DIRECTION 1
  123. *       NORMALISATION DU VECTEUR FOURNI 
  124. c       + CALCUL DE V2 = ORTHOGONAL A V1 DANS LE PLAN 2D
  125.         VNORM=SQRT(VPT1(1)*VPT1(1)+VPT1(2)*VPT1(2))
  126.         IF (VNORM.EQ.0.) THEN
  127.          CALL ERREUR (524)
  128.          SEGDES,MINTE,MELEME
  129.          RETURN
  130.         ENDIF
  131.         VGLOB1(1)=VPT1(1)/VNORM
  132.         VGLOB1(2)=VPT1(2)/VNORM
  133.         VGLOB2(1)=-VGLOB1(2)
  134.         VGLOB2(2)=VGLOB1(1)
  135. *
  136. *       OPTION PERPENDICULAIRE :
  137. *       ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR  DE L'AXE 3
  138.         IF(NUDIR2.EQ.2)THEN
  139.          VGLOB1(1)=VGLOB2(1)
  140.          VGLOB1(2)=VGLOB2(2)
  141.          VGLOB2(1)=-VGLOB1(2)
  142.          VGLOB2(2)=VGLOB1(1)
  143. *
  144. *       OPTION INCLINE :
  145. *       ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  146. *  -> \tilde{[V1,V2,V3]}=[eR eZ eT]*[cosa sina 0 ; -sina cosa 0 ; 0 0 1] 
  147. *       + eventuellement une rotation de ANG2 autour de \tilde{V1}
  148. *  -> [V1,V2,V3]=\tilde{[V1,V2,V3]}*[1 0 0; 0 cosb -sinb; 0 sinb cosb]
  149.         ELSEIF(NUDIR2.EQ.3)THEN        
  150.          COSA=COS(ANG)
  151.          SINA=SIN(ANG)
  152.          VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  153.          VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  154.          VGLOB2(1)=-VGLOB1(2)
  155.          VGLOB2(2)=VGLOB1(1)         
  156.          IF(ABS(ANG2).GT.XZPREC) THEN
  157. c          VGLOB3=(0 0 1)^T, VGLOB2=(XX1 XX2 0)^T
  158. c          VGLOB1 ne change pas, on tourne VGLOB2 (et VGLOB3 implicitement)
  159.            VGLOB2(1)=VGLOB2(1)*COS(ANG2)
  160.            VGLOB2(2)=VGLOB2(2)*COS(ANG2)
  161.            VGLOB2(3)=SIN(ANG2)
  162.          ENDIF
  163.  
  164.         ENDIF
  165.  
  166.  
  167. **********      CAS 3D      *************
  168.  
  169.       ELSEIF(IDIM.EQ.3)THEN
  170. *
  171. *       COORDONNEES DU POINT P1
  172.         CALL EXCOO1 (NUMP1,VPT1(1),VPT1(2),VPT1(3),REEL1)
  173. *
  174. *       OPTION 'DIRECTION' : VPT1 = DIRECTION 1
  175. *       NORMALISATION
  176.         IF(NUDIR1.NE.2) THEN
  177.           VNORM=SQRT(VPT1(1)*VPT1(1)+VPT1(2)*VPT1(2)+VPT1(3)*VPT1(3))
  178.           IF (VNORM.EQ.0.) THEN
  179.             CALL ERREUR (524)
  180.             SEGDES,MINTE,MELEME
  181.             RETURN
  182.           ENDIF
  183.           VGLOB1(1)=VPT1(1)/VNORM
  184.           VGLOB1(2)=VPT1(2)/VNORM
  185.           VGLOB1(3)=VPT1(3)/VNORM
  186.         ENDIF
  187. *
  188. *       COORDONNEES DU POINT P2
  189.         CALL EXCOO1 (NUMP2,VPT2(1),VPT2(2),VPT2(3),REEL1)
  190. *
  191. *       OPTION 'RADIAL' : DIRECTION 1 = P1P2 
  192. *       NORMALISATION PUIS GOTO 15
  193.         IF(NUDIR1.EQ.2) THEN
  194.          VGLOB1(1)=VPT2(1)-VPT1(1)
  195.          VGLOB1(2)=VPT2(2)-VPT1(2)
  196.          VGLOB1(3)=VPT2(3)-VPT1(3)
  197.          VNORM=SQRT(VGLOB1(1)*VGLOB1(1)+VGLOB1(2)*VGLOB1(2)+
  198.      .                VGLOB1(3)*VGLOB1(3))
  199.          VGLOB1(1)=VGLOB1(1)/VNORM
  200.          VGLOB1(2)=VGLOB1(2)/VNORM
  201.          VGLOB1(3)=VGLOB1(3)/VNORM
  202.          GO TO 15
  203.        ENDIF
  204. *
  205. *      OPTION 'DIRECTION' :
  206. *      DIRECTION 3 = PRODUIT VECTORIEL DES DEUX VECTEURS + NORMALISATION
  207.        CALL CROSS2 (VGLOB1(1),VPT2(1),VGLOB3(1),IRR)
  208.        IF (IRR.EQ.0)THEN
  209. *        DIRECTIONS D'ORTHOTROPIE MAL FOURNIES
  210.          CALL ERREUR (524)
  211.          SEGDES,MINTE,MELEME
  212.          IF(NUDIR1.EQ.2) SEGDES MINTE2
  213.          RETURN
  214.        ENDIF
  215. *      DIRECTION 2 = PRODUIT VECTORIEL DES DEUX VECTEURS + NORMALISATION
  216.        CALL CROSS2 (VGLOB3(1),VGLOB1(1),VGLOB2(1),IRR)
  217. *
  218. *      OPTION 'PERPENDICULAIRE'
  219. *      ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR  DE L'AXE 3
  220.        IF(NUDIR2.EQ.2)THEN
  221.          REEL1=VGLOB1(1)
  222.          VGLOB1(1)=VGLOB2(1)
  223.          VGLOB2(1)=-REEL1
  224.          REEL1=VGLOB1(2)
  225.          VGLOB1(2)=VGLOB2(2)
  226.          VGLOB2(2)=-REEL1
  227.          REEL1=VGLOB1(3)
  228.          VGLOB1(3)=VGLOB2(3)
  229.          VGLOB2(3)=-REEL1
  230. *
  231. *      OPTION 'INCLINE'
  232. *      ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  233.        ELSEIF(NUDIR2.EQ.3)THEN
  234.          COSA=COS(ANG)
  235.          SINA=SIN(ANG)
  236.          REEL1=VGLOB1(1)
  237.          VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  238.          VGLOB2(1)=-REEL1*SINA+VGLOB2(1)*COSA
  239.          REEL1=VGLOB1(2)
  240.          VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  241.          VGLOB2(2)=-REEL1*SINA+VGLOB2(2)*COSA
  242.          REEL1=VGLOB1(3)
  243.          VGLOB1(3)=VGLOB1(3)*COSA+VGLOB2(3)*SINA
  244.          VGLOB2(3)=-REEL1*SINA+VGLOB2(3)*COSA
  245.         ENDIF
  246.  
  247.        ENDIF
  248. ************************************************************************
  249.  
  250. *
  251. 15    CONTINUE
  252. *
  253.       IF(IFOUR.EQ.1)THEN
  254.         IDIM2=3
  255.       ELSE
  256.         IDIM2=IDIM
  257.       ENDIF
  258.       SEGINI YVAL
  259.       IPVAL=YVAL
  260.       SEGINI WORK
  261. *
  262. ************************************************************************
  263. *     BOUCLE SUR LES ELEMENTS
  264. ************************************************************************
  265. *
  266.       DO 10  IEL=1,NBELEM
  267. *
  268.       DO 20  NC=1,IDIM2
  269.       DO 20  NV=1,NPG2
  270.        VLOC1(NC,NV,IEL)=0.D0
  271.        VLOC2(NC,NV,IEL)=0.D0
  272.  20   CONTINUE
  273. *
  274. *     COORDONNEES DES NOEUDS
  275. *
  276.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  277. *
  278. *     CALCUL DU REPERE LOCAL DE L'ELEMENT
  279. *
  280.       NBSH=SHPTOT(/2)
  281.       CALL RLOCAL(XE,SHPTOT,NBSH,NBNN,TXR)
  282.         if (nbsh.eq.-1) then
  283.          call erreur(525)
  284.          return
  285.         endif
  286.       IF(IERR.NE.0)THEN
  287.        SEGSUP,WORK,YVAL
  288.        SEGDES , MELEME,MINTE
  289.        IF(NUDIR1.EQ.2) SEGDES MINTE2
  290.        RETURN
  291.       ENDIF
  292. *
  293. ********************************************************
  294. *     BOUCLE SUR LES POINTS DE GAUSS
  295. ********************************************************
  296. *
  297.       DO 80 IGAU=1,NPG2
  298.  
  299. c---- OPTION 'RADIAL' -----
  300.       IF (NUDIR1.EQ.2) THEN
  301.  
  302.         DO 25 ICO=1,IDIM
  303.           VAUX2(ICO)=0.D0
  304.           DO 25 ID=1,NBNN
  305.           VAUX2(ICO)=VAUX2(ICO)+
  306.      .          MINTE2.SHPTOT(1,ID,IGAU)*XE(ICO,ID)
  307.  25     CONTINUE
  308. *       
  309. *     --CAS 2D-- 
  310.         IF(IDIM.EQ.2) THEN
  311.  
  312.           VGLOB1(1)=VAUX1(1)-VAUX2(1)
  313.           VGLOB1(2)=VAUX1(2)-VAUX2(2)
  314. *         NORMALISATION
  315.           VNORM=SQRT(VGLOB1(1)*VGLOB1(1)+VGLOB1(2)*VGLOB1(2))
  316.           VGLOB1(1)=VGLOB1(1)/VNORM
  317.           VGLOB1(2)=VGLOB1(2)/VNORM
  318.           VGLOB2(1)=-VGLOB1(2)
  319.           VGLOB2(2)=VGLOB1(1)
  320.  
  321. *         OPTION 'PERPENDICULAIRE'
  322. *         ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR  DE L'AXE 3
  323.           IF(NUDIR2.EQ.2)THEN
  324.             VGLOB1(1)=VGLOB2(1)
  325.             VGLOB1(2)=VGLOB2(2)
  326.  
  327. *         OPTION 'INCLINE'
  328. *         ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  329.           ELSEIF(NUDIR2.EQ.3)THEN
  330.             COSA=COS(ANG)
  331.             SINA=SIN(ANG)
  332.             VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  333.             VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  334.           ENDIF
  335.  
  336.           IF(IFOUR.EQ.1) THEN
  337.             VGLOB1(3)=0.D0
  338.             VGLOB2(1)=-VGLOB1(2)
  339.             VGLOB2(2)=VGLOB1(1)
  340.             VGLOB2(3)=0.D0
  341.           ENDIF
  342.  
  343.         ENDIF
  344. *     ----------
  345.  
  346. *     --CAS 3D-- 
  347.         IF(IDIM.EQ.3) THEN
  348.  
  349.           VAUX1(1)=VAUX2(1)-VPT2(1)
  350.           VAUX1(2)=VAUX2(2)-VPT2(2)
  351.           VAUX1(3)=VAUX2(3)-VPT2(3)
  352.           PSC=VAUX1(1)*VGLOB1(1)+VAUX1(2)*VGLOB1(2)+VAUX1(3)*VGLOB1(3)
  353.           VGLOB2(1)=VAUX1(1)-PSC*VGLOB1(1)
  354.           VGLOB2(2)=VAUX1(2)-PSC*VGLOB1(2)
  355.           VGLOB2(3)=VAUX1(3)-PSC*VGLOB1(3)
  356. *         NORMALISATION
  357.           VNORM=SQRT(VGLOB2(1)*VGLOB2(1)+VGLOB2(2)*VGLOB2(2)+
  358.      .                VGLOB2(3)*VGLOB2(3))
  359.           VGLOB2(1)=VGLOB2(1)/VNORM
  360.           VGLOB2(2)=VGLOB2(2)/VNORM
  361.           VGLOB2(3)=VGLOB2(3)/VNORM
  362. *
  363. *         PRODUIT VECTORIEL   DES DEUX VECTEURS  ET NORMALISATION
  364.           CALL CROSS2 (VGLOB1(1),VGLOB2(1),VGLOB3(1),IRR)
  365.           IF (IRR.EQ.0)THEN
  366.             CALL ERREUR (524)
  367.             SEGDES,MINTE,MELEME
  368.             IF(NUDIR1.EQ.2) SEGDES MINTE2
  369.             RETURN
  370.           ENDIF
  371. *
  372. *         OPTION 'PERPENDICULAIRE'
  373. *         ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR  DE L'AXE 3
  374.           IF(NUDIR2.EQ.2)THEN
  375.             REEL1=VGLOB1(1)
  376.             VGLOB1(1)=VGLOB2(1)
  377.             VGLOB2(1)=-REEL1
  378.             REEL1=VGLOB1(2)
  379.             VGLOB1(2)=VGLOB2(2)
  380.             VGLOB2(2)=-REEL1
  381.             REEL1=VGLOB1(3)
  382.             VGLOB1(3)=VGLOB2(3)
  383.             VGLOB2(3)=-REEL1
  384. *
  385. *         OPTION 'INCLINE'
  386. *         ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  387.           ELSEIF(NUDIR2.EQ.3)THEN
  388.             COSA=COS(ANG)
  389.             SINA=SIN(ANG)
  390.             REEL1=VGLOB1(1)
  391.             VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  392.             VGLOB2(1)=-REEL1*SINA+VGLOB2(1)*COSA
  393.             REEL1=VGLOB1(2)
  394.             VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  395.             VGLOB2(2)=-REEL1*SINA+VGLOB2(2)*COSA
  396.             REEL1=VGLOB1(3)
  397.             VGLOB1(3)=VGLOB1(3)*COSA+VGLOB2(3)*SINA
  398.             VGLOB2(3)=-REEL1*SINA+VGLOB2(3)*COSA
  399.           ENDIF
  400.  
  401.         ENDIF
  402. *     ----------
  403.  
  404.       ENDIF
  405. c---- FIN DE L'OPTION 'RADIAL' -----
  406. *
  407. *     CALCUL DES COS.DIRECTEURS DU REPERE LOCAL / REPERE GLOBAL
  408. *
  409.       IF(IDIM.EQ.2)THEN
  410.         VLOC1(1,IGAU,IEL)=TXR(1,1)*VGLOB1(1)+TXR(2,1)*VGLOB1(2)
  411.         VLOC1(2,IGAU,IEL)=TXR(1,2)*VGLOB1(1)+TXR(2,2)*VGLOB1(2)
  412.         IF(IFOUR.EQ.1)THEN
  413. c       bp: en 2D Fourier, TXR= [ TXR(1,1) TXR(1,2) 0 ]  dans (R,Z,T)
  414. c                               [ TXR(2,1) TXR(2,2) 0 ]
  415. c                               [   0        0      1 ]
  416.           VLOC1(3,IGAU,IEL)=VGLOB1(3)
  417.           VLOC2(1,IGAU,IEL)=TXR(1,1)*VGLOB2(1)+TXR(2,1)*VGLOB2(2)
  418.           VLOC2(2,IGAU,IEL)=TXR(1,2)*VGLOB2(1)+TXR(2,2)*VGLOB2(2)
  419.           VLOC2(3,IGAU,IEL)=VGLOB2(3)
  420.         ENDIF
  421. *
  422.       ELSEIF(IDIM.EQ.3)THEN
  423.          DO  40 J=1,3
  424.          DO  40 I=1,3
  425.            VLOC1(J,IGAU,IEL)=VLOC1(J,IGAU,IEL)+VGLOB1(I)*TXR(I,J)
  426.            VLOC2(J,IGAU,IEL)=VLOC2(J,IGAU,IEL)+VGLOB2(I)*TXR(I,J)
  427.  40      CONTINUE
  428.       ENDIF
  429.  
  430.  80   CONTINUE
  431. ********************************************************
  432. *     FIN DE BOUCLE SUR LES POINTS DE GAUSS
  433. ********************************************************
  434.  
  435.  10   CONTINUE
  436. ************************************************************************
  437. *     FIN DE BOUCLE SUR LES ELEMENTS
  438. ************************************************************************
  439. *
  440. *     DESACTIVATION DES SEGMENTS
  441. *
  442.       SEGSUP,WORK
  443.       SEGDES , MELEME,MINTE,YVAL
  444.       IF(NUDIR1.EQ.2) SEGDES MINTE2
  445. *
  446.       RETURN
  447.       END
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  

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