C PROOB4 SOURCE PV090527 24/06/13 08:10:47 11941 C CE PROGRAMME PERMET LA PROJECTION CONIQUE D'UN MAILLAGE 3D SUR C UNE SUFACE 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 IP1 : UN POINT DEFINISSANT L'OEIL DE LA PROJECTION. C DATE: SEPTEMBRE 1996. C ********************************************************************** Subroutine proob4(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 MCOOR3 REAL*8 XCOOR3(XCOOR(/1)) ENDSEGMENT C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE. SEGMENT MCOOR4 REAL*8 XCOOR4(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),V(3),W(3),Q(3),invma(3,3) ENDSEGMENT SEGMENT CRIT INTEGER ELECRI(NBELEM),TABCRIT(NBELEM+1) ENDSEGMENT segini matric segact mcoord*mod zchoix=100. w(1)=xcoor((ip1-1)*(idim+1)+1) w(2)=xcoor((ip1-1)*(idim+1)+2) w(3)=xcoor((ip1-1)*(idim+1)+3) C PASSAGE DU MAILLAGE MAI1 EN TRI3 call change(mai1,4) 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 POSITION DU CENTRE DE GRAVITE DE MAI2 sx=0. sy=0. sz=0. isd=0 Do 2875 ia=1,nbpts If (icpr2(ia).ne.0) then sx=sx+xcoor((ia-1)*(idim+1)+1) sy=sy+xcoor((ia-1)*(idim+1)+2) sz=sz+xcoor((ia-1)*(idim+1)+3) isd = isd+1 Endif 2875 continue C CREATION DE LA MATRICE DE CHANGEMENT DE BASE v(1)=(sx/isd)-w(1) v(2)=(sy/isd)-w(2) v(3)=(sz/isd)-w(3) If ((abs(v(1))).lt.(1.E-3) .and. (abs(v(2))).lt.(1.E-3) .and. $ (abs(v(3))).lt.(1.E-3)) then v(1)=v(1)+0.1 v(2)=v(2)+0.1 v(3)=v(3)+0.1 Endif segini mcoor3 segini mcoora segini mcoor4 segini mcoor2 r1=sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) v(1) = V(1) / r1 v(2) = v(2) / r1 v(3) = v(3) / r1 r2=sqrt(v(1)*v(1)+v(2)*v(2)) if ( r2 . gt. 1.e-9) then tabmat(1,1)=-v(2)/r2 tabmat(1,2)=v(1)/r2 tabmat(1,3)=0. else tabmat(1,1)=1. tabmat(1,2)=0. tabmat(1,3)=0. endif tabmat(2,1)=-v(3)*tabmat(1,2) tabmat(2,2)=v(3)*tabmat(1,1) tabmat(2,3)=v(1)*tabmat(1,2)-v(2)*tabmat(1,1) tabmat(3,1)=v(1) tabmat(3,2)=v(2) tabmat(3,3)=v(3) * 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 CALCUL DES COORDONNEES DANS LE NOUVEAU REPERE POUR MELEME nbelem=num(/2) nbnn=num(/1) segini crit Do 31 i1=1,nbelem Do 38 k=1,nbnn ia=num(k,i1) xia1=xcoor((ia-1)*(idim+1)+1)-w(1) xia2=xcoor((ia-1)*(idim+1)+2)-w(2) xia3=xcoor((ia-1)*(idim+1)+3)-w(3) Do 32 j=1,3 xcoor3((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+ $(xia2*tabmat(j,2))+(xia3*tabmat(j,3))) 32 continue 38 continue 31 continue C CALCUL DE LA TAILLE LA PLUS GRANDE DES ELEMENTS SUIVANT OZ zlong=-xgrand Do 1534 i1=1,nbelem zmin=xgrand zmax=-xgrand Do 1233 i2=1,nbnn ib=num(i2,i1) ia=(ib-1)*4 If (xcoor3(ia+3).lt.zmin) zmin=xcoor3(ia+3) If (xcoor3(ia+3).gt.zmax) zmax=xcoor3(ia+3) 1233 continue zlong=max(zlong,zmax-zmin) 1534 continue If (zlong.lt.(5.E-2)) zlong=(5.E-2) C DETERMINATION DES ELEMENTS TOUCHANT LA TRANCHE CRITIQUE Do 317 i1=1,nbelem ime=0 Do 387 k=1,nbnn ia=num(k,i1) q(1)=(xcoor3((ia-1)*(idim+1)+1)) q(2)=(xcoor3((ia-1)*(idim+1)+2)) q(3)=(xcoor3((ia-1)*(idim+1)+3)) If (abs(xcoor3((ia-1)*(idim+1)+3)).gt.(0.55*zlong)) then xcoora((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1) xcoora((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2) xcoora((ia-1)*(idim+1)+3)=zchoix else ime=ime+1 Endif 387 continue If (ime.ne.0) then tabcrit(1)=tabcrit(1)+1 tabcrit(tabcrit(1)+1)=i1 elecri(i1)=1 Endif 317 continue C CHANGEMENT DE REPERE POUR IPT1 Do 33 ia=1,nbpts If (icpr2(ia).eq.1) then xia1=xcoor((ia-1)*(idim+1)+1)-w(1) xia2=xcoor((ia-1)*(idim+1)+2)-w(2) xia3=xcoor((ia-1)*(idim+1)+3)-w(3) Do 34 j=1,3 xcoor4((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+ $(xia2*tabmat(j,2))+(xia3*tabmat(j,3))) 34 continue If (abs(xcoor4((ia-1)*(idim+1)+3)).lt.(0.55*zlong)) then icpr2(ia)=3 else q(1)=(xcoor4((ia-1)*(idim+1)+1)) q(2)=(xcoor4((ia-1)*(idim+1)+2)) q(3)=(xcoor4((ia-1)*(idim+1)+3)) xcoor2((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1) xcoor2((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2) xcoor2((ia-1)*(idim+1)+3)=zchoix Endif Endif If (icpr2(ia).eq.2) then xia1=xcoor((ia-1)*(idim+1)+1)-w(1) xia2=xcoor((ia-1)*(idim+1)+2)-w(2) xia3=xcoor((ia-1)*(idim+1)+3)-w(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 q(1)=(xcoor4((ia-1)*(idim+1)+1)) q(2)=(xcoor4((ia-1)*(idim+1)+2)) q(3)=(xcoor4((ia-1)*(idim+1)+3)) xcoor2((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1) xcoor2((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2) xcoor2((ia-1)*(idim+1)+3)=zchoix 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 If (elecri(i1).eq.0) then 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 Endif 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=(XTOMA-XTOMI)/XZA + 1 NYZO=(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=(XTOMA-XTOMI)/XZO +1 NYZO=(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 SEGINI ISEG3 SEGINI ISEG4 DO 3 I1=1,NBEL If (elecri(i1).eq.0) then 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 201 L2=NIZ1X,NIZ2X NIZA = L2 + ( L1-1) * NXZO NUMZO(NIZA) = NUMZO(NIZA) +1 201 CONTINUE Endif 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 If (elecri(i1).eq.0) then 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 203 L2=NIZ1X,NIZ2X NIZA = L2 + ( L1-1) * NXZO IAD=NIZO(NIZA)+IDEJ(NIZA) NNMEL(IAD)=I1 IDEJ(NIZA)=IDEJ(NIZA)+1 203 CONTINUE Endif 5 CONTINUE C *********************** FIN DU ZONAGE DE MAI1 ************************** prec = 1.E-5 prec2 = -prec prec3 = -5.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 (tabcrit(1).eq.0) then 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 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 If ((k2.gt.0.) .and. (k2.lt.(nxzo+1)) .and. $ (k1.gt.0.) .and. (k1.lt.(nyzo+1))) then 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 Endif C PASSAGE TRANCHE CRITIQUE POUR UN POINT DU ZONAGE If (k2.lt.2 .or. k2.gt.(nxzo-1)) then Do 1462 iz=1,tabcrit(1) i1=tabcrit(iz+1) m1=num(1,i1) m2=num(2,i1) m3=num(3,i1) C CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2) xm1=xcoor3((m1-1)*(idim+1)+1) ym1=xcoor3((m1-1)*(idim+1)+2) zm1=xcoor3((m1-1)*(idim+1)+3) xm2=xcoor3((m2-1)*(idim+1)+1) ym2=xcoor3((m2-1)*(idim+1)+2) zm2=xcoor3((m2-1)*(idim+1)+3) xm3=xcoor3((m3-1)*(idim+1)+1) ym3=xcoor3((m3-1)*(idim+1)+2) zm3=xcoor3((m3-1)*(idim+1)+3) q(1)=(xcoor4((ia-1)*(idim+1)+1)) q(2)=(xcoor4((ia-1)*(idim+1)+2)) q(3)=(xcoor4((ia-1)*(idim+1)+3)) 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) zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3) xnew=q(1)*zk ynew=q(2)*zk znew=q(3)*zk C PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1 r1=sqrt(o1*o1+o2*o2+o3*o3) if( r1 . EQ. 0. ) then call erreur ( 21 ) endif o1 = o1/r1 o2 = o2/r1 o3 = o3/r1 r2=sqrt(o1*o1+o2*o2) if ( r2 . gt. 1.e-5) then tabmat(1,1)=-o2/r2 tabmat(1,2)=o1/r2 tabmat(1,3)=0. else tabmat(1,1)=1. tabmat(1,2)=0. tabmat(1,3)=0. endif tabmat(2,1)=-o3*tabmat(1,2) tabmat(2,2)=o3*tabmat(1,1) tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1) tabmat(3,1)=o1 tabmat(3,2)=o2 tabmat(3,3)=o3 x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3))) y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3))) x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3))) y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3))) x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3))) y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3))) x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3))) y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3))) C CALCUL DES COORDONNEES BARYCENTRIQUES 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 798 k=1,nbnn al(k)=(a(k)+b(k)*x+c(k)*y)/detam 798 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 1462 continue Endif If (icor(1).eq.0) then interr(1) = ia call erreur(782) return Endif Endif C FIN DU IF ICPR2(IA).EQ.1 If (icpr2(ia).eq.3) then C PASSAGE SUR LES ELEMENTS DE LA TRANCHE CRITIQUE Do 1463 iz=1,tabcrit(1) i1=tabcrit(iz+1) m1=num(1,i1) m2=num(2,i1) m3=num(3,i1) C CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2) xm1=xcoor3((m1-1)*(idim+1)+1) ym1=xcoor3((m1-1)*(idim+1)+2) zm1=xcoor3((m1-1)*(idim+1)+3) xm2=xcoor3((m2-1)*(idim+1)+1) ym2=xcoor3((m2-1)*(idim+1)+2) zm2=xcoor3((m2-1)*(idim+1)+3) xm3=xcoor3((m3-1)*(idim+1)+1) ym3=xcoor3((m3-1)*(idim+1)+2) zm3=xcoor3((m3-1)*(idim+1)+3) q(1)=(xcoor4((ia-1)*(idim+1)+1)) q(2)=(xcoor4((ia-1)*(idim+1)+2)) q(3)=(xcoor4((ia-1)*(idim+1)+3)) 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) zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3) xnew=q(1)*zk ynew=q(2)*zk znew=q(3)*zk C PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1 r1=sqrt(o1*o1+o2*o2+o3*o3) if( r1 . EQ. 0. ) then call erreur ( 21 ) endif o1 = o1/r1 o2 = o2/r1 o3 = o3/r1 r2=sqrt(o1*o1+o2*o2) if ( r2 . gt. 1.e-5) then tabmat(1,1)=-o2/r2 tabmat(1,2)=o1/r2 tabmat(1,3)=0. else tabmat(1,1)=1. tabmat(1,2)=0. tabmat(1,3)=0. endif tabmat(2,1)=-o3*tabmat(1,2) tabmat(2,2)=o3*tabmat(1,1) tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1) tabmat(3,1)=o1 tabmat(3,2)=o2 tabmat(3,3)=o3 x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3))) y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3))) x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3))) y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3))) x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3))) y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3))) x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3))) y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3))) C CALCUL DES COORDONNEES BARYCENTRIQUES 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 799 k=1,nbnn al(k)=(a(k)+b(k)*x+c(k)*y)/detam 799 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 1463 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) .or. (icpr2(ia).eq.3)) then dmin=1.e30 nbpts=nbpts+1 segadj mcoord n=0 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=xcoor3((m1-1)*(idim+1)+1) ym1=xcoor3((m1-1)*(idim+1)+2) zm1=xcoor3((m1-1)*(idim+1)+3) xm2=xcoor3((m2-1)*(idim+1)+1) ym2=xcoor3((m2-1)*(idim+1)+2) zm2=xcoor3((m2-1)*(idim+1)+3) xm3=xcoor3((m3-1)*(idim+1)+1) ym3=xcoor3((m3-1)*(idim+1)+2) zm3=xcoor3((m3-1)*(idim+1)+3) q(1)=(xcoor4((ia-1)*(idim+1)+1)) q(2)=(xcoor4((ia-1)*(idim+1)+2)) q(3)=(xcoor4((ia-1)*(idim+1)+3)) 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) zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3) xnew=q(1)*zk ynew=q(2)*zk znew=q(3)*zk delta1=sqrt(xnew*xnew+ynew*ynew+znew*znew) delta2=sqrt(q(1)*q(1)+q(2)* q(2)+q(3)*q(3)) pds=(xnew*q(1))+(ynew*q(2))+(znew*q(3)) If (pds.gt.0.) then If ((delta1-delta2).gt.0.) then If ((delta1-delta2).lt.dmin) then icorre(ia) = nbpts dmin=(delta1-delta2) xcoor((nbpts-1)*(idim+1)+1)=xnew xcoor((nbpts-1)*(idim+1)+2)=ynew xcoor((nbpts-1)*(idim+1)+3)=znew Endif n=n+1 Endif Endif 500 continue If (n.eq.0) then call erreur(21) return Endif 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))+w(J) 3200 continue icor(1)=0 n=0 Endif 405 continue segsup iseg1,iseg3,iseg4,iseg5,iseg6,mcoora,mcoor2 segsup icp1,icp2,mcoor3,mcoor4 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 segsup mcor,mcorre,matric Return End