C PROOB1    SOURCE    SP204843  24/08/05    21:15:02     11978          
C CE PROGRAMME PERMET LA PROJECTION D'UN MAILLAGE 3D SUR UNE SUFACE
C CARACTERISEE PAR UN AUTRE MAILLAGE 3D.

C ENTREES  MAI2  : UN MAILLAGE DE TYPE MELEME A PROJETER.
C          MAI1  : UN MAILLAGE REPRESENTANT LA SURFACE SUR LAQUELLE MAI2
C                  EST PROJETE.
C          IPOIN : UN POINT S' IL S'AGIT D'UN POINT A PROJETER.
C          IP1   : UN VECTEUR DEFINISSANT LA DIRECTION DE PROJECTION.

C DATE: AOUT 1996.

C   **********************************************************************
        Subroutine proob1(mai2,mai1,ip1)
      IMPLICIT INTEGER(I-N)
        IMPLICIT REAL*8(A-H,O-Z)


-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMCOORD
-INC CCREEL
        SEGMENT ISEG1
        REAL*8 XLIM(2,NBEL),YLIM(2,NBEL),ZLIM(2,NBEL)
        ENDSEGMENT

        SEGMENT ISEG3
        INTEGER NIZO(NZO+1)
        ENDSEGMENT

        SEGMENT ISEG4
        INTEGER NUMZO(NZO)
        ENDSEGMENT

        SEGMENT ISEG5
        INTEGER NNMEL(ILON),IDEJ(NZO)
        ENDSEGMENT

        SEGMENT ISEG6
        REAL*8 AM(4,4),AL(4),A(4),B(4),C(4)
        REAL*8 XPU(2),P(3)
        ENDSEGMENT

        SEGMENT MCOORA
        REAL*8  XCOORA(XCOOR(/1))
        ENDSEGMENT
C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.

        SEGMENT MCOOR2
        REAL*8  XCOOR2(XCOOR(/1))
        ENDSEGMENT
C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.

        SEGMENT MCOR
        INTEGER ICOR(IONMAX+1)
        ENDSEGMENT

        SEGMENT MCORRE
        INTEGER ICORRE((nbpts))
        ENDSEGMENT

        SEGMENT ICP1
        INTEGER ICPR1(nbpts)
        ENDSEGMENT

        SEGMENT ICP2
        INTEGER  ICPR2(nbpts)
        ENDSEGMENT

        SEGMENT MATRIC
        REAL*8 TABMAT(3,3),invma(3,3)
        ENDSEGMENT

C PASSAGE DU MAILLAGE MAI1 EN TRI3
        call change(mai1,4)
        if (ierr.ne.0) return
C CREATION DE LA MATRICE DE CHANGEMENT DE BASE
        segini matric
        v1=xcoor((ip1-1)*(idim+1)+1)
        v2=xcoor((ip1-1)*(idim+1)+2)
        v3=xcoor((ip1-1)*(idim+1)+3)
        r1=sqrt(v1*v1+v2*v2+v3*v3)
        if( r1 . EQ. 0. ) then
          call erreur ( 21 )
        endif
        v1 = V1 / r1
        v2 = v2 / r1
        v3 = v3 / r1
        r2=sqrt(v1*v1+v2*v2)
        if ( r2 . gt. 1.e-5) then
        tabmat(1,1)=-v2/r2
        tabmat(1,2)=v1/r2
        tabmat(1,3)=0.
        else
        tabmat(1,1)=1.
        tabmat(1,2)=0.
        tabmat(1,3)=0.
        endif
        tabmat(2,1)=-v3*tabmat(1,2)
        tabmat(2,2)=v3*tabmat(1,1)
        tabmat(2,3)=v1*tabmat(1,2)-v2*tabmat(1,1)
        tabmat(3,1)=v1
        tabmat(3,2)=v2
        tabmat(3,3)=v3
