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    CB215821  26/08/24    21:16:48     12622          
  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 81  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.    81 CONTINUE
  274. *
  275. *     COORDONNEES DES NOEUDS
  276. *
  277.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  278. *
  279. *     CALCUL DU REPERE LOCAL DE L'ELEMENT
  280. *
  281.       NBSH=SHPTOT(/2)
  282.       CALL RLOCAL(XE,SHPTOT,NBSH,NBNN,TXR)
  283.         if (nbsh.eq.-1) then
  284.          call erreur(525)
  285.          return
  286.         endif
  287.       IF(IERR.NE.0)THEN
  288.        SEGSUP,WORK,YVAL
  289.        SEGDES , MELEME,MINTE
  290.        IF(NUDIR1.EQ.2) SEGDES MINTE2
  291.        RETURN
  292.       ENDIF
  293. *
  294. ********************************************************
  295. *     BOUCLE SUR LES POINTS DE GAUSS
  296. ********************************************************
  297. *
  298.       DO 80 IGAU=1,NPG2
  299.  
  300. c---- OPTION 'RADIAL' -----
  301.       IF (NUDIR1.EQ.2) THEN
  302.  
  303.         DO 82 ICO=1,IDIM
  304.           VAUX2(ICO)=0.D0
  305.           DO 25 ID=1,NBNN
  306.           VAUX2(ICO)=VAUX2(ICO)+
  307.      .          MINTE2.SHPTOT(1,ID,IGAU)*XE(ICO,ID)
  308.  25     CONTINUE
  309.    82   CONTINUE
  310. *       
  311. *     --CAS 2D-- 
  312.         IF(IDIM.EQ.2) THEN
  313.  
  314.           VGLOB1(1)=VAUX1(1)-VAUX2(1)
  315.           VGLOB1(2)=VAUX1(2)-VAUX2(2)
  316. *         NORMALISATION
  317.           VNORM=SQRT(VGLOB1(1)*VGLOB1(1)+VGLOB1(2)*VGLOB1(2))
  318.           VGLOB1(1)=VGLOB1(1)/VNORM
  319.           VGLOB1(2)=VGLOB1(2)/VNORM
  320.           VGLOB2(1)=-VGLOB1(2)
  321.           VGLOB2(2)=VGLOB1(1)
  322.  
  323. *         OPTION 'PERPENDICULAIRE'
  324. *         ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR  DE L'AXE 3
  325.           IF(NUDIR2.EQ.2)THEN
  326.             VGLOB1(1)=VGLOB2(1)
  327.             VGLOB1(2)=VGLOB2(2)
  328.  
  329. *         OPTION 'INCLINE'
  330. *         ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  331.           ELSEIF(NUDIR2.EQ.3)THEN
  332.             COSA=COS(ANG)
  333.             SINA=SIN(ANG)
  334.             VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  335.             VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  336.           ENDIF
  337.  
  338.           IF(IFOUR.EQ.1) THEN
  339.             VGLOB1(3)=0.D0
  340.             VGLOB2(1)=-VGLOB1(2)
  341.             VGLOB2(2)=VGLOB1(1)
  342.             VGLOB2(3)=0.D0
  343.           ENDIF
  344.  
  345.         ENDIF
  346. *     ----------
  347.  
  348. *     --CAS 3D-- 
  349.         IF(IDIM.EQ.3) THEN
  350.  
  351.           VAUX1(1)=VAUX2(1)-VPT2(1)
  352.           VAUX1(2)=VAUX2(2)-VPT2(2)
  353.           VAUX1(3)=VAUX2(3)-VPT2(3)
  354.           PSC=VAUX1(1)*VGLOB1(1)+VAUX1(2)*VGLOB1(2)+VAUX1(3)*VGLOB1(3)
  355.           VGLOB2(1)=VAUX1(1)-PSC*VGLOB1(1)
  356.           VGLOB2(2)=VAUX1(2)-PSC*VGLOB1(2)
  357.           VGLOB2(3)=VAUX1(3)-PSC*VGLOB1(3)
  358. *         NORMALISATION
  359.           VNORM=SQRT(VGLOB2(1)*VGLOB2(1)+VGLOB2(2)*VGLOB2(2)+
  360.      .                VGLOB2(3)*VGLOB2(3))
  361.           VGLOB2(1)=VGLOB2(1)/VNORM
  362.           VGLOB2(2)=VGLOB2(2)/VNORM
  363.           VGLOB2(3)=VGLOB2(3)/VNORM
  364. *
  365. *         PRODUIT VECTORIEL   DES DEUX VECTEURS  ET NORMALISATION
  366.           CALL CROSS2 (VGLOB1(1),VGLOB2(1),VGLOB3(1),IRR)
  367.           IF (IRR.EQ.0)THEN
  368.             CALL ERREUR (524)
  369.             SEGDES,MINTE,MELEME
  370.             IF(NUDIR1.EQ.2) SEGDES MINTE2
  371.             RETURN
  372.           ENDIF
  373. *
  374. *         OPTION 'PERPENDICULAIRE'
  375. *         ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR  DE L'AXE 3
  376.           IF(NUDIR2.EQ.2)THEN
  377.             REEL1=VGLOB1(1)
  378.             VGLOB1(1)=VGLOB2(1)
  379.             VGLOB2(1)=-REEL1
  380.             REEL1=VGLOB1(2)
  381.             VGLOB1(2)=VGLOB2(2)
  382.             VGLOB2(2)=-REEL1
  383.             REEL1=VGLOB1(3)
  384.             VGLOB1(3)=VGLOB2(3)
  385.             VGLOB2(3)=-REEL1
  386. *
  387. *         OPTION 'INCLINE'
  388. *         ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  389.           ELSEIF(NUDIR2.EQ.3)THEN
  390.             COSA=COS(ANG)
  391.             SINA=SIN(ANG)
  392.             REEL1=VGLOB1(1)
  393.             VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  394.             VGLOB2(1)=-REEL1*SINA+VGLOB2(1)*COSA
  395.             REEL1=VGLOB1(2)
  396.             VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  397.             VGLOB2(2)=-REEL1*SINA+VGLOB2(2)*COSA
  398.             REEL1=VGLOB1(3)
  399.             VGLOB1(3)=VGLOB1(3)*COSA+VGLOB2(3)*SINA
  400.             VGLOB2(3)=-REEL1*SINA+VGLOB2(3)*COSA
  401.           ENDIF
  402.  
  403.         ENDIF
  404. *     ----------
  405.  
  406.       ENDIF
  407. c---- FIN DE L'OPTION 'RADIAL' -----
  408. *
  409. *     CALCUL DES COS.DIRECTEURS DU REPERE LOCAL / REPERE GLOBAL
  410. *
  411.       IF(IDIM.EQ.2)THEN
  412.         VLOC1(1,IGAU,IEL)=TXR(1,1)*VGLOB1(1)+TXR(2,1)*VGLOB1(2)
  413.         VLOC1(2,IGAU,IEL)=TXR(1,2)*VGLOB1(1)+TXR(2,2)*VGLOB1(2)
  414.         IF(IFOUR.EQ.1)THEN
  415. c       bp: en 2D Fourier, TXR= [ TXR(1,1) TXR(1,2) 0 ]  dans (R,Z,T)
  416. c                               [ TXR(2,1) TXR(2,2) 0 ]
  417. c                               [   0        0      1 ]
  418.           VLOC1(3,IGAU,IEL)=VGLOB1(3)
  419.           VLOC2(1,IGAU,IEL)=TXR(1,1)*VGLOB2(1)+TXR(2,1)*VGLOB2(2)
  420.           VLOC2(2,IGAU,IEL)=TXR(1,2)*VGLOB2(1)+TXR(2,2)*VGLOB2(2)
  421.           VLOC2(3,IGAU,IEL)=VGLOB2(3)
  422.         ENDIF
  423. *
  424.       ELSEIF(IDIM.EQ.3)THEN
  425.          DO  83 J=1,3
  426.          DO  40 I=1,3
  427.            VLOC1(J,IGAU,IEL)=VLOC1(J,IGAU,IEL)+VGLOB1(I)*TXR(I,J)
  428.            VLOC2(J,IGAU,IEL)=VLOC2(J,IGAU,IEL)+VGLOB2(I)*TXR(I,J)
  429.  40      CONTINUE
  430.    83    CONTINUE
  431.       ENDIF
  432.  
  433.  80   CONTINUE
  434. ********************************************************
  435. *     FIN DE BOUCLE SUR LES POINTS DE GAUSS
  436. ********************************************************
  437.  
  438.  10   CONTINUE
  439. ************************************************************************
  440. *     FIN DE BOUCLE SUR LES ELEMENTS
  441. ************************************************************************
  442. *
  443. *     DESACTIVATION DES SEGMENTS
  444. *
  445.       SEGSUP,WORK
  446.       SEGDES , MELEME,MINTE,YVAL
  447.       IF(NUDIR1.EQ.2) SEGDES MINTE2
  448. *
  449.       RETURN
  450.       END
  451.  
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  

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