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    CB215821  26/08/24    21:18:01     12622          
  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 MCOORA
  46.         REAL*8  XCOORA(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 mcoora
  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.                 xcoora((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1)
  230.                 xcoora((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2)
  231.                 xcoora((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(XCOORA(IA+1).LT.XMI) XMI=XCOORA(IA+1)
  316.             IF(XCOORA(IA+1).GT.XMA) XMA=XCOORA(IA+1)
  317.             IF(XCOORA(IA+2).LT.YMI) YMI=XCOORA(IA+2)
  318.             IF(XCOORA(IA+2).GT.YMA) YMA=XCOORA(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 3201 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.  3201          CONTINUE
  393.       Endif
  394.     3 CONTINUE
  395. C
  396. C  CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES
  397. C  ELEMENTS CONCERNEES PAR UNE ZONE
  398. C
  399.       ILON=0
  400.       NIZO(1)=1
  401.       DO 202 L1=1,NZO
  402.       NIZO(L1+1)=NIZO(L1)+NUMZO(L1)
  403.       ILON=ILON+ NUMZO(L1)
  404.  202  CONTINUE
  405.  110  FORMAT(16I5)
  406.       SEGINI ISEG5
  407.       DO 5 I1=1,NBEL
  408.           If (elecri(i1).eq.0) then
  409.           NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
  410.           NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
  411.           NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
  412.           NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
  413.                DO 3202 L1=NIZ1Y,NIZ2Y
  414.                DO 203 L2=NIZ1X,NIZ2X
  415.                    NIZA = L2 + ( L1-1) * NXZO
  416.                    IAD=NIZO(NIZA)+IDEJ(NIZA)
  417.                    NNMEL(IAD)=I1
  418.                    IDEJ(NIZA)=IDEJ(NIZA)+1
  419.   203          CONTINUE
  420.  3202          CONTINUE
  421.        Endif
  422.     5  CONTINUE
  423. C *********************** FIN DU ZONAGE DE MAI1 **************************
  424.         prec = 1.E-5
  425.         prec2 = -prec
  426.         prec3 = -5.E-3
  427. C ******RECHERCHE DU NOMBRE D'ELEMENT MAX CONTENU DANS UNE ZONE***********
  428.         la=0
  429.         ls=0
  430.         Do 420 i=1,nzo
  431.           ls=ls+numzo(i)
  432.           If (numzo(i).gt.la) then
  433.           la=numzo(i)
  434.           Endif
  435. 420     continue
  436.         ionmax=la
  437.  
  438. C ****************** TRAITEMENT SUR LES POINTS DE MAI2 *******************
  439.       nbelem=num(/2)
  440.       nbnn=num(/1)
  441.       nbpts2=nbpts
  442.       segini mcor,mcorre,iseg6
  443.       Do 405 ia=1,nbpts2
  444.          If (icpr2(ia).eq.2) then
  445.            icorre(ia)=ia
  446.          Endif
  447.          If (icpr2(ia).eq.1) then
  448.               xpu(1)=xcoor2((ia-1)*(idim+1)+1)
  449.               xpu(2)=xcoor2((ia-1)*(idim+1)+2)
  450.                If (tabcrit(1).eq.0) then
  451.                    if (((xpu(1)-xtomi).lt.prec.) .or.
  452.      $              ((xpu(2)-ytomi).lt.prec.) .or.
  453.      $              ((xpu(1)-xtoma).gt.prec2.) .or.
  454.      $              ((xpu(2)-ytoma).gt.prec2.)) then
  455.                       call erreur(21)
  456.                       return
  457.                    endif
  458.              Endif
  459. C    RECHERCHE DU NUMERO DE LA ZONE CORRESPONDANTE
  460.              k2=int(((xcoor2((ia-1)*(idim+1)+1))-xtomi-xpr)/xzo)+1
  461.              k1=int(((xcoor2((ia-1)*(idim+1)+2))-ytomi-ypr)/yzo)+1
  462.              If ((k2.gt.0.) .and. (k2.lt.(nxzo+1)) .and.
  463.      $       (k1.gt.0.) .and. (k1.lt.(nyzo+1))) then
  464.                   niza=k2+(k1-1)*nxzo
  465. C    INVENTAIRE DES ELEMENTS CONCERNES PAR LE POINT IA
  466.                  iad=nizo(niza)
  467.                  Do 406 i=1,numzo(niza)
  468.                   i1=nnmel(iad)
  469. C    CALCUL DES COORDONNEES BARYCENTRIQUES DE IA DANS L'ELEMENT I1
  470.                   j1=num(1,i1)
  471.                   j2=num(2,i1)
  472.                   j3=num(3,i1)
  473.                   j1idim=(j1-1)*(idim+1)
  474.                   j2idim=(j2-1)*(idim+1)
  475.                   j3idim=(j3-1)*(idim+1)
  476.                   Do 408 l=1,2
  477.                     p(l+1)=xpu(l)
  478.                      am(l+1,1)=xcoora(j1idim+l)
  479.                     am(l+1,2)=xcoora(j2idim+l)
  480.                     am(l+1,3)=xcoora(j3idim+l)
  481. 408              continue
  482.                  x1=am(2,1)
  483.                  x2=am(2,2)
  484.                  x3=am(2,3)
  485.                  y1=am(3,1)
  486.                  y2=am(3,2)
  487.                  y3=am(3,3)
  488.                  x=p(2)
  489.                  y=p(3)
  490.                  detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  491.                  a(1)=x2*y3-y2*x3
  492.                  a(2)=x3*y1-x1*y3
  493.                  a(3)=x1*y2-x2*y1
  494.                  b(1)=y2-y3
  495.                  b(2)=y3-y1
  496.                  b(3)=y1-y2
  497.                  c(1)=x3-x2
  498.                  c(2)=x1-x3
  499.                  c(3)=x2-x1
  500.                  Do 409 k=1,nbnn
  501.                    al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  502. 409              continue
  503. C    TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  504.                 If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  505.      $          al(3).gt.prec3) then
  506. C    LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  507.                     icor(1)=icor(1)+1
  508.                     icor(icor(1)+1)=i1
  509.                Endif
  510.                iad=iad+1
  511.  
  512. 406            continue
  513.              Endif
  514. C   PASSAGE TRANCHE CRITIQUE POUR UN POINT DU ZONAGE
  515.              If (k2.lt.2 .or. k2.gt.(nxzo-1)) then
  516.                Do 1462 iz=1,tabcrit(1)
  517.                i1=tabcrit(iz+1)
  518.                m1=num(1,i1)
  519.                m2=num(2,i1)
  520.                m3=num(3,i1)
  521. C     CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2)
  522.                xm1=xcoor3((m1-1)*(idim+1)+1)
  523.                ym1=xcoor3((m1-1)*(idim+1)+2)
  524.                zm1=xcoor3((m1-1)*(idim+1)+3)
  525.                xm2=xcoor3((m2-1)*(idim+1)+1)
  526.                ym2=xcoor3((m2-1)*(idim+1)+2)
  527.                zm2=xcoor3((m2-1)*(idim+1)+3)
  528.                xm3=xcoor3((m3-1)*(idim+1)+1)
  529.                ym3=xcoor3((m3-1)*(idim+1)+2)
  530.                zm3=xcoor3((m3-1)*(idim+1)+3)
  531.                q(1)=(xcoor4((ia-1)*(idim+1)+1))
  532.                q(2)=(xcoor4((ia-1)*(idim+1)+2))
  533.                q(3)=(xcoor4((ia-1)*(idim+1)+3))
  534. C            Calcul du zib
  535.                o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  536.                o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  537.                o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  538.                o4=-(o1*xm1+o2*ym1+o3*zm1)
  539.                zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  540.                xnew=q(1)*zk
  541.                ynew=q(2)*zk
  542.                znew=q(3)*zk
  543. C   PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1
  544.                r1=sqrt(o1*o1+o2*o2+o3*o3)
  545.                if( r1 . EQ. 0. ) then
  546.                  call erreur ( 21 )
  547.               endif
  548.               o1 = o1/r1
  549.               o2 = o2/r1
  550.               o3 = o3/r1
  551.               r2=sqrt(o1*o1+o2*o2)
  552.               if ( r2 . gt. 1.e-5) then
  553.                 tabmat(1,1)=-o2/r2
  554.                 tabmat(1,2)=o1/r2
  555.                 tabmat(1,3)=0.
  556.               else
  557.                 tabmat(1,1)=1.
  558.                 tabmat(1,2)=0.
  559.                 tabmat(1,3)=0.
  560.               endif
  561.             tabmat(2,1)=-o3*tabmat(1,2)
  562.             tabmat(2,2)=o3*tabmat(1,1)
  563.             tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1)
  564.             tabmat(3,1)=o1
  565.             tabmat(3,2)=o2
  566.             tabmat(3,3)=o3
  567.         x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3)))
  568.         y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3)))
  569.         x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3)))
  570.         y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3)))
  571.         x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3)))
  572.         y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3)))
  573.         x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3)))
  574.         y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3)))
  575. C    CALCUL DES COORDONNEES BARYCENTRIQUES
  576.             detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  577.             a(1)=x2*y3-y2*x3
  578.             a(2)=x3*y1-x1*y3
  579.             a(3)=x1*y2-x2*y1
  580.             b(1)=y2-y3
  581.             b(2)=y3-y1
  582.             b(3)=y1-y2
  583.             c(1)=x3-x2
  584.             c(2)=x1-x3
  585.             c(3)=x2-x1
  586.             Do 798 k=1,nbnn
  587.                al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  588. 798         continue
  589. C    TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  590.            If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  591.      $     al(3).gt.prec3) then
  592. C    LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  593.                 icor(1)=icor(1)+1
  594.                 icor(icor(1)+1)=i1
  595.            Endif
  596. 1462       continue
  597.          Endif
  598.  
  599.          If (icor(1).eq.0) then
  600.            interr(1) = ia
  601.            call erreur(782)
  602.            return
  603.          Endif
  604.  
  605.       Endif
  606. C   FIN DU IF ICPR2(IA).EQ.1
  607.  
  608.          If (icpr2(ia).eq.3) then
  609.  
  610. C   PASSAGE SUR LES ELEMENTS DE LA TRANCHE CRITIQUE
  611.              Do 1463 iz=1,tabcrit(1)
  612.              i1=tabcrit(iz+1)
  613.              m1=num(1,i1)
  614.              m2=num(2,i1)
  615.              m3=num(3,i1)
  616. C     CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2)
  617.              xm1=xcoor3((m1-1)*(idim+1)+1)
  618.              ym1=xcoor3((m1-1)*(idim+1)+2)
  619.              zm1=xcoor3((m1-1)*(idim+1)+3)
  620.              xm2=xcoor3((m2-1)*(idim+1)+1)
  621.              ym2=xcoor3((m2-1)*(idim+1)+2)
  622.              zm2=xcoor3((m2-1)*(idim+1)+3)
  623.              xm3=xcoor3((m3-1)*(idim+1)+1)
  624.              ym3=xcoor3((m3-1)*(idim+1)+2)
  625.              zm3=xcoor3((m3-1)*(idim+1)+3)
  626.              q(1)=(xcoor4((ia-1)*(idim+1)+1))
  627.              q(2)=(xcoor4((ia-1)*(idim+1)+2))
  628.              q(3)=(xcoor4((ia-1)*(idim+1)+3))
  629. C            Calcul du zib
  630.              o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  631.              o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  632.              o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  633.              o4=-(o1*xm1+o2*ym1+o3*zm1)
  634.              zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  635.              xnew=q(1)*zk
  636.              ynew=q(2)*zk
  637.              znew=q(3)*zk
  638. C   PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1
  639.         r1=sqrt(o1*o1+o2*o2+o3*o3)
  640.         if( r1 . EQ. 0. ) then
  641.           call erreur ( 21 )
  642.         endif
  643.         o1 = o1/r1
  644.         o2 = o2/r1
  645.         o3 = o3/r1
  646.         r2=sqrt(o1*o1+o2*o2)
  647.         if ( r2 . gt. 1.e-5) then
  648.         tabmat(1,1)=-o2/r2
  649.         tabmat(1,2)=o1/r2
  650.         tabmat(1,3)=0.
  651.         else
  652.         tabmat(1,1)=1.
  653.         tabmat(1,2)=0.
  654.         tabmat(1,3)=0.
  655.         endif
  656.         tabmat(2,1)=-o3*tabmat(1,2)
  657.         tabmat(2,2)=o3*tabmat(1,1)
  658.         tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1)
  659.         tabmat(3,1)=o1
  660.         tabmat(3,2)=o2
  661.         tabmat(3,3)=o3
  662.         x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3)))
  663.         y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3)))
  664.         x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3)))
  665.         y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3)))
  666.         x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3)))
  667.         y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3)))
  668.         x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3)))
  669.         y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3)))
  670. C    CALCUL DES COORDONNEES BARYCENTRIQUES
  671.         detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  672.                a(1)=x2*y3-y2*x3
  673.                a(2)=x3*y1-x1*y3
  674.                a(3)=x1*y2-x2*y1
  675.                b(1)=y2-y3
  676.                b(2)=y3-y1
  677.                b(3)=y1-y2
  678.                c(1)=x3-x2
  679.                c(2)=x1-x3
  680.                c(3)=x2-x1
  681.                Do 799 k=1,nbnn
  682.                   al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  683. 799            continue
  684. C    TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  685.              If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  686.      $al(3).gt.prec3) then
  687. C    LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  688.              icor(1)=icor(1)+1
  689.              icor(icor(1)+1)=i1
  690.              Endif
  691. 1463       continue
  692.          If (icor(1).eq.0) then
  693.            interr(1)=ia
  694.            call erreur(782)
  695.             return
  696.          Endif
  697.          Endif
  698. C    FIN DE L'INVENTAIRE
  699.  
  700. C    CALCUL DES PROJETES DE MAI2 SUR MAI1
  701.          If  ((icpr2(ia).eq.1) .or. (icpr2(ia).eq.3)) then
  702.  
  703.               dmin=1.e30
  704.               nbpts=nbpts+1
  705.               segadj mcoord
  706.               n=0
  707.               Do 500 m=1,icor(1)
  708.                 m1=num(1,icor(m+1))
  709.                 m2=num(2,icor(m+1))
  710.                 m3=num(3,icor(m+1))
  711. C               Calcul des coordonnees du projete
  712.                 xm1=xcoor3((m1-1)*(idim+1)+1)
  713.                 ym1=xcoor3((m1-1)*(idim+1)+2)
  714.                 zm1=xcoor3((m1-1)*(idim+1)+3)
  715.                 xm2=xcoor3((m2-1)*(idim+1)+1)
  716.                 ym2=xcoor3((m2-1)*(idim+1)+2)
  717.                 zm2=xcoor3((m2-1)*(idim+1)+3)
  718.                 xm3=xcoor3((m3-1)*(idim+1)+1)
  719.                 ym3=xcoor3((m3-1)*(idim+1)+2)
  720.                 zm3=xcoor3((m3-1)*(idim+1)+3)
  721.                 q(1)=(xcoor4((ia-1)*(idim+1)+1))
  722.                 q(2)=(xcoor4((ia-1)*(idim+1)+2))
  723.                 q(3)=(xcoor4((ia-1)*(idim+1)+3))
  724. C               Calcul du zib
  725.                 o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  726.                 o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  727.                 o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  728.                 o4=-(o1*xm1+o2*ym1+o3*zm1)
  729.                 zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  730.                 xnew=q(1)*zk
  731.                 ynew=q(2)*zk
  732.                 znew=q(3)*zk
  733.                 delta1=sqrt(xnew*xnew+ynew*ynew+znew*znew)
  734.                 delta2=sqrt(q(1)*q(1)+q(2)* q(2)+q(3)*q(3))
  735.                 pds=(xnew*q(1))+(ynew*q(2))+(znew*q(3))
  736.                 If (pds.gt.0.) then
  737.                  If ((delta1-delta2).gt.0.) then
  738.                   If ((delta1-delta2).lt.dmin) then
  739.                        icorre(ia) = nbpts
  740.                        dmin=(delta1-delta2)
  741.                        xcoor((nbpts-1)*(idim+1)+1)=xnew
  742.                        xcoor((nbpts-1)*(idim+1)+2)=ynew
  743.                        xcoor((nbpts-1)*(idim+1)+3)=znew
  744.                   Endif
  745.                   n=n+1
  746.                 Endif
  747.                 Endif
  748. 500           continue
  749.            If (n.eq.0) then
  750.            call erreur(21)
  751.            return
  752.            Endif
  753.  
  754. C *********************RETOUR A L'ANCIENNE BASE*********************
  755.         ib=nbpts
  756.         xa1=xcoor((ib-1)*(idim+1)+1)
  757.         xa2=xcoor((ib-1)*(idim+1)+2)
  758.         xa3=xcoor((ib-1)*(idim+1)+3)
  759.         do 3200 j=1,3
  760.         xcoor((ib-1)*(idim+1)+j)=(xa1*invma(j,1))+
  761.      $(xa2*invma(j,2))+(xa3*invma(j,3))+w(J)
  762. 3200               continue
  763.         icor(1)=0
  764.         n=0
  765.         Endif
  766. 405     continue
  767.         segsup iseg1,iseg3,iseg4,iseg5,iseg6,mcoora,mcoor2
  768.         segsup icp1,icp2,mcoor3,mcoor4
  769. C   CREATION DU MAILLAGE PROJETE MAI3
  770.         segini,ipt3=ipt1
  771.         ipt4 = IPT3
  772.         do 2456 IU = 1, max( 1,IPT3.LISOUS(/1))
  773.         If( IPT3.LISOUS(/1).NE.0) THEN
  774.           IPT5=IPT3.LISOUS(IU)
  775.           SEGINI,IPT4=IPT5
  776.           NBREF = 0
  777.           NBSOUS=0
  778.           SEGADJ IPT4
  779.         ENDIF
  780.         nbel1=ipt4.num(/2)
  781.         nbnn1=ipt4.num(/1)
  782.         Do 600 i=1,nbel1
  783.             Do 601 j=1,nbnn1
  784.             ia=ipt4.num(j,i)
  785.             ib=icorre(ia)
  786.             ipt4.num(j,i)=ib
  787. 601         continue
  788. 600     continue
  789. 2456    continue
  790.         do 2457 IU = 1, max( 0,IPT3.LISREF(/1))
  791.         If( IPT3.LISREF(/1).NE.0) THEN
  792.          IPT5 = IPT3.LISREF(IU)
  793.          SEGINI,IPT4=IPT5
  794.         ENDIF
  795.         nbel1=ipt4.num(/2)
  796.         nbnn1=ipt4.num(/1)
  797.         Do 602 i=1,nbel1
  798.             Do 603 j=1,nbnn1
  799.             ia=ipt4.num(j,i)
  800.             ib=icorre(ia)
  801.             ipt4.num(j,i)=ib
  802. 603         continue
  803. 602     continue
  804.         SEGDES IPT4
  805. 2457    continue
  806.         call ecrobj('MAILLAGE',ipt3)
  807.         segdes ipt3
  808.         segsup mcor,mcorre,matric
  809.         Return
  810.         End
  811.  
  812.  
  813.  
  814.  
  815.  
  816.  
  817.  
  818.  
  819.  
  820.  
  821.  
  822.  
  823.  
  824.  
  825.  

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