* pas besoin de diviser par le determinant car egal à 1.
        INVMA(1,1)= TABMAT(2,2)*TABMAT(3,3)-TABMAT(3,2)*TABMAT(2,3)
        INVMA(1,2)=-TABMAT(1,2)*TABMAT(3,3)+TABMAT(3,2)*TABMAT(1,3)
        INVMA(1,3)= TABMAT(1,2)*TABMAT(2,3)-TABMAT(2,2)*TABMAT(1,3)
        INVMA(2,1)=-TABMAT(2,1)*TABMAT(3,3)+TABMAT(3,1)*TABMAT(2,3)
        INVMA(2,2)= TABMAT(1,1)*TABMAT(3,3)-TABMAT(3,1)*TABMAT(1,3)
        INVMA(2,3)=-TABMAT(1,1)*TABMAT(2,3)+TABMAT(2,1)*TABMAT(1,3)
        INVMA(3,1)= TABMAT(2,1)*TABMAT(3,2)-TABMAT(2,2)*TABMAT(3,1)
        INVMA(3,2)=-TABMAT(1,1)*TABMAT(3,2)+TABMAT(3,1)*TABMAT(1,2)
        INVMA(3,3)= TABMAT(1,1)*TABMAT(2,2)-TABMAT(2,1)*TABMAT(1,2)

C LISTAGE DES POINTS APPARTENANT A MAI1
        segini icp1
        meleme=mai1
        segact meleme
        nbelem=num(/2)
        nbnn = num(/1)
        Do 10 i=1,nbelem
           Do 20 j=1,nbnn
              ia=num(j,i)
              icpr1(ia)=1
20          continue
10      continue

C LISTAGE DES POINTS APPARTENANT A MAI2
        ipt1=mai2
        segact ipt1
        segini icp2
        nbelem= ipt1.num(/2)
        nbnn = ipt1.num(/1)
        Do 30 i=1,nbelem
           Do 40 j=1,nbnn
              ia=ipt1.num(j,i)
              If (icpr1(ia).eq.1) then
              icpr2(ia)=2
              else
              icpr2(ia)=1
              Endif
40          continue
30      continue

C CALCUL DES COORDONNEES DANS LE NOUVEAU REPERE
        segini mcoora
        Do 31 ia=1,nbpts
              If (icpr1(ia).eq.1) then
                 xia1=xcoor((ia-1)*(idim+1)+1)
                 xia2=xcoor((ia-1)*(idim+1)+2)
                 xia3=xcoor((ia-1)*(idim+1)+3)
                 Do 32 j=1,3
                    xcoora((ia-1)*(idim+1)+j)=(xia1*tabmat(j,1))+
     $(xia2*tabmat(j,2))+(xia3*tabmat(j,3))
32               continue
              Endif
31      continue

        segini mcoor2
        Do 33 ia=1,nbpts
              If (icpr2(ia).eq.1)  then
                 xia1=xcoor((ia-1)*(idim+1)+1)
                 xia2=xcoor((ia-1)*(idim+1)+2)
                 xia3=xcoor((ia-1)*(idim+1)+3)
                 Do 34 j=1,3
                    xcoor2((ia-1)*(idim+1)+j)=(xia1*tabmat(j,1))+
     $(xia2*tabmat(j,2))+(xia3*tabmat(j,3))
34               continue
              Endif
              If (icpr2(ia).eq.2)  then
                 xia1=xcoor((ia-1)*(idim+1)+1)
                 xia2=xcoor((ia-1)*(idim+1)+2)
                 xia3=xcoor((ia-1)*(idim+1)+3)
                 Do 35 j=1,3
                    xcoor2((ia-1)*(idim+1)+j)=(xia1*tabmat(j,1))+
     $(xia2*tabmat(j,2))+(xia3*tabmat(j,3))
35               continue
              Endif
33      continue

C******************** DEBUT DU ZONAGE DE MAI1 **********************

