Télécharger la plaquette Cast3M


fr en

Télécharger rtens1.eso

Retour à la liste

Numérotation des lignes : 

rtens1

  1. C RTENS1    SOURCE    OF166741  25/02/21    21:18:24     12166          
  2.       SUBROUTINE RTENS1(IPCHE1,IFOMEM,IMOT,IPTV2,IELEME,IVAVEC,IVACOM,
  3.      &                  IVARES,IDEFO,IINTE,MELE,NPINT,NVEC,V1,V2,W2,W3,
  4.      &                  CENTR1,CENTR2,AXEI1,IER1)
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8(A-H,O-Z)
  7. *-----------------------------------------------------------------------*
  8. *     Operateur RTENS : cas de la formulation massive                   *
  9. *                                                                       *
  10. *     IPCHE1    (e)    pointeur sur un MCHAML de caracteristiques       *
  11. *                      = 0 si isotropie                                 *
  12. *     IFOMEM    (e)    = IFOUR de CCOPTIO                               *
  13. *     IMOT      (e)    indique le type de repere desire (cf RTENS)      *
  14. *     IPTV2     (e)    pointeur sur le 2nd point repere                 *
  15. *     IELEME    (e)    pointeur sur le segment MELEME (actif)           *
  16. *     IVAVEC    (e/s)  pointeur sur un segment MPTVAL (actif)           *
  17. *     IVACOM    (e/s)  pointeur sur un segment MPTVAL (actif)           *
  18. *     IVARES    (e/s)  pointeur sur un segment MPTVAL (actif)           *
  19. *     IDEFO     (e)    =1 : tenseur de deformations (contraintes sinon) *
  20. *     IINTE     (e)    pointeur sur le segment MINTE  (actif)           *
  21. *     MELE      (e)    numero de l'element-fini dans NOMTP              *
  22. *     NPINT     (e)    nombre de points d'integration (coques)          *
  23. *     NVEC      (e)    nombre de composantes du futur MCHAML            *
  24. *     V1        (e)    coordonnees et norme du 1er vecteur              *
  25. *     V2        (e)    coordonnees et norme du 2nd vecteur              *
  26. *     W2        (e)    coordonnees d'un 1er vecteur de travail          *
  27. *     W3        (e)    coordonnees d'un 2nd vecteur de travail          *
  28. *     CENTR1    (e)    coordonnees du 1er point repere                  *
  29. *     CENTR2    (e)    coordonnees du 2nd point repere                  *
  30. *     AXEI1     (e)    coordonnees du vecteur de l'axe de symetrie      *
  31. *     IER1      (s)    code d'erreur pour desactivation dans RTENS      *
  32. *     D.R.-M. le 17/3/94                                                *
  33. *-----------------------------------------------------------------------*
  34.  
  35. -INC PPARAM
  36. -INC CCOPTIO
  37. -INC CCHAMP
  38.  
  39. -INC SMCHAML
  40. -INC SMINTE
  41. -INC SMCOORD
  42. -INC SMELEME
  43.  
  44. -INC TMPTVAL
  45. *
  46. *     MWRK1,3,4 initialises dans RTENS1
  47. *
  48.       SEGMENT MWRK1
  49.         REAL*8 XEL(3,NBNN),XEL2(3,NBNN)
  50.       ENDSEGMENT
  51. *
  52.       SEGMENT MWRK3
  53.         REAL*8 A(NDIM,NDIM),R(NDIM,NDIM),RT(NDIM,NDIM),TRAV(NDIM,NDIM)
  54.       ENDSEGMENT
  55. *
  56.       SEGMENT MWRK4
  57.         REAL*8 XLOC(3,3),XGLOB(3,3)
  58.         REAL*8 TXR1(IDIM,IDIM),VALVEC(NVEC)
  59.       ENDSEGMENT
  60. *
  61.       DIMENSION VECWRK(3),V1(4),V2(4),W2(3),W3(3)
  62.       DIMENSION CENTR1(3),CENTR2(3),AXEI1(3),VECX(3),VECY(3)
  63.       DIMENSION UR(3),UTHETA(3),UPHI(3),UN(3),UT(3),XIGAU(3)
  64. *
  65.       IER1   = 0
  66.       BIDON  = 0.D0
  67.       ERRAXI = 0
  68.       MELEME = IELEME
  69.       NBNN   = NUM(/1)
  70.       NBELEM = NUM(/2)
  71.       MINTE  = IINTE
  72.       NBPGAU = POIGAU(/1)
  73. *
  74.       NDIM=IDIM
  75.       IF (IFOMEM.EQ.1) NDIM=IDIM+1
  76.       SEGINI MWRK3
  77. *
  78.       IF (IPCHE1.EQ.0.AND.IMOT.EQ.0) THEN
  79. *
  80. *        Repere cartesien : on veut le tenseur dans un repere defini
  81. *        par un ou deux vecteurs. On construit la matrice de rotation
  82. *        qui fait passer du repere general au repere defini par V1
  83. *        (et V2 en 3D)
  84. *
  85.          IF (IDIM.EQ.2) THEN
  86.             R(1,1)=V1(1)/V1(4)
  87.             R(2,1)=V1(2)/V1(4)
  88.             R(1,2)= -1.D0 * V1(2)/V1(4)
  89.             R(2,2)=V1(1)/V1(4)
  90.             IF (IFOMEM.EQ.1) THEN
  91.                R(3,3) = 1.D0
  92.             ENDIF
  93.          ELSE IF (IDIM.EQ.3) THEN
  94.             IF (IPTV2.EQ.0) THEN
  95.                CALL ERREUR(338)
  96.                SEGSUP MWRK3
  97.                IER1 = 1
  98.                RETURN
  99.             ENDIF
  100.             R(1,1)=V1(1)/V1(4)
  101.             R(2,1)=V1(2)/V1(4)
  102.             R(3,1)=V1(3)/V1(4)
  103.             R(1,2)=W2(1)
  104.             R(2,2)=W2(2)
  105.             R(3,2)=W2(3)
  106.             R(1,3)=W3(1)
  107.             R(2,3)=W3(2)
  108.             R(3,3)=W3(3)
  109.          ENDIF
  110.          CALL TRSPOD(R,NDIM,NDIM,RT)
  111.       ELSE
  112. *
  113. *     Le repere choisi n'est pas cartesien. On recupere les fonctions
  114. *     de forme et leurs derivees au centre de l'element pour calculer
  115. *     les axes locaux
  116. *
  117.          NLG=NUMGEO(MELE)
  118.          CALL RESHPT(1,NBNN,NLG,MELE,NPINT,IPT1,IRT1)
  119.          MINTE2=IPT1
  120.          SEGACT MINTE2
  121.          SEGINI MWRK4,MWRK1
  122.       ENDIF
  123. *
  124. *     Boucle sur les elements
  125. *
  126.       DO 1010  IB=1,NBELEM
  127. *
  128. *        Recherche des coordonnees des noeuds de l'element IB
  129. *
  130.          IF (IMOT.NE.0.OR.IPCHE1.NE.0)
  131.      $      CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XEL)
  132. *
  133.          IF (IPCHE1.NE.0) THEN
  134. *
  135. *           >>> Repere d'Orthotropie <<<
  136. *
  137. *           Calcul des axes locaux pour les materiaux orthotropes,
  138. *           anisotropes et unidirectionnels
  139. *
  140.             NBSH=MINTE2.SHPTOT(/2)
  141.             CALL RLOCAL(XEL,MINTE2.SHPTOT,NBSH,NBNN,TXR1)
  142.         if (nbsh.eq.-1) then
  143.          call erreur(525)
  144.          return
  145.         endif
  146.             IF (IERR.NE.0) THEN
  147.                SEGSUP MWRK1,MWRK3,MWRK4
  148.                IER1 = 1
  149.                RETURN
  150.             ENDIF
  151. *
  152. *           RECHERCHE DE NBGMAX
  153. *
  154.             NBGMAX=1
  155.             MPTVAL=IVAVEC
  156.             DO 1311 IV=1,NVEC
  157.                IF (IVAL(IV).NE.0) THEN
  158.                   MELVAL=IVAL(IV)
  159.                   NBGMAX=MAX(NBGMAX,VELCHE(/1))
  160.                ENDIF
  161.  1311       CONTINUE
  162.  
  163.          ENDIF
  164.  
  165.  
  166. *
  167. *        Boucle sur les points de Gauss
  168. *
  169.  
  170.          DO 1010 IGAU=1,NBPGAU
  171.  
  172.  
  173.  
  174. *
  175. *           >>> Repere d'Orthotropie <<<
  176. *
  177. *------------------------------------------------------------
  178. *     MLR  13/8/99  ON MET LE CALCUL DES AXES DANS LA BOUCLE
  179. *    SUR LES POINTS DE GAUSS ET NON PAS EN DEHORS
  180. *------------------------------------------------------------
  181.  
  182.          IF (IPCHE1.NE.0) THEN
  183.            IF(IGAU.EQ.1.OR.NBGMAX.GT.1) THEN
  184.  
  185.             MPTVAL=IVAVEC
  186.  
  187.             DO 1011 IV=1,NVEC
  188.                IF (IVAL(IV).NE.0) THEN
  189.                   MELVAL=IVAL(IV)
  190.                   IBMN=MIN(IB,VELCHE(/2))
  191.                   IGMN=MIN(IGAU,VELCHE(/1))
  192.                   VALVEC(IV)=VELCHE(IGMN,IBMN)
  193.                ELSE
  194.                   VALVEC(IV)=0.D0
  195.                ENDIF
  196.  1011       CONTINUE
  197.             CALL RGLOB(VALVEC,IDIM,TXR1,XLOC,XGLOB,IFOMEM)
  198.             IF (IERR.NE.0) THEN
  199.                SEGSUP MWRK1,MWRK3,MWRK4
  200.                 IER1 = 1
  201.                 RETURN
  202.             ENDIF
  203.             DO 1012 IC=1,IDIM
  204.             DO 1012 IL=1,IDIM
  205.                R(IL,IC)=XGLOB(IL,IC)
  206.  1012       CONTINUE
  207.             IF (IDIM.EQ.2.AND.IFOMEM.EQ.1) R(3,3)=1.D0
  208.             CALL TRSPOD(R,NDIM,NDIM,RT)
  209.  
  210.            ENDIF
  211.          ENDIF
  212.  
  213.  
  214. *------------------------------------------------------------
  215.  
  216.  
  217. *
  218. *           Sous-zones du MCHAML avant rotation
  219. *
  220.             MPTVAL=IVACOM
  221. *
  222. *           Initialisations
  223. *
  224.             IF (IMOT.NE.0) THEN
  225.                DO 1013 IC = 1,IDIM
  226.                XIGAU(IC) = 0.D0
  227.                DO 1013 IL=1,XEL(/2)
  228.                XIGAU(IC)=XIGAU(IC)+(SHPTOT(1,IL,IGAU)*XEL(IC,IL))
  229.  1013          CONTINUE
  230. *
  231.                SCAL=0.D0
  232.                DO 1014 IL=1,IDIM
  233.                   UTHETA(IL)=0.D0
  234.                   UPHI(IL)=0.D0
  235.                   UT(IL)=0.D0
  236.                   VECX(IL)=0.D0
  237.                   VECY(IL)=0.D0
  238.                   UR(IL)=XIGAU(IL)-CENTR1(IL)
  239.                   UN(IL)=UR(IL)
  240.  1014          CONTINUE
  241.                DO 1015 IL=1,IDIM
  242.                   SCAL=SCAL+UR(IL)*V1(IL)
  243.  1015          CONTINUE
  244.                DO 1016 IL=1,IDIM
  245.                   UR(IL)=UR(IL)-SCAL*V1(IL)
  246.  1016          CONTINUE
  247. *
  248.                SCAL=0.D0
  249.                DO 1017 IL=1,IDIM
  250.                   SCAL=SCAL+UR(IL)**2
  251.  1017          CONTINUE
  252.                SCAL = SQRT(SCAL)
  253.                IF (SCAL.EQ.0.D0) THEN
  254.                   CALL ERREUR(642)
  255.                   IER1 = 1
  256.                   GOTO 1010
  257.                ENDIF
  258.                IF (IDIM.EQ.3) THEN
  259.                   CALL NORMER(UR)
  260.                   CALL NORMER(UN)
  261.                   CALL PROVEC(V1,UR,UTHETA)
  262.                ELSE
  263. *
  264. *                 Dimension 2 : 'POLA'
  265. *
  266.                   UR(1)=UR(1)/SCAL
  267.                   UR(2)=UR(2)/SCAL
  268.                   R(1,1)=UR(1)
  269.                   R(1,2)=-UR(2)
  270.                   R(2,2)=UR(1)
  271.                   R(2,1)=UR(2)
  272.                   IF (IFOMEM.EQ.1) THEN
  273.                      R(3,3)=1D0
  274.                   ENDIF
  275.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  276.                ENDIF
  277. *
  278. *              Debut des calculs de R pour les autres reperes
  279. *
  280. *              -- Cas CYLINDRIQUE --
  281. *
  282.                IF (IMOT.EQ.2) THEN
  283.                   DO 1019 IL=1,IDIM
  284.                      R(IL,1)=UR(IL)
  285.                      R(IL,2)=UTHETA(IL)
  286.                      R(IL,3)=V1(IL)
  287.  1019             CONTINUE
  288.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  289.                ELSE
  290.  
  291.                ENDIF
  292. *
  293. *              -- Cas SPHERIQUE --
  294. *
  295.                IF (IMOT.EQ.3) THEN
  296.                   UR(1)=UN(1)
  297.                   UR(2)=UN(2)
  298.                   UR(3)=UN(3)
  299.                   UPHI(1)=UTHETA(1)
  300.                   UPHI(2)=UTHETA(2)
  301.                   UPHI(3)=UTHETA(3)
  302.                   CALL PROVEC (UPHI,UR,UTHETA)
  303.                   DO 1021 IL=1,IDIM
  304.                      R(IL,1)=UR(IL)
  305.                      R(IL,2)=UTHETA(IL)
  306.                      R(IL,3)=UPHI(IL)
  307.  1021             CONTINUE
  308.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  309.                ENDIF
  310. *
  311. *              -- Cas TORIQUE CIRCULAIRE --
  312. *
  313.                IF (IMOT.EQ.4) THEN
  314.                   VECWRK(1)=CENTR2(1)-CENTR1(1)
  315.                   VECWRK(2)=CENTR2(2)-CENTR1(2)
  316.                   VECWRK(3)=CENTR2(3)-CENTR1(3)
  317.                   CALL NORME(VECWRK,SCAL)
  318.                   UN(1)=UN(1)-SCAL*UR(1)
  319.                   UN(2)=UN(2)-SCAL*UR(2)
  320.                   UN(3)=UN(3)-SCAL*UR(3)
  321.                   CALL NORMER(UN)
  322.                   CALL PROVEC(UN,UTHETA,UT)
  323.                   DO 1022 IL=1,IDIM
  324.                      R(IL,1)=UTHETA(IL)
  325.                      R(IL,2)=UT(IL)
  326.                      R(IL,3)=UN(IL)
  327.  1022             CONTINUE
  328.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  329.                ENDIF
  330. *
  331. *              -- Cas TORIQUE CARTESIEN --
  332. *
  333.                IF (IMOT.EQ.5) THEN
  334.                   DO 1023 IL=1,IDIM
  335.                      R(IL,1)=UR(IL)
  336.                      R(IL,2)=UTHETA(IL)
  337.                      R(IL,3)=V1(IL)
  338.  1023             CONTINUE
  339.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  340.                ENDIF
  341.             ENDIF
  342. *
  343. *           Tenseur avant changement de repere
  344. *
  345.             MELVAL=IVAL(1)
  346.             IGMN = MIN(IGAU,VELCHE(/1))
  347.             IBMN = MIN(IB,  VELCHE(/2))
  348.             A(1,1) =  VELCHE(IGMN,IBMN)
  349. *
  350.             MELVAL=IVAL(2)
  351.             IGMN = MIN(IGAU,VELCHE(/1))
  352.             IBMN = MIN(IB,  VELCHE(/2))
  353.             A(2,2) = VELCHE(IGMN,IBMN)
  354. *
  355.             MELVAL=IVAL(4)
  356.             IGMN = MIN(IGAU,VELCHE(/1))
  357.             IBMN = MIN(IB,  VELCHE(/2))
  358.             A(1,2) = VELCHE(IGMN,IBMN)
  359. *
  360.             IF (IDEFO.EQ.1) A(1,2)=A(1,2)/2.D0
  361.             A(2,1)=A(1,2)
  362. *
  363.             IF (IFOMEM.LT.1) GOTO 6610
  364. *
  365.             MELVAL=IVAL(3)
  366.             IGMN = MIN(IGAU,VELCHE(/1))
  367.             IBMN = MIN(IB,  VELCHE(/2))
  368.             A(3,3) = VELCHE(IGMN,IBMN)
  369. *
  370.             MELVAL=IVAL(5)
  371.             IGMN = MIN(IGAU,VELCHE(/1))
  372.             IBMN = MIN(IB,  VELCHE(/2))
  373.             A(3,1) = VELCHE(IGMN,IBMN)
  374. *
  375.             MELVAL=IVAL(6)
  376.             IGMN = MIN(IGAU,VELCHE(/1))
  377.             IBMN = MIN(IB,  VELCHE(/2))
  378.             A(3,2) = VELCHE(IGMN,IBMN)
  379. *
  380.             IF (IDEFO.EQ.1) A(3,1)=A(3,1)/2.D0
  381.             IF (IDEFO.EQ.1) A(3,2)=A(3,2)/2.D0
  382.             A(1,3)=A(3,1)
  383.             A(2,3)=A(3,2)
  384. *
  385.             MELVAL=IVAL(3)
  386.             IGMN = MIN(IGAU,VELCHE(/1))
  387.             IBMN = MIN(IB,  VELCHE(/2))
  388.             A(3,3) = VELCHE(IGMN,IBMN)
  389. *
  390.  6610       CONTINUE
  391. *
  392.             MELVAL=IVAL(3)
  393.             IGMN = MIN(IGAU,VELCHE(/1))
  394.             IBMN = MIN(IB,  VELCHE(/2))
  395.             AUX = VELCHE(IGMN,IBMN)
  396. *                                           t
  397. *           >>>  Rotation du tenseur : A = R  A  R  <<<
  398. *
  399.             CALL MULMAT(TRAV,A,R,NDIM,NDIM,NDIM)
  400.             CALL MULMAT(A,RT,TRAV,NDIM,NDIM,NDIM)
  401. *
  402. *           Tenseur apres changement de repere
  403. *           Sous-zones du MCHAML resultat
  404. *
  405.             MPTVAL=IVARES
  406. *
  407.             MELVAL=IVAL(1)
  408.             VELCHE(IGAU,IB) = A(1,1)
  409. *
  410.             MELVAL=IVAL(2)
  411.             VELCHE(IGAU,IB) = A(2,2)
  412. *
  413.             IF (IDEFO.EQ.1) A(1,2)=A(1,2)*2.D0
  414. *
  415.             MELVAL=IVAL(4)
  416.             VELCHE(IGAU,IB) = A(1,2)
  417. *
  418.             IF (IFOMEM.LT.1) THEN
  419. *
  420.                MELVAL=IVAL(3)
  421.                VELCHE(IGAU,IB)= AUX
  422. *
  423.             ELSE
  424. *
  425.                MELVAL=IVAL(3)
  426.                VELCHE(IGAU,IB)=A(3,3)
  427. *
  428.                IF (IDEFO.EQ.1) A(3,1)=A(3,1)*2.D0
  429.                IF (IDEFO.EQ.1) A(3,2)=A(3,2)*2.D0
  430. *
  431.                MELVAL=IVAL(5)
  432.                VELCHE(IGAU,IB)= A(3,1)
  433. *
  434.                MELVAL=IVAL(6)
  435.                VELCHE(IGAU,IB)=A(3,2)
  436. *
  437.             ENDIF
  438. *
  439.  1010 CONTINUE
  440.       SEGSUP MWRK3
  441.       IF (IPCHE1.NE.0.OR.IMOT.NE.0) THEN
  442.          SEGSUP MWRK1,MWRK4
  443.          SEGDES MINTE2
  444.       ENDIF
  445.  
  446.       RETURN
  447.       END
  448.  
  449.  
  450.  

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