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    CB215821  26/08/24    21:18:19     12622          
  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 6611  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 6612 IC=1,IDIM
  204.             DO 1012 IL=1,IDIM
  205.                R(IL,IC)=XGLOB(IL,IC)
  206.  1012       CONTINUE
  207.  6612       CONTINUE
  208.             IF (IDIM.EQ.2.AND.IFOMEM.EQ.1) R(3,3)=1.D0
  209.             CALL TRSPOD(R,NDIM,NDIM,RT)
  210.  
  211.            ENDIF
  212.          ENDIF
  213.  
  214.  
  215. *------------------------------------------------------------
  216.  
  217.  
  218. *
  219. *           Sous-zones du MCHAML avant rotation
  220. *
  221.             MPTVAL=IVACOM
  222. *
  223. *           Initialisations
  224. *
  225.             IF (IMOT.NE.0) THEN
  226.                DO 6613 IC = 1,IDIM
  227.                XIGAU(IC) = 0.D0
  228.                DO 1013 IL=1,XEL(/2)
  229.                XIGAU(IC)=XIGAU(IC)+(SHPTOT(1,IL,IGAU)*XEL(IC,IL))
  230.  1013          CONTINUE
  231.  6613          CONTINUE
  232. *
  233.                SCAL=0.D0
  234.                DO 1014 IL=1,IDIM
  235.                   UTHETA(IL)=0.D0
  236.                   UPHI(IL)=0.D0
  237.                   UT(IL)=0.D0
  238.                   VECX(IL)=0.D0
  239.                   VECY(IL)=0.D0
  240.                   UR(IL)=XIGAU(IL)-CENTR1(IL)
  241.                   UN(IL)=UR(IL)
  242.  1014          CONTINUE
  243.                DO 1015 IL=1,IDIM
  244.                   SCAL=SCAL+UR(IL)*V1(IL)
  245.  1015          CONTINUE
  246.                DO 1016 IL=1,IDIM
  247.                   UR(IL)=UR(IL)-SCAL*V1(IL)
  248.  1016          CONTINUE
  249. *
  250.                SCAL=0.D0
  251.                DO 1017 IL=1,IDIM
  252.                   SCAL=SCAL+UR(IL)**2
  253.  1017          CONTINUE
  254.                SCAL = SQRT(SCAL)
  255.                IF (SCAL.EQ.0.D0) THEN
  256.                   CALL ERREUR(642)
  257.                   IER1 = 1
  258.                   GOTO 1010
  259.                ENDIF
  260.                IF (IDIM.EQ.3) THEN
  261.                   CALL NORMER(UR)
  262.                   CALL NORMER(UN)
  263.                   CALL PROVEC(V1,UR,UTHETA)
  264.                ELSE
  265. *
  266. *                 Dimension 2 : 'POLA'
  267. *
  268.                   UR(1)=UR(1)/SCAL
  269.                   UR(2)=UR(2)/SCAL
  270.                   R(1,1)=UR(1)
  271.                   R(1,2)=-UR(2)
  272.                   R(2,2)=UR(1)
  273.                   R(2,1)=UR(2)
  274.                   IF (IFOMEM.EQ.1) THEN
  275.                      R(3,3)=1D0
  276.                   ENDIF
  277.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  278.                ENDIF
  279. *
  280. *              Debut des calculs de R pour les autres reperes
  281. *
  282. *              -- Cas CYLINDRIQUE --
  283. *
  284.                IF (IMOT.EQ.2) THEN
  285.                   DO 1019 IL=1,IDIM
  286.                      R(IL,1)=UR(IL)
  287.                      R(IL,2)=UTHETA(IL)
  288.                      R(IL,3)=V1(IL)
  289.  1019             CONTINUE
  290.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  291.                ELSE
  292.  
  293.                ENDIF
  294. *
  295. *              -- Cas SPHERIQUE --
  296. *
  297.                IF (IMOT.EQ.3) THEN
  298.                   UR(1)=UN(1)
  299.                   UR(2)=UN(2)
  300.                   UR(3)=UN(3)
  301.                   UPHI(1)=UTHETA(1)
  302.                   UPHI(2)=UTHETA(2)
  303.                   UPHI(3)=UTHETA(3)
  304.                   CALL PROVEC (UPHI,UR,UTHETA)
  305.                   DO 1021 IL=1,IDIM
  306.                      R(IL,1)=UR(IL)
  307.                      R(IL,2)=UTHETA(IL)
  308.                      R(IL,3)=UPHI(IL)
  309.  1021             CONTINUE
  310.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  311.                ENDIF
  312. *
  313. *              -- Cas TORIQUE CIRCULAIRE --
  314. *
  315.                IF (IMOT.EQ.4) THEN
  316.                   VECWRK(1)=CENTR2(1)-CENTR1(1)
  317.                   VECWRK(2)=CENTR2(2)-CENTR1(2)
  318.                   VECWRK(3)=CENTR2(3)-CENTR1(3)
  319.                   CALL NORME(VECWRK,SCAL)
  320.                   UN(1)=UN(1)-SCAL*UR(1)
  321.                   UN(2)=UN(2)-SCAL*UR(2)
  322.                   UN(3)=UN(3)-SCAL*UR(3)
  323.                   CALL NORMER(UN)
  324.                   CALL PROVEC(UN,UTHETA,UT)
  325.                   DO 1022 IL=1,IDIM
  326.                      R(IL,1)=UTHETA(IL)
  327.                      R(IL,2)=UT(IL)
  328.                      R(IL,3)=UN(IL)
  329.  1022             CONTINUE
  330.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  331.                ENDIF
  332. *
  333. *              -- Cas TORIQUE CARTESIEN --
  334. *
  335.                IF (IMOT.EQ.5) THEN
  336.                   DO 1023 IL=1,IDIM
  337.                      R(IL,1)=UR(IL)
  338.                      R(IL,2)=UTHETA(IL)
  339.                      R(IL,3)=V1(IL)
  340.  1023             CONTINUE
  341.                   CALL TRSPOD (R,NDIM,NDIM,RT)
  342.                ENDIF
  343.             ENDIF
  344. *
  345. *           Tenseur avant changement de repere
  346. *
  347.             MELVAL=IVAL(1)
  348.             IGMN = MIN(IGAU,VELCHE(/1))
  349.             IBMN = MIN(IB,  VELCHE(/2))
  350.             A(1,1) =  VELCHE(IGMN,IBMN)
  351. *
  352.             MELVAL=IVAL(2)
  353.             IGMN = MIN(IGAU,VELCHE(/1))
  354.             IBMN = MIN(IB,  VELCHE(/2))
  355.             A(2,2) = VELCHE(IGMN,IBMN)
  356. *
  357.             MELVAL=IVAL(4)
  358.             IGMN = MIN(IGAU,VELCHE(/1))
  359.             IBMN = MIN(IB,  VELCHE(/2))
  360.             A(1,2) = VELCHE(IGMN,IBMN)
  361. *
  362.             IF (IDEFO.EQ.1) A(1,2)=A(1,2)/2.D0
  363.             A(2,1)=A(1,2)
  364. *
  365.             IF (IFOMEM.LT.1) GOTO 6610
  366. *
  367.             MELVAL=IVAL(3)
  368.             IGMN = MIN(IGAU,VELCHE(/1))
  369.             IBMN = MIN(IB,  VELCHE(/2))
  370.             A(3,3) = VELCHE(IGMN,IBMN)
  371. *
  372.             MELVAL=IVAL(5)
  373.             IGMN = MIN(IGAU,VELCHE(/1))
  374.             IBMN = MIN(IB,  VELCHE(/2))
  375.             A(3,1) = VELCHE(IGMN,IBMN)
  376. *
  377.             MELVAL=IVAL(6)
  378.             IGMN = MIN(IGAU,VELCHE(/1))
  379.             IBMN = MIN(IB,  VELCHE(/2))
  380.             A(3,2) = VELCHE(IGMN,IBMN)
  381. *
  382.             IF (IDEFO.EQ.1) A(3,1)=A(3,1)/2.D0
  383.             IF (IDEFO.EQ.1) A(3,2)=A(3,2)/2.D0
  384.             A(1,3)=A(3,1)
  385.             A(2,3)=A(3,2)
  386. *
  387.             MELVAL=IVAL(3)
  388.             IGMN = MIN(IGAU,VELCHE(/1))
  389.             IBMN = MIN(IB,  VELCHE(/2))
  390.             A(3,3) = VELCHE(IGMN,IBMN)
  391. *
  392.  6610       CONTINUE
  393. *
  394.             MELVAL=IVAL(3)
  395.             IGMN = MIN(IGAU,VELCHE(/1))
  396.             IBMN = MIN(IB,  VELCHE(/2))
  397.             AUX = VELCHE(IGMN,IBMN)
  398. *                                           t
  399. *           >>>  Rotation du tenseur : A = R  A  R  <<<
  400. *
  401.             CALL MULMAT(TRAV,A,R,NDIM,NDIM,NDIM)
  402.             CALL MULMAT(A,RT,TRAV,NDIM,NDIM,NDIM)
  403. *
  404. *           Tenseur apres changement de repere
  405. *           Sous-zones du MCHAML resultat
  406. *
  407.             MPTVAL=IVARES
  408. *
  409.             MELVAL=IVAL(1)
  410.             VELCHE(IGAU,IB) = A(1,1)
  411. *
  412.             MELVAL=IVAL(2)
  413.             VELCHE(IGAU,IB) = A(2,2)
  414. *
  415.             IF (IDEFO.EQ.1) A(1,2)=A(1,2)*2.D0
  416. *
  417.             MELVAL=IVAL(4)
  418.             VELCHE(IGAU,IB) = A(1,2)
  419. *
  420.             IF (IFOMEM.LT.1) THEN
  421. *
  422.                MELVAL=IVAL(3)
  423.                VELCHE(IGAU,IB)= AUX
  424. *
  425.             ELSE
  426. *
  427.                MELVAL=IVAL(3)
  428.                VELCHE(IGAU,IB)=A(3,3)
  429. *
  430.                IF (IDEFO.EQ.1) A(3,1)=A(3,1)*2.D0
  431.                IF (IDEFO.EQ.1) A(3,2)=A(3,2)*2.D0
  432. *
  433.                MELVAL=IVAL(5)
  434.                VELCHE(IGAU,IB)= A(3,1)
  435. *
  436.                MELVAL=IVAL(6)
  437.                VELCHE(IGAU,IB)=A(3,2)
  438. *
  439.             ENDIF
  440. *
  441.  1010 CONTINUE
  442.  6611 CONTINUE
  443.       SEGSUP MWRK3
  444.       IF (IPCHE1.NE.0.OR.IMOT.NE.0) THEN
  445.          SEGSUP MWRK1,MWRK4
  446.          SEGDES MINTE2
  447.       ENDIF
  448.  
  449.       RETURN
  450.       END
  451.  
  452.  
  453.  
  454.  

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