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/01/09    21:15:37     12442          
  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.          IF (MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9 .OR. MFR.EQ.77) THEN
  135.             N1PTEL=NBPGAU
  136.             N1EL=NBELEM
  137.          ELSEIF(MFR.EQ.7 .OR. MFR.EQ.13) THEN
  138.             N1PTEL=NBPGAU
  139.             N1EL=NBELEM
  140.          ELSE
  141.             N1PTEL = 0
  142.             N1EL = 0
  143.          ENDIF
  144.          N2PTEL=0
  145.          N2EL  =0
  146. C
  147. C     CREATION DU MCHAML DE LA SOUS ZONE
  148. C
  149.          NJAC=IDIM
  150.          N2 = NJAC
  151.          SEGINI MCHAML
  152.          ICHAML(ISOUS)=MCHAML
  153.          NSR   =1
  154.          NCOSOR=NJAC
  155.          SEGINI MPTVAL
  156.          IVAJAC=MPTVAL
  157. C
  158. C     2 OU 3 COMPOSANTES
  159. C
  160.          ICOMP=1
  161.          IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN
  162.             NOMCHE(ICOMP)='VR  '
  163.          ELSE
  164.             NOMCHE(ICOMP)='VX  '
  165.          ENDIF
  166.          TYPCHE(ICOMP)='REAL*8'
  167.          SEGINI MELVA1
  168.          IELVAL(ICOMP)=MELVA1
  169.          IVAL(ICOMP)=MELVA1
  170. C
  171.          ICOMP=2
  172.          IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN
  173.             NOMCHE(ICOMP)='VZ  '
  174.          ELSE
  175.             NOMCHE(ICOMP)='VY  '
  176.          ENDIF
  177.          TYPCHE(ICOMP)='REAL*8'
  178.          SEGINI MELVA2
  179.          IELVAL(ICOMP)=MELVA2
  180.          IVAL(ICOMP)=MELVA2
  181. C
  182.          MELVA3 = 0
  183.          IF (IDIM .EQ. 3) THEN
  184.             ICOMP=3
  185.             NOMCHE(ICOMP)='VZ  '
  186.             TYPCHE(ICOMP)='REAL*8'
  187.             SEGINI MELVA3
  188.             IELVAL(ICOMP)=MELVA3
  189.             IVAL(ICOMP)=MELVA3
  190.          ENDIF
  191. C
  192.          SEGINI TRA
  193. C
  194. C ================ FORMULATION MASSIVE =======================
  195. C
  196.          IF(MFR.EQ.1.OR.MFR.EQ.33) THEN
  197.             GOTO 520
  198. C
  199. C ================ FORMULATION COQUE  MINCE =====================
  200. C
  201. C
  202.          ELSE IF(MFR.EQ.3.OR.MFR.EQ.9 .OR. MFR.EQ.77) THEN
  203.           IDI2=IDIM-1
  204.           DO 3000 IB=1,NBELEM
  205. *--------------Calcul de la normale a l'élément
  206.             INU1=MELEME.NUM(1,IB) 
  207.             INU2=MELEME.NUM(2,IB) 
  208.             IREF1 = IDIMP1*(INU1 - 1)
  209.             IREF2 = IDIMP1*(INU2 - 1)
  210.             XNORU = 0.D0
  211.             DO IC = 1, IDIM
  212.                r_z = XCOOR(IREF2+IC)-XCOOR(IREF1+IC)
  213.                XU(IC) = r_z
  214.                XNORU = XNORU + (r_z * r_z)
  215.             ENDDO
  216.             XNORU = SQRT(XNORU)
  217.             DO IC = 1, IDIM
  218.                XU(IC) = XU(IC)/XNORU
  219.             ENDDO
  220.             IF (IDIM .EQ. 2) THEN
  221.                XW(1) = -XU(2)
  222.                XW(2) = XU(1)
  223.             ELSE
  224.                IN = 3
  225.  33            CONTINUE
  226.                INUIN=MELEME.NUM(IN,IB)
  227.                IREF3 = IDIMP1*(INUIN - 1)
  228.                XNORV = 0.D0
  229.                DO IC = 1, IDIM
  230.                   r_z = XCOOR(IREF3+IC)-XCOOR(IREF1+IC)
  231.                   XV(IC) = r_z
  232.                   XNORV = XNORV + (r_z * r_z)
  233.                ENDDO
  234.                XNORV = SQRT(XNORV)
  235.                DO IC = 1, IDIM
  236.                  XV(IC) = XV(IC)/XNORV
  237.                ENDDO
  238.                XW(1) = XU(2)*XV(3) - XU(3)*XV(2)
  239.                XW(2) = XU(3)*XV(1) - XU(1)*XV(3)
  240.                XW(3) = XU(1)*XV(2) - XU(2)*XV(1)
  241.                XNORW = 0.
  242.                DO IC = 1, IDIM
  243.                  XNORW = XNORW + XW(IC)*XW(IC)
  244.                ENDDO
  245.                IF (XNORW .LE. XTOL) THEN
  246.                   write(ioimp,*) '  IN=',IN
  247.                   write(ioimp,*) '  W1,W2,W3=',XW(1),XW(2),XW(3)
  248.                   write(ioimp,*) '  XNORW2=',XNORW
  249.                   IN = IN + 1
  250.                   if(IN.le.NBNN) GOTO 33
  251.                ENDIF
  252.                XNORW = SQRT(XNORW)
  253.                IF(XNORW .LE.XTOL) then
  254.                   write(ioimp,*) 'INU1,INU2,INUIN=',INU1,INU2,INUIN
  255.                   write(ioimp,*) 'U1,U2,U3=',XU(1),XU(2),XU(3)
  256.                   write(ioimp,*) 'V1,V2,V3=',XV(1),XV(2),XV(3)
  257.                   write(ioimp,*) 'W1,W2,W3=',XW(1),XW(2),XW(3)
  258.                   write(ioimp,*) 'XNORW=',XNORW
  259.                   call erreur(345)
  260.                   return
  261.                ENDIF
  262.                DO IC = 1, IDIM
  263.                   XW(IC) = XW(IC)/XNORW
  264.                ENDDO
  265.             ENDIF
  266. *--------------Fin du calcul de la normale a l'élément
  267.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  268. C
  269.             IF(IDIM.EQ.2)THEN
  270.                CALL VPAST2(XE,BPSS)
  271.             ELSE IF(IDIM.EQ.3) THEN
  272.                CALL VPAST(XE,BPSS)
  273.             ENDIF
  274.             CALL VCORL1(XE,XEL,BPSS,NBNN)
  275.             MPTVAL=IVAJAC
  276.             DO 3002 IC=1,NBPGAU
  277.                IF (INORM .EQ. 1) THEN
  278.                   DJAC = 1.D0
  279.                ELSE
  280.                   DO IE=1,NBNN
  281.                      DO ID=1,6
  282.                         SHP(ID,IE)=SHPTOT(ID,IE,IC)
  283.                      ENDDO
  284.                   ENDDO
  285.                   CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC)
  286.                ENDIF
  287.                MELVAL = IVAL(1)
  288.                IBMN=MIN(IB, MELVAL.VELCHE(/2))
  289.                IGMN=MIN(IC,MELVAL.VELCHE(/1))
  290.                MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(1)
  291.                MELVAL = IVAL(2)
  292.                IBMN=MIN(IB, MELVAL.VELCHE(/2))
  293.                IGMN=MIN(IC,MELVAL.VELCHE(/1))
  294.                MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(2)
  295.                IF (IDIM .EQ. 3) THEN
  296.                   MELVAL = IVAL(3)
  297.                   IBMN=MIN(IB, MELVAL.VELCHE(/2))
  298.                   IGMN=MIN(IC,MELVAL.VELCHE(/1))
  299.                   MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(3)
  300.                ENDIF
  301.  3002       CONTINUE
  302.  3000    CONTINUE
  303.          GOTO 520
  304. C
  305. C ================ FORMULATION POUTRE ET TUYAU ====================
  306. C
  307.       ELSE IF(MFR.EQ.7.OR.MFR.EQ.13) THEN
  308.          IDI2=IDIM-1
  309.          DO 4000 IB=1,NBELEM
  310. *-----------Calcul de la tangente a l'élément
  311.             IREF1 = IDIMP1*(MELEME.NUM(1, IB) - 1)
  312.             IREF2 = IDIMP1*(MELEME.NUM(2, IB) - 1)
  313.             XNORU = 0.D0
  314.             DO IC = 1, IDIM
  315.                r_z = XCOOR(IREF2+IC)-XCOOR(IREF1+IC)
  316.                XU(IC) = r_z
  317.                XNORU = XNORU + (r_z * r_z)
  318.             ENDDO
  319.             XNORU = SQRT(XNORU)
  320.             DO IC = 1, IDIM
  321.                XU(IC) = XU(IC)/XNORU
  322.             ENDDO
  323. *-----------Fin du calcul de la tangente a l'élément
  324. *-----------Calcul de la tangente a l'élément en chaque point de Gauss
  325. c  BP : On suppose le jacobien constant dans l'element (idem POUJAC.eso)
  326. C  => on sort le calcul du jacobien de la boucle sur les points de Gauss.
  327. C  Cela implique que la POUTre de Bernoulli n'est pas isoparamétrique...
  328.             IF (INORM .EQ. 1) THEN
  329.                DJAC = 1.D0
  330.             ELSE
  331.                DJAC = 1.D0/DBLE(NBPGAU)
  332.             ENDIF
  333.             MPTVAL=IVAJAC
  334.             DO 4002 IC=1,NBPGAU
  335.                MELVAL = IVAL(1)
  336.                IBMN=MIN(IB, MELVAL.VELCHE(/2))
  337.                IGMN=MIN(IC,MELVAL.VELCHE(/1))
  338.                MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(1)
  339.                MELVAL = IVAL(2)
  340.                IBMN=MIN(IB, MELVAL.VELCHE(/2))
  341.                IGMN=MIN(IC,MELVAL.VELCHE(/1))
  342.                MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(2)
  343.                IF (IDIM .EQ. 3) THEN
  344.                   MELVAL = IVAL(3)
  345.                   IBMN=MIN(IB, MELVAL.VELCHE(/2))
  346.                   IGMN=MIN(IC,MELVAL.VELCHE(/1))
  347.                   MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(3)
  348.                ENDIF
  349.  4002       CONTINUE
  350.  4000    CONTINUE
  351.  
  352.          GOTO 520
  353. C
  354. C ================ FORMULATION COQUE  EPAISSE ====================
  355. C
  356.          ELSE IF(MFR.EQ.5) THEN
  357.             SEGACT MINTE1
  358.             NBPGA1=MINTE1.POIGAU(/1)
  359.             NBN1  =MINTE1.SHPTOT(/2)
  360.             SEGINI TR1
  361. C
  362. C     UNE PETITE HORREUR ON CONSIDERE LES EPAISSEURS CONSTANTES
  363. C
  364.             DO 5010 IC=1,NBNN
  365.                TH(IC)=UN
  366.  5010       CONTINUE
  367.             DO 5000 IB=1,NBELEM
  368. *--------------Calcul de la normale a l'élément
  369.                IREF1 = IDIMP1*(MELEME.NUM(1, IB) - 1)
  370. *               IREF2 = IDIMP1*(MELEME.NUM(2, IB) - 1)
  371. * bp : les EF de coque epaisse etant quadratiques (coq6 et coq8), on
  372. * prend les noeuds "coins" pour eviter pb avec les noeuds 1,2,3 si courbures
  373.                IREF2 = IDIMP1*(MELEME.NUM(3, IB) - 1)
  374.                XNORU = 0.
  375.                DO IC = 1, IDIM
  376.                   r_z = XCOOR(IREF2+IC)-XCOOR(IREF1+IC)
  377.                   XU(IC) = r_z
  378.                   XNORU = XNORU + (r_z * r_z)
  379.                ENDDO
  380.                XNORU = SQRT(XNORU)
  381.                DO IC = 1, IDIM
  382.                   XU(IC) = XU(IC)/XNORU
  383.                ENDDO
  384.                IF (IDIM .EQ. 2) THEN
  385.                   XW(1) = -XU(2)
  386.                   XW(2) = XU(1)
  387.                ELSE
  388. *                  IN = 3
  389.                   IN = 5
  390.    13             IREF3 = IDIMP1*(MELEME.NUM(IN, IB) - 1)
  391.                   XNORV = 0.D0
  392.                   DO IC = 1, IDIM
  393.                     r_z = XCOOR(IREF3 + IC) - XCOOR(IREF1 + IC)
  394.                     XV(IC) = r_z
  395.                     XNORV = XNORV + (r_z * r_z)
  396.                   ENDDO
  397.                   XNORV = SQRT(XNORV)
  398.                   DO IC = 1, IDIM
  399.                     XV(IC) = XV(IC)/XNORV
  400.                   ENDDO
  401.                   XW(1) = XU(2)*XV(3) - XU(3)*XV(2)
  402.                   XW(2) = XU(3)*XV(1) - XU(1)*XV(3)
  403.                   XW(3) = XU(1)*XV(2) - XU(2)*XV(1)
  404.                   XNORW = 0.
  405.                   DO IC = 1, IDIM
  406.                      XNORW = XNORW + XW(IC)*XW(IC)
  407.                   ENDDO
  408.                   XNORW = SQRT(XNORW)
  409.                   IF (XNORW .LT. 1.E-4) THEN
  410.                    if(IN.LT.NBNN) then
  411.                      IN = IN + 1
  412.                      GOTO 13
  413.                    else
  414.                     write(6,*) 'Difficultes pour etablir la normale de'
  415.                     write(6,*) 'l element',IB
  416.                     write(6,*) 'Verifiez votre maillage'
  417.                     GOTO 9990
  418.                    endif
  419.                   ENDIF
  420.                   DO IC = 1, IDIM
  421.                      XW(IC) = XW(IC)/XNORW
  422.                   ENDDO
  423.                ENDIF
  424. *--------------Fin du calcul de la normale a l'élément
  425.                IF (INORM .EQ. 0) THEN
  426.                   CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  427. C
  428.                   CALL CQ8LOC(XE,NBN1,MINTE1.SHPTOT,TXR,IRR)
  429.                ENDIF
  430. C
  431.                MPTVAL=IVAJAC
  432.                DO 5002 IC=1,NBPGAU
  433.                   IF (INORM .EQ. 1) THEN
  434.                      DJAC = 1.D0
  435.                   ELSE
  436.                      E=DZEGAU(IC)
  437.          CALL COQ8JC(IC,NBN1,E,XE,TH,TXR,SHPTOT,XJ,DJAC,IRR)
  438.                   ENDIF
  439.                   MELVAL = IVAL(1)
  440.                   IBMN=MIN(IB, MELVAL.VELCHE(/2))
  441.                   IGMN=MIN(IC, MELVAL.VELCHE(/1))
  442.                   MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(1)
  443.                   MELVAL = IVAL(2)
  444.                   IBMN=MIN(IB, MELVAL.VELCHE(/2))
  445.                   IGMN=MIN(IC, MELVAL.VELCHE(/1))
  446.                   MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(2)
  447.                   IF (IDIM .EQ. 3) THEN
  448.                      MELVAL = IVAL(3)
  449.                      IBMN=MIN(IB, MELVAL.VELCHE(/2))
  450.                      IGMN=MIN(IC, MELVAL.VELCHE(/1))
  451.                      MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(3)
  452.                   ENDIF
  453.  5002          CONTINUE
  454.  5000       CONTINUE
  455.             SEGSUP TR1
  456.             GOTO 520
  457.          ENDIF
  458. C---------------------------------------------------------------------
  459. C  DESACTIVATION DES SEGMENTS PROPRES A LA ZONE GEOMETRIQUE ISOUS
  460. C---------------------------------------------------------------------
  461.  520     CONTINUE
  462.          MPTVAL=IVAJAC
  463.          SEGSUP MPTVAL
  464.          SEGSUP TRA
  465.  
  466.   500 CONTINUE
  467.       RETURN
  468. C
  469. C-------------------------------------------------------------------
  470. C         ERREUR DANS UNE ZONE , DESACTIVATION ET RETOUR
  471. C-------------------------------------------------------------------
  472.  9990 CONTINUE
  473.       IRET = 0
  474.       MPTVAL=IVAJAC
  475.       SEGSUP MPTVAL
  476.  
  477. *      RETURN
  478.       END
  479.  
  480.  

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