C ON CALCULE LA TAILLE MAXI D'UN ELEMENT DANS TOUTES LES DIRECTIONS
C AFIN DE CREER UN ZONAGE DE L'ESPACE. EN MEME TEMPS ON CALCULE
C LA DIMENSION HORS TOUT DU MAILLAGE
C
      IDIM1=4
      NBEL = NUM(/2)
      NBNN=NUM(/1)
      SEGINI ISEG1
      ILOC=0
      XZO=0.D0
      YZO=0.D0
      ZZO=0.D0
      XZA=XGRAND
      YZA=XGRAND
      ZZA=XGRAND
      XTOMI=XGRAND
      XTOMA=-XGRAND
      YTOMI=XGRAND
      YTOMA=-XGRAND
      ZTOMI=XGRAND
      ZTOMA=-XGRAND
      DO 1 I1=1,NBEL
         XMI=XGRAND
         YMI=XGRAND
         ZMI=XGRAND
         YMA=-XGRAND
         XMA=-XGRAND
         ZMA=-XGRAND
         DO 2 I2 = 1,NBNN
            IB=NUM(I2,I1)
            IA=(IB-1)*IDIM1
            IF(XCOORA(IA+1).LT.XMI) XMI=XCOORA(IA+1)
            IF(XCOORA(IA+1).GT.XMA) XMA=XCOORA(IA+1)
            IF(XCOORA(IA+2).LT.YMI) YMI=XCOORA(IA+2)
            IF(XCOORA(IA+2).GT.YMA) YMA=XCOORA(IA+2)
    2    CONTINUE
         XLIM(1,I1)=XMI
         XLIM(2,I1)=XMA
         YLIM(1,I1)=YMI
         YLIM(2,I1)=YMA
         XZO=MAX (XZO,XMA-XMI)
         YZO=MAX (YZO,YMA-YMI)
         XZA=MIN(XZA,XMA-XMI)
         YZA=MIN(YZA,YMA-YMI)
         IF(XMI.LT.XTOMI) XTOMI=XMI
         IF(XMA.GT.XTOMA) XTOMA=XMA
         IF(YMI.LT.YTOMI) YTOMI=YMI
         IF(YMA.GT.YTOMA) YTOMA=YMA

    1 CONTINUE
      XPR=MIN(XZO*1.D-2,(XTOMA-XTOMI)/2.D+4)
      YPR=MIN(YZO*1.D-2,(YTOMA-YTOMI)/2.D+4)
      XZA=XZA*0.97
      YZA=YZA*0.97
      XTOMI= XTOMI - (XTOMA-XTOMI)/1.D+4
      XTOMA= XTOMA + (XTOMA-XTOMI)/1.D+4
      YTOMI= YTOMI - (YTOMA-YTOMI)/1.D+4
      YTOMA= YTOMA + (YTOMA-YTOMI)/1.D+4
      XZA=MIN ( XZA, XTOMA-XTOMI)
      YZA=MIN ( YZA, YTOMA-YTOMI)
      XZA = max ( xza , (XTOMA-XTOMI) / 50)
      YZA = max ( yza , (YTOMA-YTOMI) / 50)
      NXZO=1
      if (abs(XZA).GT.((XMI+XMA)*XZPREC)) NXZO=int((XTOMA-XTOMI)/XZA)+1
      if (abs(YZA).GT.((YMI+YMA)*XZPREC)) NYZO=int((YTOMA-YTOMI)/YZA)+1
      XZO=XZA
      YZO=YZA
      NZZO=1
       NXDEP=MIN(NXZO,10)
       NYDEP=MIN(NYZO,10)
         IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN
          XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/90
          NXZO=MAX(INT(NXZO/XY),NXDEP)
          NYZO=MAX(INT(NYZO/XY),NYDEP)
          IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN
            XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/60
            NXZO=MAX(INT(NXZO/XY),NXDEP)
            NYZO=MAX(INT(NYZO/XY),NYDEP)
          ENDIF
          XZO=(XTOMA-XTOMI)/NXZO
          YZO=(YTOMA-YTOMI)/NYZO
          NXZO=int((XTOMA-XTOMI)/XZO) +1
          NYZO=int((YTOMA-YTOMI)/YZO) +1
         ENDIF

C
C  ON VEUT CONSTRUIRE LA LISTE DES ELEMENTS TOUCHANT UNE ZONE
C  POUR CELA ON COMMENCE PAR COMPTER COMBIEN D'ELEMENT TOUCHENT
C  CHAQUE ZONE ET EN MEME TEMPS ON STOCKE LES ZONES TOUCHEES
C  PAR CHAQUE ELEMENT ET LEUR NOMBRE
C

      NZO=NXZO*NYZO
      IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NZO NXZO NYZO NZZO  ''
     $,4I7) ') NZO,NXZO,NYZO,NZZO
      NXYZO=NXZO*NYZO
      if (nzo.lt.0) then
        call erreur(21)
        return
      endif
      SEGINI ISEG3
      SEGINI ISEG4
      DO 3 I1=1,NBEL
          NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
          NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
          NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
          NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
               DO 201 L1=NIZ1Y,NIZ2Y
               DO 2011 L2=NIZ1X,NIZ2X
                   NIZA = L2 + ( L1-1) * NXZO
                   NUMZO(NIZA) = NUMZO(NIZA) +1
 2011          CONTINUE
  201          CONTINUE
    3 CONTINUE
