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    PV090527  26/07/20    21:15:04     12601          
  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.          if (ierr.ne.0) return
  186.       ELSE
  187. * SG Ceci est un residu de l'ancienne programmation en 2D, 
  188. * ou on oriente par rapport a un vecteur (0. 0. 1.). En fait, comme
  189. * on n'orientait que les TRI et les QUA, ceux-ci ont deja ete
  190. * traites dans le cas massif. Du coup, on ne passera normalement ici
  191. * que pour les elements que l'on ne sait pas orienter (SEG en 2D par ex)
  192. * Voir IPT2=IPT1 plus loin
  193.          ICLE=1
  194.       ENDIF
  195.  
  196.       IF (ICLE.EQ.1) THEN
  197. *  NORMALISATION DU VECTEUR
  198.          DP=SQRT(XP(1)**2+XP(2)**2+XP(3)**2)
  199.          IF (DP.LE.XPETIT) THEN
  200.             CALL ERREUR(277)
  201.             DP=1.D0
  202.          ENDIF
  203.          XQ(1)=XP(1)/DP
  204.          XQ(2)=XP(2)/DP
  205.          XQ(3)=XP(3)/DP
  206.       ENDIF
  207.       NBREF=IPT1.LISREF(/1)
  208.       NBSOUS=0
  209.       NBELEM=IPT1.NUM(/2)
  210.       NBNN=IPT1.NUM(/1)
  211.       ITYP=KSURF(IPT1.ITYPEL)
  212.       IF (ITYP.NE.0.AND.ITYP.EQ.IPT1.ITYPEL) GOTO 1
  213. *
  214. * ce n'est pas des éléments de surface
  215. * est-on en présence de cub8 (shb8)
  216. *     
  217.       if(ipt1.itypel.eq.14) then
  218.          nbref=2
  219.          segini ipt2
  220.          nbnn=4
  221.          nbref=0
  222.          segini ipt3,ipt4
  223.          ipt2.lisref(1)=ipt3
  224.          ipt2.lisref(2)=ipt4
  225.          ipt3.itypel=8
  226.          ipt4.itypel=8
  227.          ipt2.itypel=14
  228.          idim1=idim+1
  229.          do i=1,ipt1.num(/2)
  230.             ia=ipt1.num(1,i) - 1
  231.             ib=ipt1.num(5,i) - 1
  232.             ic=ipt1.num(2,i) - 1
  233.             id=ipt1.num(3,i) - 1
  234.             xab=xcoor(ib*idim1+1)-xcoor(ia*idim1+1)
  235.             yab=xcoor(ib*idim1+2)-xcoor(ia*idim1+2)
  236.             zab=xcoor(ib*idim1+3)-xcoor(ia*idim1+3)
  237.             xac=xcoor(ic*idim1+1)-xcoor(ia*idim1+1)
  238.             yac=xcoor(ic*idim1+2)-xcoor(ia*idim1+2)
  239.             zac=xcoor(ic*idim1+3)-xcoor(ia*idim1+3)
  240.             xad=xcoor(id*idim1+1)-xcoor(ia*idim1+1)
  241.             yad=xcoor(id*idim1+2)-xcoor(ia*idim1+2)
  242.             zad=xcoor(id*idim1+3)-xcoor(ia*idim1+3)
  243.             xpvec= yac*zad - zac*yad
  244.             ypvec= zac*xad - xac*zad
  245.             zpvec= xac*yad - yac*xad
  246.             pmix=xpvec*xab+ypvec*yab+zpvec*zab
  247.             isens=0
  248.             if(pmix.ge.0.d0) isens=1
  249.             if(icle.eq.1) then
  250.                xsc=xab*xq(1)+yab*xq(2)+zab*xq(3)
  251.             else
  252.                xsc= xab*(xq(1)-xcoor(ib*idim1+1))+
  253.      $              yab*(xq(2)-xcoor(ib*idim1+2))+
  254.      $              zab*(xq(3)-xcoor(ib*idim1+3))
  255.             endif
  256.             if(xsc.ge.0.d0) then
  257.                do j=1,8
  258.                   ipt2.num(j,i)=ipt1.num(j,i)
  259.                enddo
  260.                if(isens.eq.2) then
  261.                   iu=ipt2.num(2,i)
  262.                   ipt2.num(2,i)=ipt2.num(4,i)
  263.                   ipt2.num(4,i)=iu
  264.                   iu=ipt2.num(6,i)
  265.                   ipt2.num(6,i)=ipt2.num(8,i)
  266.                   ipt2.num(8,i)=iu
  267.                endif
  268.             else
  269.                ipt2.num(1,i)=ipt1.num(5,i)
  270.                ipt2.num(2,i)=ipt1.num(8,i)
  271.                ipt2.num(3,i)=ipt1.num(7,i)
  272.                ipt2.num(4,i)=ipt1.num(6,i)
  273.                ipt2.num(5,i)=ipt1.num(1,i)
  274.                ipt2.num(6,i)=ipt1.num(4,i)
  275.                ipt2.num(7,i)=ipt1.num(3,i)
  276.                ipt2.num(8,i)=ipt1.num(2,i)
  277.                if(isens.eq.2) then
  278.                   iu=ipt2.num(2,i)
  279.                   ipt2.num(2,i)=ipt2.num(4,i)
  280.                   ipt2.num(4,i)=iu
  281.                   iu=ipt2.num(6,i)
  282.                   ipt2.num(6,i)=ipt2.num(8,i)
  283.                   ipt2.num(8,i)=iu
  284.                endif
  285.             endif
  286.             ipt2.icolor(i)=ipt1.icolor(i)
  287.             ipt3.num(1,i)=ipt2.num(1,i)
  288.             ipt3.num(2,i)=ipt2.num(2,i)
  289.             ipt3.num(3,i)=ipt2.num(3,i)
  290.             ipt3.num(4,i)=ipt2.num(4,i)
  291.             ipt3.icolor(i)=ipt2.icolor(i)
  292.             ipt4.num(1,i)=ipt2.num(5,i)
  293.             ipt4.num(2,i)=ipt2.num(6,i)
  294.             ipt4.num(3,i)=ipt2.num(7,i)
  295.             ipt4.num(4,i)=ipt2.num(8,i)
  296.             ipt4.icolor(i)=ipt2.icolor(i)
  297.          enddo
  298.          segdes ipt3,ipt4
  299.          return
  300.       endif
  301. *
  302. * Sinon, on ne savait pas orienter les éléments de ipt1
  303. *
  304.       IPT2=IPT1
  305.       RETURN
  306. *     
  307. *  on a des vrais éléments de surfaces
  308. *     
  309.  1    CONTINUE
  310. * Cas QUAF TRI7 ou QUA9
  311.       LQUAF=(ITYP.EQ.11.OR.ITYP.EQ.7)
  312.       SEGINI IPT2
  313.       IPT2.ITYPEL=ITYP
  314.       IF (NBREF.EQ.0) GOTO 3
  315.       DO 2 I=1,NBREF
  316.          IPT2.LISREF(I)=IPT1.LISREF(I)
  317.  2    CONTINUE
  318.  3    CONTINUE
  319.       SEGACT MCOORD
  320.       IB=NSPOS(ITYP)-1
  321.       IS1=IBSOM(IB+1)
  322.       IS2=IBSOM(IB+2)
  323.       IS3=IBSOM(IB+3)
  324.  
  325.  
  326.       DO 4 J=1,NBELEM
  327.          IP1=IPT1.NUM(IS1,J)
  328.          IP2=IPT1.NUM(IS2,J)
  329.          IP3=IPT1.NUM(IS3,J)
  330.          IREF1=(IDIM+1)*(IP1-1)
  331.          IREF2=(IDIM+1)*(IP2-1)
  332.          IREF3=(IDIM+1)*(IP3-1)
  333.          XV1=XCOOR(IREF2+1)-XCOOR(IREF1+1)
  334.          YV1=XCOOR(IREF2+2)-XCOOR(IREF1+2)
  335.          ZV1=XCOOR(IREF2+3)-XCOOR(IREF1+3)
  336.          IF (IDIM.NE.3) ZV1=0.D0
  337.          XV2=XCOOR(IREF2+1)-XCOOR(IREF3+1)
  338.          YV2=XCOOR(IREF2+2)-XCOOR(IREF3+2)
  339.          ZV2=XCOOR(IREF2+3)-XCOOR(IREF3+3)
  340.          IF (IDIM.NE.3) ZV2=0.D0
  341.          XV3=ZV1*YV2-ZV2*YV1
  342.          YV3=XV1*ZV2-XV2*ZV1
  343.          ZV3=YV1*XV2-YV2*XV1
  344.          XVN=SQRT(XV3**2+YV3**2+ZV3**2)
  345.          IF (ICLE.EQ.2) THEN
  346.             XQ(1)=0.D0
  347.             XQ(2)=0.D0
  348.             XQ(3)=0.D0
  349.             DO 10 L=1,NBNN
  350.                IREF=(IPT1.NUM(L,J)-1)*(IDIM+1)
  351.                DO 11 K=1,3
  352.                   XQ(K)=XQ(K)+XCOOR(IREF+K)
  353.  11            CONTINUE
  354.  10         CONTINUE
  355.             XQ(1)=XQ(1)/NBNN
  356.             XQ(2)=XQ(2)/NBNN
  357.             XQ(3)=XQ(3)/NBNN
  358.             XQ(1)=XP(1)-XQ(1)
  359.             XQ(2)=XP(2)-XQ(2)
  360.             XQ(3)=XP(3)-XQ(3)
  361.             DP=SQRT(XQ(1)**2+XQ(2)**2+XQ(3)**2)
  362.             IF (DP.LE.XPETIT) THEN
  363.                CALL ERREUR(277)
  364.                DP=1.D0
  365.             ENDIF
  366.             XQ(1)=XQ(1)/DP
  367.             XQ(2)=XQ(2)/DP
  368.             XQ(3)=XQ(3)/DP
  369.          ENDIF
  370.          TEST=XV3*XQ(1)+YV3*XQ(2)+ZV3*XQ(3)
  371.          IF (ABS(TEST).LE.0.0175D0*XVN) CALL ERREUR(278)
  372.          IF (TEST.GE.0.D0) THEN
  373.             DO 6 I=1,NBNN
  374.                IPT2.NUM(I,J)=IPT1.NUM(I,J)
  375.  6          CONTINUE
  376.          ELSE
  377.             IF (.NOT.LQUAF) THEN
  378.                DO 7 I=1,NBNN
  379.                   IPT2.NUM(MOD(NBNN+1-I,NBNN)+1,J)=IPT1.NUM(I,J)
  380.  7             CONTINUE
  381.             ELSE
  382.                NBNN1=NBNN-1
  383.                DO 8 I=1,NBNN1
  384.                   IPT2.NUM(MOD(NBNN1+1-I,NBNN1)+1,J)=IPT1.NUM(I,J)
  385.  8             CONTINUE
  386.                IPT2.NUM(NBNN,J)=IPT1.NUM(NBNN,J)
  387.             ENDIF
  388.          ENDIF
  389.          IPT2.ICOLOR(J)=IPT1.ICOLOR(J)
  390.  4    CONTINUE
  391.       RETURN
  392.       END
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  

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