Télécharger la plaquette Cast3M


fr en

Télécharger proob4.eso

Retour à la liste

Numérotation des lignes : 

proob4

  1. C PROOB4    SOURCE    PV        20/04/01    21:16:15     10569          
  2. C CE PROGRAMME PERMET LA PROJECTION CONIQUE D'UN MAILLAGE 3D SUR
  3. C UNE SUFACE CARACTERISEE PAR UN AUTRE MAILLAGE 3D.
  4.  
  5. C ENTREES/ MAI2  : UN MAILLAGE DE TYPE MELEME A PROJETER.
  6. C          MAI1  : UN MAILLAGE REPRESENTANT LA SURFACE SUR LAQUELLE MAI2
  7. C                  EST PROJETE.
  8. C          IP1   : UN POINT DEFINISSANT L'OEIL DE LA PROJECTION.
  9.  
  10. C DATE: SEPTEMBRE 1996.
  11.  
  12. C   **********************************************************************
  13.         Subroutine proob4(mai2,mai1,ip1)
  14.       IMPLICIT INTEGER(I-N)
  15.         IMPLICIT REAL*8(A-H,O-Z)
  16.  
  17.  
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC SMELEME
  21. -INC SMCOORD
  22. -INC CCREEL
  23. *-
  24.         SEGMENT ISEG1
  25.         REAL*8 XLIM(2,NBEL),YLIM(2,NBEL),ZLIM(2,NBEL)
  26.         ENDSEGMENT
  27.  
  28.         SEGMENT ISEG3
  29.         INTEGER NIZO(NZO+1)
  30.         ENDSEGMENT
  31.  
  32.         SEGMENT ISEG4
  33.         INTEGER NUMZO(NZO)
  34.         ENDSEGMENT
  35.  
  36.         SEGMENT ISEG5
  37.         INTEGER NNMEL(ILON),IDEJ(NZO)
  38.         ENDSEGMENT
  39.  
  40.         SEGMENT ISEG6
  41.         REAL*8 AM(4,4),AL(4),A(4),B(4),C(4)
  42.         REAL*8 XPU(2),P(3)
  43.         ENDSEGMENT
  44.  
  45.         SEGMENT MCOOR1
  46.         REAL*8  XCOOR1(XCOOR(/1))
  47.         ENDSEGMENT
  48. C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.
  49.  
  50.         SEGMENT MCOOR2
  51.         REAL*8  XCOOR2(XCOOR(/1))
  52.         ENDSEGMENT
  53. C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.
  54.  
  55.         SEGMENT MCOOR3
  56.         REAL*8  XCOOR3(XCOOR(/1))
  57.         ENDSEGMENT
  58. C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.
  59.  
  60.         SEGMENT MCOOR4
  61.         REAL*8  XCOOR4(XCOOR(/1))
  62.         ENDSEGMENT
  63. C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.
  64.  
  65.         SEGMENT MCOR
  66.         INTEGER ICOR(IONMAX+1)
  67.         ENDSEGMENT
  68.  
  69.         SEGMENT MCORRE
  70.         INTEGER ICORRE((nbpts))
  71.         ENDSEGMENT
  72.  
  73.         SEGMENT ICP1
  74.         INTEGER ICPR1(nbpts)
  75.         ENDSEGMENT
  76.  
  77.         SEGMENT ICP2
  78.         INTEGER  ICPR2(nbpts)
  79.         ENDSEGMENT
  80.  
  81.         SEGMENT MATRIC
  82.         REAL*8 TABMAT(3,3),V(3),W(3),Q(3),invma(3,3)
  83.         ENDSEGMENT
  84.  
  85.         SEGMENT CRIT
  86.         INTEGER ELECRI(NBELEM),TABCRIT(NBELEM+1)
  87.         ENDSEGMENT
  88.  
  89.  
  90.         segini matric
  91.         segact mcoord*mod
  92.         zchoix=100.
  93.         w(1)=xcoor((ip1-1)*(idim+1)+1)
  94.         w(2)=xcoor((ip1-1)*(idim+1)+2)
  95.         w(3)=xcoor((ip1-1)*(idim+1)+3)
  96.  
  97. C PASSAGE DU MAILLAGE MAI1 EN TRI3
  98.         call change(mai1,4)
  99.  
  100. C LISTAGE DES POINTS APPARTENANT A MAI1
  101.         segini icp1
  102.         meleme=mai1
  103.         segact meleme
  104.         nbelem=num(/2)
  105.         nbnn = num(/1)
  106.         Do 10 i=1,nbelem
  107.            Do 20 j=1,nbnn
  108.               ia=num(j,i)
  109.               icpr1(ia)=1
  110. 20          continue
  111. 10      continue
  112.  
  113. C LISTAGE DES POINTS APPARTENANT A MAI2
  114.         ipt1=mai2
  115.         segact ipt1
  116.         segini icp2
  117.         nbelem= ipt1.num(/2)
  118.         nbnn = ipt1.num(/1)
  119.         Do 30 i=1,nbelem
  120.            Do 40 j=1,nbnn
  121.               ia=ipt1.num(j,i)
  122.               If (icpr1(ia).eq.1) then
  123.               icpr2(ia)=2
  124.               else
  125.               icpr2(ia)=1
  126.               Endif
  127. 40          continue
  128. 30      continue
  129.  
  130. C POSITION DU CENTRE DE GRAVITE DE MAI2
  131.         sx=0.
  132.         sy=0.
  133.         sz=0.
  134.         isd=0
  135.         Do 2875 ia=1,nbpts
  136.            If (icpr2(ia).ne.0) then
  137.            sx=sx+xcoor((ia-1)*(idim+1)+1)
  138.            sy=sy+xcoor((ia-1)*(idim+1)+2)
  139.            sz=sz+xcoor((ia-1)*(idim+1)+3)
  140.            isd = isd+1
  141.            Endif
  142. 2875    continue
  143.  
  144. C CREATION DE LA MATRICE DE CHANGEMENT DE BASE
  145.         v(1)=(sx/isd)-w(1)
  146.         v(2)=(sy/isd)-w(2)
  147.         v(3)=(sz/isd)-w(3)
  148.         If ((abs(v(1))).lt.(1.E-3) .and. (abs(v(2))).lt.(1.E-3) .and.
  149.      $ (abs(v(3))).lt.(1.E-3)) then
  150.         v(1)=v(1)+0.1
  151.         v(2)=v(2)+0.1
  152.         v(3)=v(3)+0.1
  153.         Endif
  154.         segini mcoor3
  155.         segini mcoor1
  156.         segini mcoor4
  157.         segini mcoor2
  158.         r1=sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
  159.         v(1) = V(1) / r1
  160.         v(2) = v(2) / r1
  161.         v(3) = v(3) / r1
  162.         r2=sqrt(v(1)*v(1)+v(2)*v(2))
  163.         if ( r2 . gt. 1.e-9) then
  164.         tabmat(1,1)=-v(2)/r2
  165.         tabmat(1,2)=v(1)/r2
  166.         tabmat(1,3)=0.
  167.         else
  168.         tabmat(1,1)=1.
  169.         tabmat(1,2)=0.
  170.         tabmat(1,3)=0.
  171.         endif
  172.         tabmat(2,1)=-v(3)*tabmat(1,2)
  173.         tabmat(2,2)=v(3)*tabmat(1,1)
  174.         tabmat(2,3)=v(1)*tabmat(1,2)-v(2)*tabmat(1,1)
  175.         tabmat(3,1)=v(1)
  176.         tabmat(3,2)=v(2)
  177.         tabmat(3,3)=v(3)
  178. * pas besoin de diviser par le determinant car egal à 1.
  179.         INVMA(1,1)= TABMAT(2,2)*TABMAT(3,3)-TABMAT(3,2)*TABMAT(2,3)
  180.         INVMA(1,2)=-TABMAT(1,2)*TABMAT(3,3)+TABMAT(3,2)*TABMAT(1,3)
  181.         INVMA(1,3)= TABMAT(1,2)*TABMAT(2,3)-TABMAT(2,2)*TABMAT(1,3)
  182.         INVMA(2,1)=-TABMAT(2,1)*TABMAT(3,3)+TABMAT(3,1)*TABMAT(2,3)
  183.         INVMA(2,2)= TABMAT(1,1)*TABMAT(3,3)-TABMAT(3,1)*TABMAT(1,3)
  184.         INVMA(2,3)=-TABMAT(1,1)*TABMAT(2,3)+TABMAT(2,1)*TABMAT(1,3)
  185.         INVMA(3,1)= TABMAT(2,1)*TABMAT(3,2)-TABMAT(2,2)*TABMAT(3,1)
  186.         INVMA(3,2)=-TABMAT(1,1)*TABMAT(3,2)+TABMAT(3,1)*TABMAT(1,2)
  187.         INVMA(3,3)= TABMAT(1,1)*TABMAT(2,2)-TABMAT(2,1)*TABMAT(1,2)
  188. C   CALCUL DES COORDONNEES DANS LE NOUVEAU REPERE POUR MELEME
  189.         nbelem=num(/2)
  190.         nbnn=num(/1)
  191.         segini crit
  192.         Do 31 i1=1,nbelem
  193.             Do 38 k=1,nbnn
  194.                  ia=num(k,i1)
  195.                  xia1=xcoor((ia-1)*(idim+1)+1)-w(1)
  196.                  xia2=xcoor((ia-1)*(idim+1)+2)-w(2)
  197.                  xia3=xcoor((ia-1)*(idim+1)+3)-w(3)
  198.                  Do 32 j=1,3
  199.                     xcoor3((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+
  200.      $(xia2*tabmat(j,2))+(xia3*tabmat(j,3)))
  201. 32               continue
  202. 38         continue
  203. 31      continue
  204.  
  205. C   CALCUL DE LA TAILLE LA PLUS GRANDE DES ELEMENTS SUIVANT OZ
  206.         zlong=-xgrand
  207.         Do 1534 i1=1,nbelem
  208.             zmin=xgrand
  209.             zmax=-xgrand
  210.             Do 1233 i2=1,nbnn
  211.               ib=num(i2,i1)
  212.               ia=(ib-1)*4
  213.               If (xcoor3(ia+3).lt.zmin)  zmin=xcoor3(ia+3)
  214.               If (xcoor3(ia+3).gt.zmax)  zmax=xcoor3(ia+3)
  215. 1233        continue
  216.             zlong=max(zlong,zmax-zmin)
  217. 1534    continue
  218.         If (zlong.lt.(5.E-2))  zlong=(5.E-2)
  219.  
  220. C DETERMINATION DES ELEMENTS TOUCHANT LA TRANCHE CRITIQUE
  221.         Do 317 i1=1,nbelem
  222.             ime=0
  223.             Do 387 k=1,nbnn
  224.                  ia=num(k,i1)
  225.                   q(1)=(xcoor3((ia-1)*(idim+1)+1))
  226.                   q(2)=(xcoor3((ia-1)*(idim+1)+2))
  227.                   q(3)=(xcoor3((ia-1)*(idim+1)+3))
  228.                 If (abs(xcoor3((ia-1)*(idim+1)+3)).gt.(0.55*zlong)) then
  229.                 xcoor1((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1)
  230.                 xcoor1((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2)
  231.                 xcoor1((ia-1)*(idim+1)+3)=zchoix
  232.                 else
  233.                 ime=ime+1
  234.                 Endif
  235. 387         continue
  236.            If (ime.ne.0) then
  237.            tabcrit(1)=tabcrit(1)+1
  238.            tabcrit(tabcrit(1)+1)=i1
  239.            elecri(i1)=1
  240.            Endif
  241. 317      continue
  242.  
  243. C   CHANGEMENT DE REPERE POUR IPT1
  244.         Do 33 ia=1,nbpts
  245.               If (icpr2(ia).eq.1)  then
  246.                  xia1=xcoor((ia-1)*(idim+1)+1)-w(1)
  247.                  xia2=xcoor((ia-1)*(idim+1)+2)-w(2)
  248.                  xia3=xcoor((ia-1)*(idim+1)+3)-w(3)
  249.                  Do 34 j=1,3
  250.                     xcoor4((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+
  251.      $(xia2*tabmat(j,2))+(xia3*tabmat(j,3)))
  252. 34               continue
  253.             If (abs(xcoor4((ia-1)*(idim+1)+3)).lt.(0.55*zlong)) then
  254.                  icpr2(ia)=3
  255.                  else
  256.                  q(1)=(xcoor4((ia-1)*(idim+1)+1))
  257.                  q(2)=(xcoor4((ia-1)*(idim+1)+2))
  258.                  q(3)=(xcoor4((ia-1)*(idim+1)+3))
  259.                xcoor2((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1)
  260.                xcoor2((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2)
  261.                xcoor2((ia-1)*(idim+1)+3)=zchoix
  262.                  Endif
  263.                 Endif
  264.               If (icpr2(ia).eq.2)  then
  265.                  xia1=xcoor((ia-1)*(idim+1)+1)-w(1)
  266.                  xia2=xcoor((ia-1)*(idim+1)+2)-w(2)
  267.                  xia3=xcoor((ia-1)*(idim+1)+3)-w(3)
  268.                  Do 35 j=1,3
  269.                     xcoor2((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+
  270.      $(xia2*tabmat(j,2))+(xia3*tabmat(j,3)))
  271. 35               continue
  272.                  q(1)=(xcoor4((ia-1)*(idim+1)+1))
  273.                  q(2)=(xcoor4((ia-1)*(idim+1)+2))
  274.                  q(3)=(xcoor4((ia-1)*(idim+1)+3))
  275.                xcoor2((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1)
  276.                xcoor2((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2)
  277.                  xcoor2((ia-1)*(idim+1)+3)=zchoix
  278.               Endif
  279. 33      continue
  280.  
  281. C******************** DEBUT DU ZONAGE DE MAI1 **********************
  282.  
  283. C ON CALCULE LA TAILLE MAXI D'UN ELEMENT DANS TOUTES LES DIRECTIONS
  284. C AFIN DE CREER UN ZONAGE DE L'ESPACE. EN MEME TEMPS ON CALCULE
  285. C LA DIMENSION HORS TOUT DU MAILLAGE
  286. C
  287.       IDIM1=4
  288.       NBEL = NUM(/2)
  289.       NBNN=NUM(/1)
  290.       SEGINI ISEG1
  291.       ILOC=0
  292.       XZO=0.D0
  293.       YZO=0.D0
  294.       ZZO=0.D0
  295.       XZA=XGRAND
  296.       YZA=XGRAND
  297.       ZZA=XGRAND
  298.       XTOMI=XGRAND
  299.       XTOMA=-XGRAND
  300.       YTOMI=XGRAND
  301.       YTOMA=-XGRAND
  302.       ZTOMI=XGRAND
  303.       ZTOMA=-XGRAND
  304.       DO 1 I1=1,NBEL
  305.          If (elecri(i1).eq.0) then
  306.          XMI=XGRAND
  307.          YMI=XGRAND
  308.          ZMI=XGRAND
  309.          YMA=-XGRAND
  310.          XMA=-XGRAND
  311.          ZMA=-XGRAND
  312.          DO 2 I2 = 1,NBNN
  313.             IB=NUM(I2,I1)
  314.             IA=(IB-1)*IDIM1
  315.             IF(XCOOR1(IA+1).LT.XMI) XMI=XCOOR1(IA+1)
  316.             IF(XCOOR1(IA+1).GT.XMA) XMA=XCOOR1(IA+1)
  317.             IF(XCOOR1(IA+2).LT.YMI) YMI=XCOOR1(IA+2)
  318.             IF(XCOOR1(IA+2).GT.YMA) YMA=XCOOR1(IA+2)
  319.     2    CONTINUE
  320.          XLIM(1,I1)=XMI
  321.          XLIM(2,I1)=XMA
  322.          YLIM(1,I1)=YMI
  323.          YLIM(2,I1)=YMA
  324.          XZO=MAX (XZO,XMA-XMI)
  325.          YZO=MAX (YZO,YMA-YMI)
  326.          XZA=MIN(XZA,XMA-XMI)
  327.          YZA=MIN(YZA,YMA-YMI)
  328.          IF(XMI.LT.XTOMI) XTOMI=XMI
  329.          IF(XMA.GT.XTOMA) XTOMA=XMA
  330.          IF(YMI.LT.YTOMI) YTOMI=YMI
  331.          IF(YMA.GT.YTOMA) YTOMA=YMA
  332.        Endif
  333.     1 CONTINUE
  334.       XPR=MIN(XZO*1.D-2,(XTOMA-XTOMI)/2.D+4)
  335.       YPR=MIN(YZO*1.D-2,(YTOMA-YTOMI)/2.D+4)
  336.       XZA=XZA*0.97
  337.       YZA=YZA*0.97
  338.       XTOMI= XTOMI - (XTOMA-XTOMI)/1.D+4
  339.       XTOMA= XTOMA + (XTOMA-XTOMI)/1.D+4
  340.       YTOMI= YTOMI - (YTOMA-YTOMI)/1.D+4
  341.       YTOMA= YTOMA + (YTOMA-YTOMI)/1.D+4
  342.       XZA=MIN ( XZA, XTOMA-XTOMI)
  343.       YZA=MIN ( YZA, YTOMA-YTOMI)
  344.       XZA = max ( xza , (XTOMA-XTOMI) / 50)
  345.       YZA = max ( yza , (YTOMA-YTOMI) / 50)
  346.       NXZO=(XTOMA-XTOMI)/XZA + 1
  347.       NYZO=(YTOMA-YTOMI)/YZA + 1
  348.       XZO=XZA
  349.       YZO=YZA
  350.       NZZO=1
  351.        NXDEP=MIN(NXZO,10)
  352.        NYDEP=MIN(NYZO,10)
  353.          IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN
  354.           XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/90
  355.           NXZO=MAX(INT(NXZO/XY),NXDEP)
  356.           NYZO=MAX(INT(NYZO/XY),NYDEP)
  357.           IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN
  358.             XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/60
  359.             NXZO=MAX(INT(NXZO/XY),NXDEP)
  360.             NYZO=MAX(INT(NYZO/XY),NYDEP)
  361.           ENDIF
  362.           XZO=(XTOMA-XTOMI)/NXZO
  363.           YZO=(YTOMA-YTOMI)/NYZO
  364.           NXZO=(XTOMA-XTOMI)/XZO +1
  365.           NYZO=(YTOMA-YTOMI)/YZO +1
  366.          ENDIF
  367.  
  368. C
  369. C  ON VEUT CONSTRUIRE LA LISTE DES ELEMENTS TOUCHANT UNE ZONE
  370. C  POUR CELA ON COMMENCE PAR COMPTER COMBIEN D'ELEMENT TOUCHENT
  371. C  CHAQUE ZONE ET EN MEME TEMPS ON STOCKE LES ZONES TOUCHEES
  372. C  PAR CHAQUE ELEMENT ET LEUR NOMBRE
  373. C
  374.  
  375.       NZO=NXZO*NYZO
  376.       IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NZO NXZO NYZO NZZO  ''
  377.      $,4I7) ') NZO,NXZO,NYZO,NZZO
  378.       NXYZO=NXZO*NYZO
  379.       SEGINI ISEG3
  380.       SEGINI ISEG4
  381.       DO 3 I1=1,NBEL
  382.           If (elecri(i1).eq.0) then
  383.           NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
  384.           NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
  385.           NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
  386.           NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
  387.                DO 201 L1=NIZ1Y,NIZ2Y
  388.                DO 201 L2=NIZ1X,NIZ2X
  389.                    NIZA = L2 + ( L1-1) * NXZO
  390.                    NUMZO(NIZA) = NUMZO(NIZA) +1
  391.   201          CONTINUE
  392.       Endif
  393.     3 CONTINUE
  394. C
  395. C  CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES
  396. C  ELEMENTS CONCERNEES PAR UNE ZONE
  397. C
  398.       ILON=0
  399.       NIZO(1)=1
  400.       DO 202 L1=1,NZO
  401.       NIZO(L1+1)=NIZO(L1)+NUMZO(L1)
  402.       ILON=ILON+ NUMZO(L1)
  403.  202  CONTINUE
  404.  110  FORMAT(16I5)
  405.       SEGINI ISEG5
  406.       DO 5 I1=1,NBEL
  407.           If (elecri(i1).eq.0) then
  408.           NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
  409.           NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
  410.           NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
  411.           NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
  412.                DO 203 L1=NIZ1Y,NIZ2Y
  413.                DO 203 L2=NIZ1X,NIZ2X
  414.                    NIZA = L2 + ( L1-1) * NXZO
  415.                    IAD=NIZO(NIZA)+IDEJ(NIZA)
  416.                    NNMEL(IAD)=I1
  417.                    IDEJ(NIZA)=IDEJ(NIZA)+1
  418.   203          CONTINUE
  419.        Endif
  420.     5  CONTINUE
  421. C *********************** FIN DU ZONAGE DE MAI1 **************************
  422.         prec = 1.E-5
  423.         prec2 = -prec
  424.         prec3 = -5.E-3
  425. C ******RECHERCHE DU NOMBRE D'ELEMENT MAX CONTENU DANS UNE ZONE***********
  426.         la=0
  427.         ls=0
  428.         Do 420 i=1,nzo
  429.           ls=ls+numzo(i)
  430.           If (numzo(i).gt.la) then
  431.           la=numzo(i)
  432.           Endif
  433. 420     continue
  434.         ionmax=la
  435.  
  436. C ****************** TRAITEMENT SUR LES POINTS DE MAI2 *******************
  437.       nbelem=num(/2)
  438.       nbnn=num(/1)
  439.       nbpts2=nbpts
  440.       segini mcor,mcorre,iseg6
  441.       Do 405 ia=1,nbpts2
  442.          If (icpr2(ia).eq.2) then
  443.            icorre(ia)=ia
  444.          Endif
  445.          If (icpr2(ia).eq.1) then
  446.               xpu(1)=xcoor2((ia-1)*(idim+1)+1)
  447.               xpu(2)=xcoor2((ia-1)*(idim+1)+2)
  448.                If (tabcrit(1).eq.0) then
  449.                    if (((xpu(1)-xtomi).lt.prec.) .or.
  450.      $              ((xpu(2)-ytomi).lt.prec.) .or.
  451.      $              ((xpu(1)-xtoma).gt.prec2.) .or.
  452.      $              ((xpu(2)-ytoma).gt.prec2.)) then
  453.                       call erreur(21)
  454.                       return
  455.                    endif
  456.              Endif
  457. C    RECHERCHE DU NUMERO DE LA ZONE CORRESPONDANTE
  458.              k2=int(((xcoor2((ia-1)*(idim+1)+1))-xtomi-xpr)/xzo)+1
  459.              k1=int(((xcoor2((ia-1)*(idim+1)+2))-ytomi-ypr)/yzo)+1
  460.              If ((k2.gt.0.) .and. (k2.lt.(nxzo+1)) .and.
  461.      $       (k1.gt.0.) .and. (k1.lt.(nyzo+1))) then
  462.                   niza=k2+(k1-1)*nxzo
  463. C    INVENTAIRE DES ELEMENTS CONCERNES PAR LE POINT IA
  464.                  iad=nizo(niza)
  465.                  Do 406 i=1,numzo(niza)
  466.                   i1=nnmel(iad)
  467. C    CALCUL DES COORDONNEES BARYCENTRIQUES DE IA DANS L'ELEMENT I1
  468.                   j1=num(1,i1)
  469.                   j2=num(2,i1)
  470.                   j3=num(3,i1)
  471.                   j1idim=(j1-1)*(idim+1)
  472.                   j2idim=(j2-1)*(idim+1)
  473.                   j3idim=(j3-1)*(idim+1)
  474.                   Do 408 l=1,2
  475.                     p(l+1)=xpu(l)
  476.                      am(l+1,1)=xcoor1(j1idim+l)
  477.                     am(l+1,2)=xcoor1(j2idim+l)
  478.                     am(l+1,3)=xcoor1(j3idim+l)
  479. 408              continue
  480.                  x1=am(2,1)
  481.                  x2=am(2,2)
  482.                  x3=am(2,3)
  483.                  y1=am(3,1)
  484.                  y2=am(3,2)
  485.                  y3=am(3,3)
  486.                  x=p(2)
  487.                  y=p(3)
  488.                  detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  489.                  a(1)=x2*y3-y2*x3
  490.                  a(2)=x3*y1-x1*y3
  491.                  a(3)=x1*y2-x2*y1
  492.                  b(1)=y2-y3
  493.                  b(2)=y3-y1
  494.                  b(3)=y1-y2
  495.                  c(1)=x3-x2
  496.                  c(2)=x1-x3
  497.                  c(3)=x2-x1
  498.                  Do 409 k=1,nbnn
  499.                    al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  500. 409              continue
  501. C    TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  502.                 If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  503.      $          al(3).gt.prec3) then
  504. C    LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  505.                     icor(1)=icor(1)+1
  506.                     icor(icor(1)+1)=i1
  507.                Endif
  508.                iad=iad+1
  509.  
  510. 406            continue
  511.              Endif
  512. C   PASSAGE TRANCHE CRITIQUE POUR UN POINT DU ZONAGE
  513.              If (k2.lt.2 .or. k2.gt.(nxzo-1)) then
  514.                Do 1462 iz=1,tabcrit(1)
  515.                i1=tabcrit(iz+1)
  516.                m1=num(1,i1)
  517.                m2=num(2,i1)
  518.                m3=num(3,i1)
  519. C     CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2)
  520.                xm1=xcoor3((m1-1)*(idim+1)+1)
  521.                ym1=xcoor3((m1-1)*(idim+1)+2)
  522.                zm1=xcoor3((m1-1)*(idim+1)+3)
  523.                xm2=xcoor3((m2-1)*(idim+1)+1)
  524.                ym2=xcoor3((m2-1)*(idim+1)+2)
  525.                zm2=xcoor3((m2-1)*(idim+1)+3)
  526.                xm3=xcoor3((m3-1)*(idim+1)+1)
  527.                ym3=xcoor3((m3-1)*(idim+1)+2)
  528.                zm3=xcoor3((m3-1)*(idim+1)+3)
  529.                q(1)=(xcoor4((ia-1)*(idim+1)+1))
  530.                q(2)=(xcoor4((ia-1)*(idim+1)+2))
  531.                q(3)=(xcoor4((ia-1)*(idim+1)+3))
  532. C            Calcul du zib
  533.                o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  534.                o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  535.                o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  536.                o4=-(o1*xm1+o2*ym1+o3*zm1)
  537.                zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  538.                xnew=q(1)*zk
  539.                ynew=q(2)*zk
  540.                znew=q(3)*zk
  541. C   PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1
  542.                r1=sqrt(o1*o1+o2*o2+o3*o3)
  543.                if( r1 . EQ. 0. ) then
  544.                  call erreur ( 21 )
  545.               endif
  546.               o1 = o1/r1
  547.               o2 = o2/r1
  548.               o3 = o3/r1
  549.               r2=sqrt(o1*o1+o2*o2)
  550.               if ( r2 . gt. 1.e-5) then
  551.                 tabmat(1,1)=-o2/r2
  552.                 tabmat(1,2)=o1/r2
  553.                 tabmat(1,3)=0.
  554.               else
  555.                 tabmat(1,1)=1.
  556.                 tabmat(1,2)=0.
  557.                 tabmat(1,3)=0.
  558.               endif
  559.             tabmat(2,1)=-o3*tabmat(1,2)
  560.             tabmat(2,2)=o3*tabmat(1,1)
  561.             tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1)
  562.             tabmat(3,1)=o1
  563.             tabmat(3,2)=o2
  564.             tabmat(3,3)=o3
  565.         x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3)))
  566.         y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3)))
  567.         x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3)))
  568.         y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3)))
  569.         x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3)))
  570.         y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3)))
  571.         x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3)))
  572.         y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3)))
  573. C    CALCUL DES COORDONNEES BARYCENTRIQUES
  574.             detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  575.             a(1)=x2*y3-y2*x3
  576.             a(2)=x3*y1-x1*y3
  577.             a(3)=x1*y2-x2*y1
  578.             b(1)=y2-y3
  579.             b(2)=y3-y1
  580.             b(3)=y1-y2
  581.             c(1)=x3-x2
  582.             c(2)=x1-x3
  583.             c(3)=x2-x1
  584.             Do 798 k=1,nbnn
  585.                al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  586. 798         continue
  587. C    TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  588.            If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  589.      $     al(3).gt.prec3) then
  590. C    LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  591.                 icor(1)=icor(1)+1
  592.                 icor(icor(1)+1)=i1
  593.            Endif
  594. 1462       continue
  595.          Endif
  596.  
  597.          If (icor(1).eq.0) then
  598.            interr(1) = ia
  599.            call erreur(782)
  600.            return
  601.          Endif
  602.  
  603.       Endif
  604. C   FIN DU IF ICPR2(IA).EQ.1
  605.  
  606.          If (icpr2(ia).eq.3) then
  607.  
  608. C   PASSAGE SUR LES ELEMENTS DE LA TRANCHE CRITIQUE
  609.              Do 1463 iz=1,tabcrit(1)
  610.              i1=tabcrit(iz+1)
  611.              m1=num(1,i1)
  612.              m2=num(2,i1)
  613.              m3=num(3,i1)
  614. C     CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2)
  615.              xm1=xcoor3((m1-1)*(idim+1)+1)
  616.              ym1=xcoor3((m1-1)*(idim+1)+2)
  617.              zm1=xcoor3((m1-1)*(idim+1)+3)
  618.              xm2=xcoor3((m2-1)*(idim+1)+1)
  619.              ym2=xcoor3((m2-1)*(idim+1)+2)
  620.              zm2=xcoor3((m2-1)*(idim+1)+3)
  621.              xm3=xcoor3((m3-1)*(idim+1)+1)
  622.              ym3=xcoor3((m3-1)*(idim+1)+2)
  623.              zm3=xcoor3((m3-1)*(idim+1)+3)
  624.              q(1)=(xcoor4((ia-1)*(idim+1)+1))
  625.              q(2)=(xcoor4((ia-1)*(idim+1)+2))
  626.              q(3)=(xcoor4((ia-1)*(idim+1)+3))
  627. C            Calcul du zib
  628.              o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  629.              o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  630.              o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  631.              o4=-(o1*xm1+o2*ym1+o3*zm1)
  632.              zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  633.              xnew=q(1)*zk
  634.              ynew=q(2)*zk
  635.              znew=q(3)*zk
  636. C   PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1
  637.         r1=sqrt(o1*o1+o2*o2+o3*o3)
  638.         if( r1 . EQ. 0. ) then
  639.           call erreur ( 21 )
  640.         endif
  641.         o1 = o1/r1
  642.         o2 = o2/r1
  643.         o3 = o3/r1
  644.         r2=sqrt(o1*o1+o2*o2)
  645.         if ( r2 . gt. 1.e-5) then
  646.         tabmat(1,1)=-o2/r2
  647.         tabmat(1,2)=o1/r2
  648.         tabmat(1,3)=0.
  649.         else
  650.         tabmat(1,1)=1.
  651.         tabmat(1,2)=0.
  652.         tabmat(1,3)=0.
  653.         endif
  654.         tabmat(2,1)=-o3*tabmat(1,2)
  655.         tabmat(2,2)=o3*tabmat(1,1)
  656.         tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1)
  657.         tabmat(3,1)=o1
  658.         tabmat(3,2)=o2
  659.         tabmat(3,3)=o3
  660.         x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3)))
  661.         y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3)))
  662.         x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3)))
  663.         y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3)))
  664.         x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3)))
  665.         y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3)))
  666.         x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3)))
  667.         y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3)))
  668. C    CALCUL DES COORDONNEES BARYCENTRIQUES
  669.         detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  670.                a(1)=x2*y3-y2*x3
  671.                a(2)=x3*y1-x1*y3
  672.                a(3)=x1*y2-x2*y1
  673.                b(1)=y2-y3
  674.                b(2)=y3-y1
  675.                b(3)=y1-y2
  676.                c(1)=x3-x2
  677.                c(2)=x1-x3
  678.                c(3)=x2-x1
  679.                Do 799 k=1,nbnn
  680.                   al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  681. 799            continue
  682. C    TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  683.              If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  684.      $al(3).gt.prec3) then
  685. C    LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  686.              icor(1)=icor(1)+1
  687.              icor(icor(1)+1)=i1
  688.              Endif
  689. 1463       continue
  690.          If (icor(1).eq.0) then
  691.            interr(1)=ia
  692.            call erreur(782)
  693.             return
  694.          Endif
  695.          Endif
  696. C    FIN DE L'INVENTAIRE
  697.  
  698. C    CALCUL DES PROJETES DE MAI2 SUR MAI1
  699.          If  ((icpr2(ia).eq.1) .or. (icpr2(ia).eq.3)) then
  700.  
  701.               dmin=1.e30
  702.               nbpts=nbpts+1
  703.               segadj mcoord
  704.               n=0
  705.               Do 500 m=1,icor(1)
  706.                 m1=num(1,icor(m+1))
  707.                 m2=num(2,icor(m+1))
  708.                 m3=num(3,icor(m+1))
  709. C               Calcul des coordonnees du projete
  710.                 xm1=xcoor3((m1-1)*(idim+1)+1)
  711.                 ym1=xcoor3((m1-1)*(idim+1)+2)
  712.                 zm1=xcoor3((m1-1)*(idim+1)+3)
  713.                 xm2=xcoor3((m2-1)*(idim+1)+1)
  714.                 ym2=xcoor3((m2-1)*(idim+1)+2)
  715.                 zm2=xcoor3((m2-1)*(idim+1)+3)
  716.                 xm3=xcoor3((m3-1)*(idim+1)+1)
  717.                 ym3=xcoor3((m3-1)*(idim+1)+2)
  718.                 zm3=xcoor3((m3-1)*(idim+1)+3)
  719.                 q(1)=(xcoor4((ia-1)*(idim+1)+1))
  720.                 q(2)=(xcoor4((ia-1)*(idim+1)+2))
  721.                 q(3)=(xcoor4((ia-1)*(idim+1)+3))
  722. C               Calcul du zib
  723.                 o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  724.                 o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  725.                 o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  726.                 o4=-(o1*xm1+o2*ym1+o3*zm1)
  727.                 zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  728.                 xnew=q(1)*zk
  729.                 ynew=q(2)*zk
  730.                 znew=q(3)*zk
  731.                 delta1=sqrt(xnew*xnew+ynew*ynew+znew*znew)
  732.                 delta2=sqrt(q(1)*q(1)+q(2)* q(2)+q(3)*q(3))
  733.                 pds=(xnew*q(1))+(ynew*q(2))+(znew*q(3))
  734.                 If (pds.gt.0.) then
  735.                  If ((delta1-delta2).gt.0.) then
  736.                   If ((delta1-delta2).lt.dmin) then
  737.                        icorre(ia) = nbpts
  738.                        dmin=(delta1-delta2)
  739.                        xcoor((nbpts-1)*(idim+1)+1)=xnew
  740.                        xcoor((nbpts-1)*(idim+1)+2)=ynew
  741.                        xcoor((nbpts-1)*(idim+1)+3)=znew
  742.                   Endif
  743.                   n=n+1
  744.                 Endif
  745.                 Endif
  746. 500           continue
  747.            If (n.eq.0) then
  748.            call erreur(21)
  749.            return
  750.            Endif
  751.  
  752. C *********************RETOUR A L'ANCIENNE BASE*********************
  753.         ib=nbpts
  754.         xa1=xcoor((ib-1)*(idim+1)+1)
  755.         xa2=xcoor((ib-1)*(idim+1)+2)
  756.         xa3=xcoor((ib-1)*(idim+1)+3)
  757.         do 3200 j=1,3
  758.         xcoor((ib-1)*(idim+1)+j)=(xa1*invma(j,1))+
  759.      $(xa2*invma(j,2))+(xa3*invma(j,3))+w(J)
  760. 3200               continue
  761.         icor(1)=0
  762.         n=0
  763.         Endif
  764. 405     continue
  765.         segsup iseg1,iseg3,iseg4,iseg5,iseg6,mcoor1,mcoor2
  766.         segsup icp1,icp2,mcoor3,mcoor4
  767. C   CREATION DU MAILLAGE PROJETE MAI3
  768.         segini,ipt3=ipt1
  769.         ipt4 = IPT3
  770.         do 2456 IU = 1, max( 1,IPT3.LISOUS(/1))
  771.         If( IPT3.LISOUS(/1).NE.0) THEN
  772.           IPT5=IPT3.LISOUS(IU)
  773.           SEGINI,IPT4=IPT5
  774.           NBREF = 0
  775.           NBSOUS=0
  776.           SEGADJ IPT4
  777.         ENDIF
  778.         nbel1=ipt4.num(/2)
  779.         nbnn1=ipt4.num(/1)
  780.         Do 600 i=1,nbel1
  781.             Do 601 j=1,nbnn1
  782.             ia=ipt4.num(j,i)
  783.             ib=icorre(ia)
  784.             ipt4.num(j,i)=ib
  785. 601         continue
  786. 600     continue
  787. 2456    continue
  788.         do 2457 IU = 1, max( 0,IPT3.LISREF(/1))
  789.         If( IPT3.LISREF(/1).NE.0) THEN
  790.          IPT5 = IPT3.LISREF(IU)
  791.          SEGINI,IPT4=IPT5
  792.         ENDIF
  793.         nbel1=ipt4.num(/2)
  794.         nbnn1=ipt4.num(/1)
  795.         Do 602 i=1,nbel1
  796.             Do 603 j=1,nbnn1
  797.             ia=ipt4.num(j,i)
  798.             ib=icorre(ia)
  799.             ipt4.num(j,i)=ib
  800. 603         continue
  801. 602     continue
  802.         SEGDES IPT4
  803. 2457    continue
  804.         call ecrobj('MAILLAGE',ipt3)
  805.         segdes ipt3
  806.         segsup mcor,mcorre,matric
  807.         Return
  808.         End
  809.  
  810.  
  811.  
  812.  
  813.  
  814.  
  815.  
  816.  
  817.  
  818.  
  819.  
  820.  
  821.  

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