C
C  CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES
C  ELEMENTS CONCERNEES PAR UNE ZONE
C
      ILON=0
      NIZO(1)=1
      DO 202 L1=1,NZO
      NIZO(L1+1)=NIZO(L1)+NUMZO(L1)
      ILON=ILON+ NUMZO(L1)
 202  CONTINUE
 110  FORMAT(16I5)
      SEGINI ISEG5
      DO 5 I1=1,NBEL
          NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) + 1
          NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) + 1
          NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) + 1
          NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) + 1
               DO 203 L1=NIZ1Y,NIZ2Y
               DO 2031 L2=NIZ1X,NIZ2X
                   NIZA = L2 + ( L1-1) * NXZO
                   IAD=NIZO(NIZA)+IDEJ(NIZA)
                   NNMEL(IAD)=I1
                   IDEJ(NIZA)=IDEJ(NIZA)+1
 2031          CONTINUE
  203          CONTINUE
    5  CONTINUE
C *********************** FIN DU ZONAGE DE MAI1 **************************
        prec = 1.E-5
        prec2 = -prec
        prec3 = -2.E-3
C ******RECHERCHE DU NOMBRE D'ELEMENT MAX CONTENU DANS UNE ZONE***********
        la=0
        ls=0
        Do 420 i=1,nzo
          ls=ls+numzo(i)
          If (numzo(i).gt.la) then
          la=numzo(i)
          Endif
420     continue
        ionmax=la

C ****************** TRAITEMENT SUR LES POINTS DE MAI2 *******************
        nbelem=num(/2)
        nbnn=num(/1)
        nbpts2=nbpts
        segini mcor,mcorre,iseg6
        Do 405 ia=1,nbpts2
           If (icpr2(ia).eq.2) then
           icorre(ia)=ia
           Endif
           If (icpr2(ia).eq.1) then
           xpu(1)=xcoor2((ia-1)*(idim+1)+1)
           xpu(2)=xcoor2((ia-1)*(idim+1)+2)
         if (((xpu(1)-xtomi).lt.prec.) .or.
     $ ((xpu(2)-ytomi).lt.prec.) .or.
     $ ((xpu(1)-xtoma).gt.prec2.) .or.
     $ ((xpu(2)-ytoma).gt.prec2.)) then
            call erreur(21)
            return
         endif
C    RECHERCHE DU NUMERO DE LA ZONE CORRESPONDANTE
         k2=int(((xcoor2((ia-1)*(idim+1)+1))-xtomi-xpr)/xzo)+1
         k1=int(((xcoor2((ia-1)*(idim+1)+2))-ytomi-ypr)/yzo)+1
           niza=k2+(k1-1)*nxzo
C    INVENTAIRE DES ELEMENTS CONCERNES PAR LE POINT IA
           iad=nizo(niza)
           Do 406 i=1,numzo(niza)
             i1=nnmel(iad)
C    CALCUL DES COORDONNEES BARYCENTRIQUES DE IA DANS L'ELEMENT I1
             j1=num(1,i1)
             j2=num(2,i1)
             j3=num(3,i1)
             j1idim=(j1-1)*(idim+1)
             j2idim=(j2-1)*(idim+1)
             j3idim=(j3-1)*(idim+1)
             Do 408 l=1,2
                p(l+1)=xpu(l)
                am(l+1,1)=xcoora(j1idim+l)
                am(l+1,2)=xcoora(j2idim+l)
                am(l+1,3)=xcoora(j3idim+l)
408          continue
               x1=am(2,1)
               x2=am(2,2)
               x3=am(2,3)
               y1=am(3,1)
               y2=am(3,2)
               y3=am(3,3)
               x=p(2)
               y=p(3)
               detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
               a(1)=x2*y3-y2*x3
               a(2)=x3*y1-x1*y3
               a(3)=x1*y2-x2*y1
               b(1)=y2-y3
               b(2)=y3-y1
               b(3)=y1-y2
               c(1)=x3-x2
               c(2)=x1-x3
               c(3)=x2-x1
               Do 409 k=1,nbnn
                  al(k)=(a(k)+b(k)*x+c(k)*y)/detam
409            continue
C    TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
             If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
     $al(3).gt.prec3) then
C    LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
                    icor(1)=icor(1)+1
                    icor(icor(1)+1)=i1
             Endif
                   iad=iad+1

