Télécharger la plaquette Cast3M


fr en

Télécharger orient.eso

Retour à la liste

Numérotation des lignes : 

orient

  1. C ORIENT    SOURCE    PV        20/04/14    21:15:04     10584          
  2. C   CET OPERATEUR ORIENTE LES ELEMENTS ORIENTABLES EN FONCTION DE XP
  3. C     
  4.       SUBROUTINE ORIENT(IPT1,IPT2,XP,ICLE)
  5. *     
  6. C   IPT1     (E)   MAILLAGE A ORIENTER (segment ACTIF)
  7. C   IPT2     (S)   MAILLAGE ORIENTE    (segment ACTIF)
  8. C     IPT2 peut être égal à IPT1 si on ne savait pas orienter ou si le
  9. C     maillage
  10. C     est deja orienté...
  11. *   
  12. *   ICLE=1 ORIEN DIRE
  13. *   ICLE=2 ORIEN POINT
  14. *
  15. C SG : 2016/05/17 ajout orientation elements massifs
  16. *
  17. *     
  18.       IMPLICIT INTEGER(I-N)
  19.       IMPLICIT REAL*8 (A-H,O-z)
  20.  
  21. -INC CCREEL
  22. -INC SMELEME
  23. -INC SMLENTI
  24. -INC CCGEOME
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC SMCOORD
  29.  
  30.       DIMENSION XP(3)
  31.       DIMENSION XQ(3)
  32.       LOGICAL LQUAF
  33.       PARAMETER (NLMASS=22)
  34. *    NLNOMA=2*2 + 8*3 + 12*4
  35.       PARAMETER (NLNOMA=76)
  36. *     Tableau donnant la liste des types d'elements susceptibles d'etre
  37. *     traites dans le cas massif
  38.       INTEGER LTYMAS(NLMASS)
  39. *     Pour chaque element susceptible d'etre traite, ce tableau donne
  40. *     l'adresse dans le tableau LNOMAS de la description des noeuds 
  41. *     servant à calculer l'orientation
  42.       INTEGER LADMAS(NLMASS)
  43. *     Ce tableau donne les idim+1 noeuds (triedre) de chaque element massif
  44. *     servant pour calculer l'orientation
  45. *     On aurait pu le simplifier sachant que les QUAD et QUAF sont
  46. *     similaires...
  47.       INTEGER LNOMAS(NLNOMA)
  48. *                  SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8 QUA9   10
  49.       DATA LTYMAS/ 2,     3,  4,   5,   6,   7,   8,   9,   10,  11,    
  50. C                  CUB8 CU20  PRI6  PR15  TET4 TE10  PYR5  PY13        18 
  51.      $             14,   15,   16,   17,   23,  24,   25,   26,   
  52. C                  CU27  PR21  TE15  PY19                              22
  53.      $             33,   34,   35,   36/
  54. *                  SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8 QUA9   10
  55.       DATA LADMAS/ 1,     3,   5,   8,   11,  14,  17,  20,  23,  26, 
  56. C                  CUB8 CU20  PRI6  PR15  TET4 TE10  PYR5  PY13        18 
  57.      $             29,    33,   37,  41,    45,  49,   53,   57,   
  58. C                  CU27  PR21  TE15  PY19                              22
  59.      $             61,    65,    69,   73/
  60. *                  SEG2   SEG3 TRI3    TRI4    TRI6   TRI7   QUA4    
  61.       DATA LNOMAS/ 1,2,   1,3,  1,2,3,  1,2,3, 1,3,5, 1,3,5,  1,2,4, 
  62. C                  QUA5   QUA8    QUA9    CUB8     CU20       PRI6    
  63.      $             1,2,4, 1,3,7,  1,3,7,  1,2,4,5,  1,3,7,13, 1,2,3,4, 
  64. C                  PR15      TET4     TE10      PYR5     PY13        
  65.      $             1,3,5,10, 1,2,3,4, 1,3,5,10, 1,2,4,5,  1,3,7,13,                  
  66. C                  CU27      PR21       TE15       PY19  
  67.      $             1,3,7,13, 1,3,5,10,  1,3,5,10,  1,3,7,13/
  68.       segact mcoord
  69. *
  70. * Cas de l'orientation des éléments massifs
  71. *     On souhaite orienter positivement les éléments (ie même
  72. *     orientation que l'élément de référence)
  73. *
  74. * Est-on dans le cas des elements massifs ?
  75.       ilmass=0
  76.       ITY=IPT1.ITYPEL
  77.       DO i=1,nlmass
  78.          if (ity.eq.ltymas(i)) then
  79.             ilmass=i
  80.             goto 666
  81.          endif
  82.       enddo
  83.  666  continue
  84. *dbg      write(ioimp,*) 'ilmass=',ilmass
  85. *dbg      write(ioimp,*) 'idim=',idim
  86.       if (ilmass.eq.0) goto 1000
  87.       if (LDLR(ity).ne.idim) goto 1000
  88. *dbg      write(ioimp,*) 'Cas massif'
  89. *     Si oui, il faut que ICLE=0 sauf pour le cub8 (shb8) qui est traité
  90. *     plus loin
  91.       if (icle.ne.0) then
  92.          if (ity.eq.14) goto 1000
  93.          if (icle.eq.1) moterr(1:8)='DIRE    '
  94.          if (icle.eq.2) moterr(1:8)='POIN    '
  95. *dbg         write(ioimp,*) 'icle=',icle
  96. * 803 2
  97. * Option %m1:8 incompatible avec les donnees
  98.          call erreur(803)
  99.          return
  100.       endif
  101. *     Ici, on teste l'orientation des éléments et on la stocke dans une
  102. *     liste d'entiers valant +1 ou -1. 
  103. *     On stocke le nombre d'orientation négative dans ninv
  104. *
  105. * oubli des references (suivant la logique de invers)
  106. *eff      nbref=0
  107. *eff      nbsous=0
  108.       nbnn=ipt1.num(/1)
  109.       nbelem=ipt1.num(/2)
  110. *eff      segini meleme
  111. *eff      itypel=ipt1.itypel
  112.       iadr=ladmas(ilmass)-1
  113. *dbg      write(ioimp,*) 'iadr=',iadr
  114.       jg=nbelem
  115.       segini mlenti
  116.       ninv=0
  117.       do il=1,nbelem
  118. * test orientation
  119.          if (idim.eq.1) then
  120.             ia=ipt1.num(lnomas(iadr+1),il)
  121.             ib=ipt1.num(lnomas(iadr+2),il)
  122. *dbg            write(ioimp,*) 'ia=',ia,' ib=',ib
  123.             ira=(idim+1)*(ia-1)
  124.             irb=(idim+1)*(ib-1)
  125.             x1=xcoor(irb+1)-xcoor(ira+1)
  126.             pmix=x1
  127.          elseif (idim.eq.2) then
  128.             ia=ipt1.num(lnomas(iadr+1),il)
  129.             ib=ipt1.num(lnomas(iadr+2),il)
  130.             ic=ipt1.num(lnomas(iadr+3),il)
  131.             ira=(idim+1)*(ia-1)
  132.             irb=(idim+1)*(ib-1)
  133.             irc=(idim+1)*(ic-1)
  134.             x1=xcoor(irb+1)-xcoor(ira+1)
  135.             y1=xcoor(irb+2)-xcoor(ira+2)
  136.             x2=xcoor(irc+1)-xcoor(ira+1)
  137.             y2=xcoor(irc+2)-xcoor(ira+2)
  138.             pmix=x1*y2-y1*x2
  139.          elseif (idim.eq.3) then
  140. *dbg            write(ioimp,*) 'lnomas1=',lnomas(iadr+1)
  141.             ia=ipt1.num(lnomas(iadr+1),il)
  142.             ib=ipt1.num(lnomas(iadr+2),il)
  143.             ic=ipt1.num(lnomas(iadr+3),il)
  144.             id=ipt1.num(lnomas(iadr+4),il)
  145. *dbg            write(ioimp,*) 'ia=',ia,' ib=',ib,' ic=',ic,' id=',id
  146.             ira=(idim+1)*(ia-1)
  147.             irb=(idim+1)*(ib-1)
  148.             irc=(idim+1)*(ic-1)
  149.             ird=(idim+1)*(id-1)
  150.             x1=xcoor(irb+1)-xcoor(ira+1)
  151.             y1=xcoor(irb+2)-xcoor(ira+2)
  152.             z1=xcoor(irb+3)-xcoor(ira+3)
  153.             x2=xcoor(irc+1)-xcoor(ira+1)
  154.             y2=xcoor(irc+2)-xcoor(ira+2)
  155.             z2=xcoor(irc+3)-xcoor(ira+3)
  156.             x3=xcoor(ird+1)-xcoor(ira+1)
  157.             y3=xcoor(ird+2)-xcoor(ira+2)
  158.             z3=xcoor(ird+3)-xcoor(ira+3)
  159.             pmix=x1*(y2*z3-y3*z2)+y1*(z2*x3-z3*x2)+z1*(x2*y3-x3*y2)
  160.          endif
  161. *dbg         write(ioimp,*) 'pmix=',pmix
  162.          if (pmix.ge.0.d0) then
  163.             lect(il)=+1
  164.          else 
  165.             ninv=ninv+1
  166.             lect(il)=-1
  167.          endif
  168.       enddo
  169.       if (ninv.gt.0) then
  170.          call inver4(ipt1,ipt2,3,mlenti)
  171.          if (ierr.ne.0) return
  172.       else
  173.          ipt2=ipt1
  174.       endif
  175.       segsup mlenti
  176. * fin du cas massif
  177.       return
  178. *
  179. * Cas de l'orientation des éléments non massifs (sauf cub8 shb8)
  180. *
  181. *
  182.  1000 continue
  183.       IF (IDIM.EQ.3) THEN
  184.          IF (ICLE.NE.1.AND.ICLE.NE.2) CALL ERREUR(1074)
  185.       ELSE
  186. * SG Ceci est un residu de l'ancienne programmation en 2D, 
  187. * ou on oriente par rapport a un vecteur (0. 0. 1.). En fait, comme
  188. * on n'orientait que les TRI et les QUA, ceux-ci ont deja ete
  189. * traites dans le cas massif. Du coup, on ne passera normalement ici
  190. * que pour les elements que l'on ne sait pas orienter (SEG en 2D par ex)
  191. * Voir IPT2=IPT1 plus loin
  192.          ICLE=1
  193.       ENDIF
  194.  
  195.       IF (ICLE.EQ.1) THEN
  196. *  NORMALISATION DU VECTEUR
  197.          DP=SQRT(XP(1)**2+XP(2)**2+XP(3)**2)
  198.          IF (DP.LE.XPETIT) THEN
  199.             CALL ERREUR(277)
  200.             DP=1.D0
  201.          ENDIF
  202.          XQ(1)=XP(1)/DP
  203.          XQ(2)=XP(2)/DP
  204.          XQ(3)=XP(3)/DP
  205.       ENDIF
  206.       NBREF=IPT1.LISREF(/1)
  207.       NBSOUS=0
  208.       NBELEM=IPT1.NUM(/2)
  209.       NBNN=IPT1.NUM(/1)
  210.       ITYP=KSURF(IPT1.ITYPEL)
  211.       IF (ITYP.NE.0.AND.ITYP.EQ.IPT1.ITYPEL) GOTO 1
  212. *
  213. * ce n'est pas des éléments de surface
  214. * est-on en présence de cub8 (shb8)
  215. *     
  216.       if(ipt1.itypel.eq.14) then
  217.          nbref=2
  218.          segini ipt2
  219.          nbnn=4
  220.          nbref=0
  221.          segini ipt3,ipt4
  222.          ipt2.lisref(1)=ipt3
  223.          ipt2.lisref(2)=ipt4
  224.          ipt3.itypel=8
  225.          ipt4.itypel=8
  226.          ipt2.itypel=14
  227.          idim1=idim+1
  228.          do i=1,ipt1.num(/2)
  229.             ia=ipt1.num(1,i) - 1
  230.             ib=ipt1.num(5,i) - 1
  231.             ic=ipt1.num(2,i) - 1
  232.             id=ipt1.num(3,i) - 1
  233.             xab=xcoor(ib*idim1+1)-xcoor(ia*idim1+1)
  234.             yab=xcoor(ib*idim1+2)-xcoor(ia*idim1+2)
  235.             zab=xcoor(ib*idim1+3)-xcoor(ia*idim1+3)
  236.             xac=xcoor(ic*idim1+1)-xcoor(ia*idim1+1)
  237.             yac=xcoor(ic*idim1+2)-xcoor(ia*idim1+2)
  238.             zac=xcoor(ic*idim1+3)-xcoor(ia*idim1+3)
  239.             xad=xcoor(id*idim1+1)-xcoor(ia*idim1+1)
  240.             yad=xcoor(id*idim1+2)-xcoor(ia*idim1+2)
  241.             zad=xcoor(id*idim1+3)-xcoor(ia*idim1+3)
  242.             xpvec= yac*zad - zac*yad
  243.             ypvec= zac*xad - xac*zad
  244.             zpvec= xac*yad - yac*xad
  245.             pmix=xpvec*xab+ypvec*yab+zpvec*zab
  246.             isens=0
  247.             if(pmix.ge.0.d0) isens=1
  248.             if(icle.eq.1) then
  249.                xsc=xab*xq(1)+yab*xq(2)+zab*xq(3)
  250.             else
  251.                xsc= xab*(xq(1)-xcoor(ib*idim1+1))+
  252.      $              yab*(xq(2)-xcoor(ib*idim1+2))+
  253.      $              zab*(xq(3)-xcoor(ib*idim1+3))
  254.             endif
  255.             if(xsc.ge.0.d0) then
  256.                do j=1,8
  257.                   ipt2.num(j,i)=ipt1.num(j,i)
  258.                enddo
  259.                if(isens.eq.2) then
  260.                   iu=ipt2.num(2,i)
  261.                   ipt2.num(2,i)=ipt2.num(4,i)
  262.                   ipt2.num(4,i)=iu
  263.                   iu=ipt2.num(6,i)
  264.                   ipt2.num(6,i)=ipt2.num(8,i)
  265.                   ipt2.num(8,i)=iu
  266.                endif
  267.             else
  268.                ipt2.num(1,i)=ipt1.num(5,i)
  269.                ipt2.num(2,i)=ipt1.num(8,i)
  270.                ipt2.num(3,i)=ipt1.num(7,i)
  271.                ipt2.num(4,i)=ipt1.num(6,i)
  272.                ipt2.num(5,i)=ipt1.num(1,i)
  273.                ipt2.num(6,i)=ipt1.num(4,i)
  274.                ipt2.num(7,i)=ipt1.num(3,i)
  275.                ipt2.num(8,i)=ipt1.num(2,i)
  276.                if(isens.eq.2) then
  277.                   iu=ipt2.num(2,i)
  278.                   ipt2.num(2,i)=ipt2.num(4,i)
  279.                   ipt2.num(4,i)=iu
  280.                   iu=ipt2.num(6,i)
  281.                   ipt2.num(6,i)=ipt2.num(8,i)
  282.                   ipt2.num(8,i)=iu
  283.                endif
  284.             endif
  285.             ipt2.icolor(i)=ipt1.icolor(i)
  286.             ipt3.num(1,i)=ipt2.num(1,i)
  287.             ipt3.num(2,i)=ipt2.num(2,i)
  288.             ipt3.num(3,i)=ipt2.num(3,i)
  289.             ipt3.num(4,i)=ipt2.num(4,i)
  290.             ipt3.icolor(i)=ipt2.icolor(i)
  291.             ipt4.num(1,i)=ipt2.num(5,i)
  292.             ipt4.num(2,i)=ipt2.num(6,i)
  293.             ipt4.num(3,i)=ipt2.num(7,i)
  294.             ipt4.num(4,i)=ipt2.num(8,i)
  295.             ipt4.icolor(i)=ipt2.icolor(i)
  296.          enddo
  297.          segdes ipt3,ipt4
  298.          return
  299.       endif
  300. *
  301. * Sinon, on ne savait pas orienter les éléments de ipt1
  302. *
  303.       IPT2=IPT1
  304.       RETURN
  305. *     
  306. *  on a des vrais éléments de surfaces
  307. *     
  308.  1    CONTINUE
  309. * Cas QUAF TRI7 ou QUA9
  310.       LQUAF=(ITYP.EQ.11.OR.ITYP.EQ.7)
  311.       SEGINI IPT2
  312.       IPT2.ITYPEL=ITYP
  313.       IF (NBREF.EQ.0) GOTO 3
  314.       DO 2 I=1,NBREF
  315.          IPT2.LISREF(I)=IPT1.LISREF(I)
  316.  2    CONTINUE
  317.  3    CONTINUE
  318.       SEGACT MCOORD
  319.       IB=NSPOS(ITYP)-1
  320.       IS1=IBSOM(IB+1)
  321.       IS2=IBSOM(IB+2)
  322.       IS3=IBSOM(IB+3)
  323.  
  324.  
  325.       DO 4 J=1,NBELEM
  326.          IP1=IPT1.NUM(IS1,J)
  327.          IP2=IPT1.NUM(IS2,J)
  328.          IP3=IPT1.NUM(IS3,J)
  329.          IREF1=(IDIM+1)*(IP1-1)
  330.          IREF2=(IDIM+1)*(IP2-1)
  331.          IREF3=(IDIM+1)*(IP3-1)
  332.          XV1=XCOOR(IREF2+1)-XCOOR(IREF1+1)
  333.          YV1=XCOOR(IREF2+2)-XCOOR(IREF1+2)
  334.          ZV1=XCOOR(IREF2+3)-XCOOR(IREF1+3)
  335.          IF (IDIM.NE.3) ZV1=0.D0
  336.          XV2=XCOOR(IREF2+1)-XCOOR(IREF3+1)
  337.          YV2=XCOOR(IREF2+2)-XCOOR(IREF3+2)
  338.          ZV2=XCOOR(IREF2+3)-XCOOR(IREF3+3)
  339.          IF (IDIM.NE.3) ZV2=0.D0
  340.          XV3=ZV1*YV2-ZV2*YV1
  341.          YV3=XV1*ZV2-XV2*ZV1
  342.          ZV3=YV1*XV2-YV2*XV1
  343.          XVN=SQRT(XV3**2+YV3**2+ZV3**2)
  344.          IF (ICLE.EQ.2) THEN
  345.             XQ(1)=0.D0
  346.             XQ(2)=0.D0
  347.             XQ(3)=0.D0
  348.             DO 10 L=1,NBNN
  349.                IREF=(IPT1.NUM(L,J)-1)*(IDIM+1)
  350.                DO 11 K=1,3
  351.                   XQ(K)=XQ(K)+XCOOR(IREF+K)
  352.  11            CONTINUE
  353.  10         CONTINUE
  354.             XQ(1)=XQ(1)/NBNN
  355.             XQ(2)=XQ(2)/NBNN
  356.             XQ(3)=XQ(3)/NBNN
  357.             XQ(1)=XP(1)-XQ(1)
  358.             XQ(2)=XP(2)-XQ(2)
  359.             XQ(3)=XP(3)-XQ(3)
  360.             DP=SQRT(XQ(1)**2+XQ(2)**2+XQ(3)**2)
  361.             IF (DP.LE.XPETIT) THEN
  362.                CALL ERREUR(277)
  363.                DP=1.D0
  364.             ENDIF
  365.             XQ(1)=XQ(1)/DP
  366.             XQ(2)=XQ(2)/DP
  367.             XQ(3)=XQ(3)/DP
  368.          ENDIF
  369.          TEST=XV3*XQ(1)+YV3*XQ(2)+ZV3*XQ(3)
  370.          IF (ABS(TEST).LE.0.0175D0*XVN) CALL ERREUR(278)
  371.          IF (TEST.GE.0.D0) THEN
  372.             DO 6 I=1,NBNN
  373.                IPT2.NUM(I,J)=IPT1.NUM(I,J)
  374.  6          CONTINUE
  375.          ELSE
  376.             IF (.NOT.LQUAF) THEN
  377.                DO 7 I=1,NBNN
  378.                   IPT2.NUM(MOD(NBNN+1-I,NBNN)+1,J)=IPT1.NUM(I,J)
  379.  7             CONTINUE
  380.             ELSE
  381.                NBNN1=NBNN-1
  382.                DO 8 I=1,NBNN1
  383.                   IPT2.NUM(MOD(NBNN1+1-I,NBNN1)+1,J)=IPT1.NUM(I,J)
  384.  8             CONTINUE
  385.                IPT2.NUM(NBNN,J)=IPT1.NUM(NBNN,J)
  386.             ENDIF
  387.          ENDIF
  388.          IPT2.ICOLOR(J)=IPT1.ICOLOR(J)
  389.  4    CONTINUE
  390.       RETURN
  391.       END
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  

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