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

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