406        continue
         If (icor(1).eq.0) then
           interr(1) = ia
           call erreur(782)
            return
         Endif
         Endif
C    FIN DE L'INVENTAIRE
C    CALCUL DES PROJETES DE MAI2 SUR MAI1
         If  (icpr2(ia).eq.1) then
C    SI IA SE PROJETTE SUR PLUSIEURS ELEMENTS
              dmin=1.e30
              nbpts=nbpts+1
              segadj mcoord
              Do 500 m=1,icor(1)
                m1=num(1,icor(m+1))
                m2=num(2,icor(m+1))
                m3=num(3,icor(m+1))
C               Calcul des coordonnees du projete
                xm1=xcoora((m1-1)*(idim+1)+1)
                ym1=xcoora((m1-1)*(idim+1)+2)
                zm1=xcoora((m1-1)*(idim+1)+3)
                xm2=xcoora((m2-1)*(idim+1)+1)
                ym2=xcoora((m2-1)*(idim+1)+2)
                zm2=xcoora((m2-1)*(idim+1)+3)
                xm3=xcoora((m3-1)*(idim+1)+1)
                ym3=xcoora((m3-1)*(idim+1)+2)
                zm3=xcoora((m3-1)*(idim+1)+3)
                xnew=xcoor2((ia-1)*(idim+1)+1)
                ynew=xcoor2((ia-1)*(idim+1)+2)
C               Calcul du zib
                o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
                o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
                o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
                o4=-(o1*xm1+o2*ym1+o3*zm1)
                znew=-(o1*xnew+o2*ynew+o4)/o3
                zpoint = xcoor2((ia-1)*(idim+1)+3)
*               If ((znew-zpoint).gt.0.) then
                  If ((znew-zpoint).lt.dmin) then
                       icorre(ia) = nbpts
                       dmin=znew-zpoint
                       xcoor((nbpts-1)*(idim+1)+1)=xnew
                       xcoor((nbpts-1)*(idim+1)+2)=ynew
                       xcoor((nbpts-1)*(idim+1)+3)=znew
                  Endif
*                Endif
500           continue

C *********************RETOUR A L'ANCIENNE BASE*********************
        ib=nbpts
        xa1=xcoor((ib-1)*(idim+1)+1)
        xa2=xcoor((ib-1)*(idim+1)+2)
        xa3=xcoor((ib-1)*(idim+1)+3)
                 Do 3200 j=1,3
                    xcoor((ib-1)*(idim+1)+j)=(xa1*invma(j,1))+
     $(xa2*invma(j,2))+(xa3*invma(j,3))
3200               continue
         Endif
        icor(1)=0
405     continue

        segsup iseg1,iseg3,iseg4,iseg5,iseg6,mcoora,mcoor2
        segsup icp1,icp2
C   CREATION DU MAILLAGE PROJETE MAI3
        segini,ipt3=ipt1
        ipt4 = IPT3
        do 2456 IU = 1, max( 1,IPT3.LISOUS(/1))
        If( IPT3.LISOUS(/1).NE.0) THEN
          IPT5=IPT3.LISOUS(IU)
          SEGINI,IPT4=IPT5
          NBREF = 0
          NBSOUS=0
          SEGADJ IPT4
        ENDIF
        nbel1=ipt4.num(/2)
        nbnn1=ipt4.num(/1)
        Do 600 i=1,nbel1
            Do 601 j=1,nbnn1
            ia=ipt4.num(j,i)
            ib=icorre(ia)
            ipt4.num(j,i)=ib
601         continue
600     continue
2456    continue
        do 2457 IU = 1, max( 0,IPT3.LISREF(/1))
        If( IPT3.LISREF(/1).NE.0) THEN
         IPT5 = IPT3.LISREF(IU)
         SEGINI,IPT4=IPT5
        ENDIF
        nbel1=ipt4.num(/2)
        nbnn1=ipt4.num(/1)
        Do 602 i=1,nbel1
            Do 603 j=1,nbnn1
            ia=ipt4.num(j,i)
            ib=icorre(ia)
            ipt4.num(j,i)=ib
603         continue
602     continue
        SEGDES IPT4
2457    continue
        call ecrobj('MAILLAGE',ipt3)
        segdes ipt3,ipt1,meleme
        segsup mcor,mcorre,matric
        Return
        End