Télécharger la plaquette Cast3M


fr en

Télécharger jacono.eso

Retour à la liste

Numérotation des lignes : 

jacono

  1. C JACONO    SOURCE    GOUNAND   26/07/30    21:15:06     12610          
  2.  
  3. C=======================================================================
  4. C    ENTREES :
  5. C    ---------
  6. C      IPMODL= pointeur sur un MMODEL
  7. C      INORM = 1 si les vecteurs doivent etre normes 0 sinon
  8. C
  9. C    SORTIES :
  10. C    --------
  11. C
  12. C      IPCHE = CHAMELEM contenant les JACOBIENS
  13. C          ( = normale aux faces des elements dans le cas des coques)
  14. C          ( = tangente a la fibre neutre dans le cas des poutres)
  15. C      IRET  = 1 si succes 0 sinon
  16. C
  17. C    Passage au nouveau Chamelem PAR S.RAMAHANDRY le 11/09/90
  18. C
  19. C    2013-01-02 (BP) : ajout zones cohesives (ZCO2,3 et 4 => coque mince)
  20. C                      + calcul de la tangente pour les poutres
  21. C
  22. C
  23. C=====================================================================
  24.       SUBROUTINE JACONO(IPMODL,INORM,IPCHE,IRET)
  25.  
  26.       IMPLICIT INTEGER(I-N)
  27.       IMPLICIT REAL*8(A-H,O-Z)
  28.  
  29. -INC PPARAM
  30. -INC CCOPTIO
  31. -INC CCHAMP
  32. -INC CCREEL
  33.  
  34. -INC SMCHAML
  35. -INC SMMODEL
  36. -INC SMELEME
  37. -INC SMCOORD
  38. -INC SMINTE
  39.  
  40. -INC TMPTVAL
  41.  
  42.       SEGMENT TRA
  43.          REAL*8 XEL(3,NBNN) ,SHP(6,NBNN) ,XE(3,NBNN)
  44.       ENDSEGMENT
  45. C
  46.       SEGMENT TR1
  47.          REAL*8 TH(NBN1) ,TXR(3,3,NBN1) ,XJ(3,3)
  48.       ENDSEGMENT
  49. C
  50.       PARAMETER(UN=1.D0,XZER=0.D0)
  51.       DIMENSION BPSS(3,3)
  52.  
  53.       DIMENSION XU(3), XV(3), XW(3)
  54.  
  55.       IDIMP1 = IDIM+1
  56.       NHRM=NIFOUR
  57.       IRET=1
  58.       XTOL=XZPREC*10.D0
  59. C
  60. C     ACTIVATION DU MODELE
  61. C
  62.       MMODEL= IPMODL
  63.       SEGACT MMODEL
  64.       NSOUS=KMODEL(/1)
  65. C
  66. C     CREATION DU MCHELM
  67. C
  68.       N1=NSOUS
  69.       N3=6
  70.       IF (INORM .EQ. 1) THEN
  71.         L1=8
  72.       ELSE
  73.         L1=16
  74.       ENDIF
  75.       SEGINI MCHELM
  76.       IF (INORM .EQ. 1) THEN
  77.          TITCHE='NORMALES'
  78.       ELSE
  79.          TITCHE='VECTEURS SURFACE'
  80.       ENDIF
  81.       IFOCHE=IFOUR
  82.       IPCHE=MCHELM
  83. C____________________________________________________________________
  84. C
  85. C     DEBUT DE LA BOUCLE SUR LES DIFFERENTES ZONES
  86. C____________________________________________________________________
  87. C
  88.       DO 500 ISOUS=1,NSOUS
  89. C
  90. C     ON RECUPERE L INFORMATION GENERALE
  91. C
  92.          IMODEL=KMODEL(ISOUS)
  93.          SEGACT IMODEL
  94.          IPMAIL=IMAMOD
  95.          IMACHE(ISOUS)=IPMAIL
  96.          CONCHE(ISOUS)=CONMOD
  97. C
  98. C     TRAITEMENT DU MODELE
  99. C
  100.          MELE=NEFMOD
  101.          MELEME=IMAMOD
  102.          NFOR=FORMOD(/2)
  103.          NMAT=MATMOD(/2)
  104. C____________________________________________________________________
  105. C
  106. C     INFORMATION SUR L'ELEMENT FINI
  107. C____________________________________________________________________
  108. C
  109.          MELE =INFELE(1)
  110.          MFR  =INFELE(13)
  111.          MINTE=INFMOD(7)
  112.          MINTE1=INFMOD(3)
  113. C
  114.          INFCHE(ISOUS,1)=0
  115.          INFCHE(ISOUS,2)=0
  116.          INFCHE(ISOUS,3)=NHRM
  117.          INFCHE(ISOUS,4)=MINTE
  118.          INFCHE(ISOUS,5)=0
  119.          INFCHE(ISOUS,6)=5
  120. C
  121. C     INITIALISATION DE MINTE
  122. C
  123.          SEGACT MINTE
  124.          NBPGAU=POIGAU(/1)
  125. C
  126. C     ACTIVATION DU MELEME
  127. C
  128.          SEGACT MELEME
  129.          NBNN  =NUM(/1)
  130.          NBELEM=NUM(/2)
  131. C
  132. C     RECHERCHE DE LA TAILLE DES MELVAL A ALLOUER
  133. C
  134. * COQUE
  135.          IF (MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9 .OR. MFR.EQ.77) THEN
  136.             N1PTEL=NBPGAU
  137.             N1EL=NBELEM
  138. * POUTRE TUYAU UNIAXIALE
  139.          ELSEIF(MFR.EQ.7 .OR. MFR.EQ.13 .OR. MFR.EQ.27) THEN
  140.             N1PTEL=NBPGAU
  141.             N1EL=NBELEM
  142.          ELSE
  143.             N1PTEL = 0
  144.             N1EL = 0
  145.          ENDIF
  146.          N2PTEL=0
  147.          N2EL  =0
  148. C
  149. C     CREATION DU MCHAML DE LA SOUS ZONE
  150. C
  151.          NJAC=IDIM
  152.          N2 = NJAC
  153.          SEGINI MCHAML
  154.          ICHAML(ISOUS)=MCHAML
  155.          NSR   =1
  156.          NCOSOR=NJAC
  157.          SEGINI MPTVAL
  158.          IVAJAC=MPTVAL
  159. C
  160. C     2 OU 3 COMPOSANTES
  161. C
  162.          ICOMP=1
  163.          IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN
  164.             NOMCHE(ICOMP)='VR  '
  165.          ELSE
  166.             NOMCHE(ICOMP)='VX  '
  167.          ENDIF
  168.          TYPCHE(ICOMP)='REAL*8'
  169.          SEGINI MELVA1
  170.          IELVAL(ICOMP)=MELVA1
  171.          IVAL(ICOMP)=MELVA1
  172. C
  173.          ICOMP=2
  174.          IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN
  175.             NOMCHE(ICOMP)='VZ  '
  176.          ELSE
  177.             NOMCHE(ICOMP)='VY  '
  178.          ENDIF
  179.          TYPCHE(ICOMP)='REAL*8'
  180.          SEGINI MELVA2
  181.          IELVAL(ICOMP)=MELVA2
  182.          IVAL(ICOMP)=MELVA2
  183. C
  184.          MELVA3 = 0
  185.          IF (IDIM .EQ. 3) THEN
  186.             ICOMP=3
  187.             NOMCHE(ICOMP)='VZ  '
  188.             TYPCHE(ICOMP)='REAL*8'
  189.             SEGINI MELVA3
  190.             IELVAL(ICOMP)=MELVA3
  191.             IVAL(ICOMP)=MELVA3
  192.          ENDIF
  193. C
  194.          SEGINI TRA
  195. C
  196. C ================ FORMULATION MASSIVE =======================
  197. C
  198.          IF(MFR.EQ.1.OR.MFR.EQ.33) THEN
  199.             GOTO 520
  200. C
  201. C ================ FORMULATION COQUE  MINCE =====================
  202. C
  203. C
  204.          ELSE IF(MFR.EQ.3.OR.MFR.EQ.9 .OR. MFR.EQ.77) THEN
  205.           IDI2=IDIM-1
  206.           DO 3000 IB=1,NBELEM
  207. *--------------Calcul de la normale a l'élément
  208.             INU1=MELEME.NUM(1,IB)
  209.             INU2=MELEME.NUM(2,IB)
  210.             IREF1 = IDIMP1*(INU1 - 1)
  211.             IREF2 = IDIMP1*(INU2 - 1)
  212.             XNORU = 0.D0
  213.             DO IC = 1, IDIM
  214.                r_z = XCOOR(IREF2+IC)-XCOOR(IREF1+IC)
  215.                XU(IC) = r_z
  216.                XNORU = XNORU + (r_z * r_z)
  217.             ENDDO
  218.             XNORU = SQRT(XNORU)
  219.             DO IC = 1, IDIM
  220.                XU(IC) = XU(IC)/XNORU
  221.             ENDDO
  222.             IF (IDIM .EQ. 2) THEN
  223.                XW(1) = -XU(2)
  224.                XW(2) = XU(1)
  225.             ELSE
  226.                IN = 3
  227.  33            CONTINUE
  228.                INUIN=MELEME.NUM(IN,IB)
  229.                IREF3 = IDIMP1*(INUIN - 1)
  230.                XNORV = 0.D0
  231.                DO IC = 1, IDIM
  232.                   r_z = XCOOR(IREF3+IC)-XCOOR(IREF1+IC)
  233.                   XV(IC) = r_z
  234.                   XNORV = XNORV + (r_z * r_z)
  235.                ENDDO
  236.                XNORV = SQRT(XNORV)
  237.                DO IC = 1, IDIM
  238.                  XV(IC) = XV(IC)/XNORV
  239.                ENDDO
  240.                XW(1) = XU(2)*XV(3) - XU(3)*XV(2)
  241.                XW(2) = XU(3)*XV(1) - XU(1)*XV(3)
  242.                XW(3) = XU(1)*XV(2) - XU(2)*XV(1)
  243.                XNORW = 0.
  244.                DO IC = 1, IDIM
  245.                  XNORW = XNORW + XW(IC)*XW(IC)
  246.                ENDDO
  247.                IF (XNORW .LE. XTOL) THEN
  248.                   write(ioimp,*) '  IN=',IN
  249.                   write(ioimp,*) '  W1,W2,W3=',XW(1),XW(2),XW(3)
  250.                   write(ioimp,*) '  XNORW2=',XNORW
  251.                   IN = IN + 1
  252.                   if(IN.le.NBNN) GOTO 33
  253.                ENDIF
  254.                XNORW = SQRT(XNORW)
  255.                IF(XNORW .LE.XTOL) then
  256.                   write(ioimp,*) 'INU1,INU2,INUIN=',INU1,INU2,INUIN
  257.                   write(ioimp,*) 'U1,U2,U3=',XU(1),XU(2),XU(3)
  258.                   write(ioimp,*) 'V1,V2,V3=',XV(1),XV(2),XV(3)
  259.                   write(ioimp,*) 'W1,W2,W3=',XW(1),XW(2),XW(3)
  260.                   write(ioimp,*) 'XNORW=',XNORW
  261.                   call erreur(345)
  262.                   return
  263.                ENDIF
  264.                DO IC = 1, IDIM
  265.                   XW(IC) = XW(IC)/XNORW
  266.                ENDDO
  267.             ENDIF
  268. *--------------Fin du calcul de la normale a l'élément
  269.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  270. C
  271.             IF(IDIM.EQ.2)THEN
  272.                CALL VPAST2(XE,BPSS)
  273.             ELSE IF(IDIM.EQ.3) THEN
  274.                CALL VPAST(XE,BPSS)
  275.             ENDIF
  276.             CALL VCORL1(XE,XEL,BPSS,NBNN)
  277.             MPTVAL=IVAJAC
  278.             DO 3002 IC=1,NBPGAU
  279.                IF (INORM .EQ. 1) THEN
  280.                   DJAC = 1.D0
  281.                ELSE
  282.                   DO IE=1,NBNN
  283.                      DO ID=1,6
  284.                         SHP(ID,IE)=SHPTOT(ID,IE,IC)
  285.                      ENDDO
  286.                   ENDDO
  287.                   CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC)
  288.                ENDIF
  289.                MELVAL = IVAL(1)
  290.                IBMN=MIN(IB, MELVAL.VELCHE(/2))
  291.                IGMN=MIN(IC,MELVAL.VELCHE(/1))
  292.                MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(1)
  293.                MELVAL = IVAL(2)
  294.                IBMN=MIN(IB, MELVAL.VELCHE(/2))
  295.                IGMN=MIN(IC,MELVAL.VELCHE(/1))
  296.                MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(2)
  297.                IF (IDIM .EQ. 3) THEN
  298.                   MELVAL = IVAL(3)
  299.                   IBMN=MIN(IB, MELVAL.VELCHE(/2))
  300.                   IGMN=MIN(IC,MELVAL.VELCHE(/1))
  301.                   MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(3)
  302.                ENDIF
  303.  3002       CONTINUE
  304.  3000    CONTINUE
  305.          GOTO 520
  306. C
  307. C ================ FORMULATION POUTRE, TUYAU et UNIAXIALE  ====================
  308. C
  309.       ELSE IF(MFR.EQ.7.OR.MFR.EQ.13.OR.MFR.EQ.27) THEN
  310.          IDI2=IDIM-1
  311.          ILAST=MELEME.NUM(/1)
  312.          DO 4000 IB=1,NBELEM
  313. *-----------Calcul de la tangente a l'élément
  314.             IREF1 = IDIMP1*(MELEME.NUM(1, IB) - 1)
  315.             IREF2 = IDIMP1*(MELEME.NUM(ILAST, IB) - 1)
  316.             XNORU = 0.D0
  317.             DO IC = 1, IDIM
  318.                r_z = XCOOR(IREF2+IC)-XCOOR(IREF1+IC)
  319.                XU(IC) = r_z
  320.                XNORU = XNORU + (r_z * r_z)
  321.             ENDDO
  322.             XNORU = SQRT(XNORU)
  323.             DO IC = 1, IDIM
  324.                XU(IC) = XU(IC)/XNORU
  325.             ENDDO
  326. *-----------Fin du calcul de la tangente a l'élément
  327. *-----------Calcul de la tangente a l'élément en chaque point de Gauss
  328. c  BP : On suppose le jacobien constant dans l'element (idem POUJAC.eso)
  329. C  => on sort le calcul du jacobien de la boucle sur les points de Gauss.
  330. C  Cela implique que la POUTre de Bernoulli n'est pas isoparamétrique...
  331.             IF (INORM .EQ. 1) THEN
  332.                DJAC = 1.D0
  333.             ELSE
  334.                DJAC = 1.D0/DBLE(NBPGAU)
  335.             ENDIF
  336.             MPTVAL=IVAJAC
  337.             DO 4002 IC=1,NBPGAU
  338.                MELVAL = IVAL(1)
  339.                IBMN=MIN(IB, MELVAL.VELCHE(/2))
  340.                IGMN=MIN(IC,MELVAL.VELCHE(/1))
  341.                MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(1)
  342.                MELVAL = IVAL(2)
  343.                IBMN=MIN(IB, MELVAL.VELCHE(/2))
  344.                IGMN=MIN(IC,MELVAL.VELCHE(/1))
  345.                MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(2)
  346.                IF (IDIM .EQ. 3) THEN
  347.                   MELVAL = IVAL(3)
  348.                   IBMN=MIN(IB, MELVAL.VELCHE(/2))
  349.                   IGMN=MIN(IC,MELVAL.VELCHE(/1))
  350.                   MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(3)
  351.                ENDIF
  352.  4002       CONTINUE
  353.  4000    CONTINUE
  354.  
  355.          GOTO 520
  356. C
  357. C ================ FORMULATION COQUE  EPAISSE ====================
  358. C
  359.          ELSE IF(MFR.EQ.5) THEN
  360.             SEGACT MINTE1
  361.             NBPGA1=MINTE1.POIGAU(/1)
  362.             NBN1  =MINTE1.SHPTOT(/2)
  363.             SEGINI TR1
  364. C
  365. C     UNE PETITE HORREUR ON CONSIDERE LES EPAISSEURS CONSTANTES
  366. C
  367.             DO 5010 IC=1,NBNN
  368.                TH(IC)=UN
  369.  5010       CONTINUE
  370.             DO 5000 IB=1,NBELEM
  371. *--------------Calcul de la normale a l'élément
  372.                IREF1 = IDIMP1*(MELEME.NUM(1, IB) - 1)
  373. *               IREF2 = IDIMP1*(MELEME.NUM(2, IB) - 1)
  374. * bp : les EF de coque epaisse etant quadratiques (coq6 et coq8), on
  375. * prend les noeuds "coins" pour eviter pb avec les noeuds 1,2,3 si courbures
  376.                IREF2 = IDIMP1*(MELEME.NUM(3, IB) - 1)
  377.                XNORU = 0.
  378.                DO IC = 1, IDIM
  379.                   r_z = XCOOR(IREF2+IC)-XCOOR(IREF1+IC)
  380.                   XU(IC) = r_z
  381.                   XNORU = XNORU + (r_z * r_z)
  382.                ENDDO
  383.                XNORU = SQRT(XNORU)
  384.                DO IC = 1, IDIM
  385.                   XU(IC) = XU(IC)/XNORU
  386.                ENDDO
  387.                IF (IDIM .EQ. 2) THEN
  388.                   XW(1) = -XU(2)
  389.                   XW(2) = XU(1)
  390.                ELSE
  391. *                  IN = 3
  392.                   IN = 5
  393.    13             IREF3 = IDIMP1*(MELEME.NUM(IN, IB) - 1)
  394.                   XNORV = 0.D0
  395.                   DO IC = 1, IDIM
  396.                     r_z = XCOOR(IREF3 + IC) - XCOOR(IREF1 + IC)
  397.                     XV(IC) = r_z
  398.                     XNORV = XNORV + (r_z * r_z)
  399.                   ENDDO
  400.                   XNORV = SQRT(XNORV)
  401.                   DO IC = 1, IDIM
  402.                     XV(IC) = XV(IC)/XNORV
  403.                   ENDDO
  404.                   XW(1) = XU(2)*XV(3) - XU(3)*XV(2)
  405.                   XW(2) = XU(3)*XV(1) - XU(1)*XV(3)
  406.                   XW(3) = XU(1)*XV(2) - XU(2)*XV(1)
  407.                   XNORW = 0.
  408.                   DO IC = 1, IDIM
  409.                      XNORW = XNORW + XW(IC)*XW(IC)
  410.                   ENDDO
  411.                   XNORW = SQRT(XNORW)
  412.                   IF (XNORW .LT. 1.E-4) THEN
  413.                    if(IN.LT.NBNN) then
  414.                      IN = IN + 1
  415.                      GOTO 13
  416.                    else
  417.                     write(6,*) 'Difficultes pour etablir la normale de'
  418.                     write(6,*) 'l element',IB
  419.                     write(6,*) 'Verifiez votre maillage'
  420.                     GOTO 9990
  421.                    endif
  422.                   ENDIF
  423.                   DO IC = 1, IDIM
  424.                      XW(IC) = XW(IC)/XNORW
  425.                   ENDDO
  426.                ENDIF
  427. *--------------Fin du calcul de la normale a l'élément
  428.                IF (INORM .EQ. 0) THEN
  429.                   CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  430. C
  431.                   CALL CQ8LOC(XE,NBN1,MINTE1.SHPTOT,TXR,IRR)
  432.                ENDIF
  433. C
  434.                MPTVAL=IVAJAC
  435.                DO 5002 IC=1,NBPGAU
  436.                   IF (INORM .EQ. 1) THEN
  437.                      DJAC = 1.D0
  438.                   ELSE
  439.                      E=DZEGAU(IC)
  440.          CALL COQ8JC(IC,NBN1,E,XE,TH,TXR,SHPTOT,XJ,DJAC,IRR)
  441.                   ENDIF
  442.                   MELVAL = IVAL(1)
  443.                   IBMN=MIN(IB, MELVAL.VELCHE(/2))
  444.                   IGMN=MIN(IC, MELVAL.VELCHE(/1))
  445.                   MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(1)
  446.                   MELVAL = IVAL(2)
  447.                   IBMN=MIN(IB, MELVAL.VELCHE(/2))
  448.                   IGMN=MIN(IC, MELVAL.VELCHE(/1))
  449.                   MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(2)
  450.                   IF (IDIM .EQ. 3) THEN
  451.                      MELVAL = IVAL(3)
  452.                      IBMN=MIN(IB, MELVAL.VELCHE(/2))
  453.                      IGMN=MIN(IC, MELVAL.VELCHE(/1))
  454.                      MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(3)
  455.                   ENDIF
  456.  5002          CONTINUE
  457.  5000       CONTINUE
  458.             SEGSUP TR1
  459.             GOTO 520
  460.          ENDIF
  461. C---------------------------------------------------------------------
  462. C  DESACTIVATION DES SEGMENTS PROPRES A LA ZONE GEOMETRIQUE ISOUS
  463. C---------------------------------------------------------------------
  464.  520     CONTINUE
  465.          MPTVAL=IVAJAC
  466.          SEGSUP MPTVAL
  467.          SEGSUP TRA
  468.  
  469.   500 CONTINUE
  470.       RETURN
  471. C
  472. C-------------------------------------------------------------------
  473. C         ERREUR DANS UNE ZONE , DESACTIVATION ET RETOUR
  474. C-------------------------------------------------------------------
  475.  9990 CONTINUE
  476.       IRET = 0
  477.       MPTVAL=IVAJAC
  478.       SEGSUP MPTVAL
  479.  
  480. *      RETURN
  481.       END
  482.  
  483.  

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