Changeset 2841


Ignore:
Timestamp:
01/14/10 15:35:28 (15 years ago)
Author:
Mathieu Morlighem
Message:

added Bamgx/objects/QuadTree.cpp

Location:
issm/trunk/src/c
Files:
1 added
1 deleted
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Bamgx/QuadTree.h

    r2800 r2841  
    1 // -*- Mode : c++ -*-
    2 //
    3 // SUMMARY  :     
    4 // USAGE    :       
    5 // ORG      :
    6 // AUTHOR   : Frederic Hecht
    7 // E-MAIL   : hecht@ann.jussieu.fr
    8 //
    9 
    10 /*
    11  
    12  This file is part of Freefem++
    13  
    14  Freefem++ is free software; you can redistribute it and/or modify
    15  it under the terms of the GNU Lesser General Public License as published by
    16  the Free Software Foundation; either version 2.1 of the License, or
    17  (at your option) any later version.
    18  
    19  Freefem++  is distributed in the hope that it will be useful,
    20  but WITHOUT ANY WARRANTY; without even the implied warranty of
    21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    22  GNU Lesser General Public License for more details.
    23  
    24  You should have received a copy of the GNU Lesser General Public License
    25  along with Freefem++; if not, write to the Free Software
    26  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
    27  */
    28 
    291namespace bamg {
    302
    … …  
    10678};
    10779}
    108 //#undef IJ

  • issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp

    r2820 r2841  
    1515        /*Constructor/Destructor*/
    1616
    17         /*Others*/
     17        /*Methods*/
    1818        /*FUNCTION GeometricalEdge::R1tg{{{1*/
    1919        Real8 GeometricalEdge::R1tg(Real8 theta,R2 & t) const // 1/R of radius of cuvature

  • issm/trunk/src/c/Bamgx/objects/Geometry.cpp

    r2820 r2841  
    1717
    1818        /*Constructors/Destructors*/
     19        /*FUNCTION  Geometry::Geometry(const Geometry & Gh){{{1*/
     20        Geometry::Geometry(const Geometry & Gh) {
     21                Int4 i;
     22                *this = Gh;
     23                NbRef =0;
     24                quadtree=0;
     25                name = new char[strlen(Gh.name)+4];
     26                strcpy(name,"cp:");
     27                strcat(name,Gh.name);
     28                vertices = nbv ? new GeometricalVertex[nbv] : NULL;
     29                triangles = nbt ? new  Triangle[nbt]:NULL;
     30                edges = nbe ? new GeometricalEdge[nbe]:NULL;
     31                curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;
     32                subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL;
     33                for (i=0;i<nbv;i++)
     34                 vertices[i].Set(Gh.vertices[i],Gh,*this);
     35                for (i=0;i<nbe;i++)
     36                 edges[i].Set(Gh.edges[i],Gh,*this);
     37                for (i=0;i<NbOfCurves;i++)
     38                 curves[i].Set(Gh.curves[i],Gh,*this);
     39                for (i=0;i<NbSubDomains;i++)
     40                 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
     41
     42                //    for (i=0;i<nbt;i++)
     43                //      triangles[i].Set(Gh.triangles[i],Gh,*this);
     44                assert(!nbt);   
     45        }
     46        /*}}}1*/
     47        /*FUNCTION Geometry::~Geometry(){{{1*/
     48        Geometry::~Geometry() {
     49                long int verbosity=0;
     50
     51                assert(NbRef<=0);
     52                if(verbosity>9)
     53                 cout << "DELETE      ~Geometry "<< this  << endl;
     54                if(vertices)  delete [] vertices;vertices=0;
     55                if(edges)     delete [] edges;edges=0;
     56                // if(edgescomponante) delete [] edgescomponante; edgescomponante=0;
     57                if(triangles) delete [] triangles;triangles=0;
     58                if(quadtree)  delete  quadtree;quadtree=0;
     59                if(curves)  delete  []curves;curves=0;NbOfCurves=0;
     60                if(name) delete [] name;name=0;
     61                if(subdomains) delete [] subdomains;subdomains=0;
     62                //  if(ordre)     delete [] ordre;
     63                EmptyGeometry();
     64        }
     65        /*}}}1*/
    1966
    2067        /*IO*/
    … …  
    433480        /*}}}1*/
    434481
    435         /*Other*/
    436 
    437 void Geometry::EmptyGeometry()  // empty geometry
    438   {
    439   OnDisk=0;
    440   NbRef=0;
    441   name =0;
    442   quadtree=0;
    443   curves=0;
    444  // edgescomponante=0;
    445   triangles=0;
    446   edges=0;
    447   vertices=0;
    448   NbSubDomains=0;
    449   //  nbtf=0;
    450 //  BeginOfCurve=0; 
    451   nbiv=nbv=nbvx=0;
    452   nbe=nbt=nbtx=0;
    453   NbOfCurves=0;
    454 //  BeginOfCurve=0;
    455   subdomains=0;
    456   MaximalAngleOfCorner = 10*Pi/180;
    457   }
    458 
    459 
    460 
    461 Geometry::Geometry(const Geometry & Gh)
    462  { Int4 i;
    463    *this = Gh;
    464    NbRef =0;
    465    quadtree=0;
    466    name = new char[strlen(Gh.name)+4];
    467    strcpy(name,"cp:");
    468    strcat(name,Gh.name);
    469    vertices = nbv ? new GeometricalVertex[nbv] : NULL;
    470    triangles = nbt ? new  Triangle[nbt]:NULL;
    471    edges = nbe ? new GeometricalEdge[nbe]:NULL;
    472    curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;
    473    subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL;
    474    for (i=0;i<nbv;i++)
    475      vertices[i].Set(Gh.vertices[i],Gh,*this);
    476    for (i=0;i<nbe;i++)
    477      edges[i].Set(Gh.edges[i],Gh,*this);
    478    for (i=0;i<NbOfCurves;i++)
    479      curves[i].Set(Gh.curves[i],Gh,*this);
    480    for (i=0;i<NbSubDomains;i++)
    481      subdomains[i].Set(Gh.subdomains[i],Gh,*this);
    482      
    483    //    for (i=0;i<nbt;i++)
    484    //      triangles[i].Set(Gh.triangles[i],Gh,*this);
    485    assert(!nbt);   
    486  }
    487 
    488 
    489 GeometricalEdge* Geometry::Contening(const R2 P,  GeometricalEdge * start) const
    490 {
    491   GeometricalEdge* on =start,* pon=0;
    492   // walk with the cos on geometry
    493   //  cout << P ;
    494   int k=0;
    495    while(pon != on)
    496      { 
    497        pon = on;
    498        assert(k++<100);
    499        R2 A= (*on)[0];
    500        R2 B= (*on)[1];
    501        R2 AB = B-A;
    502        R2 AP = P-A;
    503        R2 BP = P-B;
    504        //   cout << "::  " << on - edges << " "  <<  AB*AP  << " " <<  AB*BP << " " << A << B << endl;
    505        if ( (AB,AP) < 0)
    506          on = on->Adj[0];
    507        else if ( (AB,BP)  > 0)
    508          on = on->Adj[1];
    509        else
    510          return on;
    511      }
    512    return on;
    513 }
    514 GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const
    515  { 
    516     Real8 save_s=s;
    517     int NbTry=0;
     482        /*Methods*/
     483        /*FUNCTION  Geometry::EmptyGeometry(){{{1*/
     484        void Geometry::EmptyGeometry() {
     485                OnDisk=0;
     486                NbRef=0;
     487                name =0;
     488                quadtree=0;
     489                curves=0;
     490                // edgescomponante=0;
     491                triangles=0;
     492                edges=0;
     493                vertices=0;
     494                NbSubDomains=0;
     495                //  nbtf=0;
     496                //  BeginOfCurve=0; 
     497                nbiv=nbv=nbvx=0;
     498                nbe=nbt=nbtx=0;
     499                NbOfCurves=0;
     500                //  BeginOfCurve=0;
     501                subdomains=0;
     502                MaximalAngleOfCorner = 10*Pi/180;
     503        }
     504        /*}}}1*/
     505        /*FUNCTION  Geometry::Contening{{{1*/
     506        GeometricalEdge* Geometry::Contening(const R2 P,  GeometricalEdge * start) const {
     507                GeometricalEdge* on =start,* pon=0;
     508                // walk with the cos on geometry
     509                //  cout << P ;
     510                int k=0;
     511                while(pon != on)
     512                  { 
     513                        pon = on;
     514                        assert(k++<100);
     515                        R2 A= (*on)[0];
     516                        R2 B= (*on)[1];
     517                        R2 AB = B-A;
     518                        R2 AP = P-A;
     519                        R2 BP = P-B;
     520                        //   cout << "::  " << on - edges << " "  <<  AB*AP  << " " <<  AB*BP << " " << A << B << endl;
     521                        if ( (AB,AP) < 0)
     522                         on = on->Adj[0];
     523                        else if ( (AB,BP)  > 0)
     524                         on = on->Adj[1];
     525                        else
     526                         return on;
     527                  }
     528                return on;
     529        }
     530        /*}}}1*/
     531        /*FUNCTION  Geometry::Geometry::ProjectOnCurve {{{1*/
     532        GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const {
     533                Real8 save_s=s;
     534                int NbTry=0;
    518535retry:   
    519     s=save_s;
    520     GeometricalEdge * on = e.on;
    521     assert(on);
    522     assert( e[0].on &&  e[1].on);
    523     const Vertex &v0=e[0],&v1=e[1];
    524     V.m = Metric(1.0-s, v0,s, v1);
     536                s=save_s;
     537                GeometricalEdge * on = e.on;
     538                assert(on);
     539                assert( e[0].on &&  e[1].on);
     540                const Vertex &v0=e[0],&v1=e[1];
     541                V.m = Metric(1.0-s, v0,s, v1);
    525542#define MXE__LINE  __LINE__+1
    526     const int mxe =100;
    527     GeometricalEdge *ge[mxe+1];
    528     int    sensge[mxe+1];
    529     Real8  lge[mxe+1];
    530     int bge=mxe/2,tge=bge;
    531     ge[bge] = e.on;
    532     sensge[bge]=1;
    533 
    534     R2 V0 = v0,V1=v1,V01=V1-V0;
    535     VertexOnGeom  vg0= *v0.on,  vg1=*v1.on;
    536     if(NbTry) cout << "bug: s==== " << s << " e=" <<  V0 << " " << V1 << endl;
    537 
    538     //    GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL;
    539     GeometricalEdge * eg0=on, *eg1=on;
    540     R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
    541     if(NbTry) cout <<" G edge= " << Ag << Bg << endl << " v edge" << V01 << " v geom " << AB  <<  (V01,AB) <<endl;
    542     int OppositeSens = (V01,AB) < 0;
    543     int sens0=0,sens1=1;
    544     if (OppositeSens)
    545       s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
    546     if(NbTry) cout << "bug: edge = " << v0.r << " -> " << v1.r << endl
    547                    << "sg 0 = " << vg0
    548                    << " on = " << Number(on) << ":" << Ag << Bg << "; "
    549                    <<  " sg 1= " << vg1
    550                    << "--------------------------------------------" << endl;
    551     while (eg0 != (GeometricalEdge*) vg0  &&  (*eg0)(sens0) != (GeometricalVertex*) vg0)
    552       {
    553       if (bge<=0) {
    554         //          int kkk;
    555           // if (NbTry) cout <<"Read (int) to Show Sioux window", cin >> kkk ;
    556                if(NbTry)
    557                   {
    558                     cerr << " -- Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl;
    559                     cerr << "   The mesh of the  Geometry is to fine: ";
    560                     cerr << "     1)  a mesh edge  contening more than "<< mxe/2 << " geometrical edges." << endl;
    561                     cerr << "     2)  code bug : be sure that we call   Triangles::SetVertexFieldOn() before " << endl;
    562                     cerr << "   To solve the problem do a coarsening of the geometrical mesh " << endl;
    563                     cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl;       
    564                     MeshError(222);
     543                const int mxe =100;
     544                GeometricalEdge *ge[mxe+1];
     545                int    sensge[mxe+1];
     546                Real8  lge[mxe+1];
     547                int bge=mxe/2,tge=bge;
     548                ge[bge] = e.on;
     549                sensge[bge]=1;
     550
     551                R2 V0 = v0,V1=v1,V01=V1-V0;
     552                VertexOnGeom  vg0= *v0.on,  vg1=*v1.on;
     553                if(NbTry) cout << "bug: s==== " << s << " e=" <<  V0 << " " << V1 << endl;
     554
     555                //    GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL;
     556                GeometricalEdge * eg0=on, *eg1=on;
     557                R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
     558                if(NbTry) cout <<" G edge= " << Ag << Bg << endl << " v edge" << V01 << " v geom " << AB  <<  (V01,AB) <<endl;
     559                int OppositeSens = (V01,AB) < 0;
     560                int sens0=0,sens1=1;
     561                if (OppositeSens)
     562                 s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
     563                if(NbTry) cout << "bug: edge = " << v0.r << " -> " << v1.r << endl
     564                 << "sg 0 = " << vg0
     565                        << " on = " << Number(on) << ":" << Ag << Bg << "; "
     566                        <<  " sg 1= " << vg1
     567                        << "--------------------------------------------" << endl;
     568                while (eg0 != (GeometricalEdge*) vg0  &&  (*eg0)(sens0) != (GeometricalVertex*) vg0)
     569                  {
     570                        if (bge<=0) {
     571                                //          int kkk;
     572                                // if (NbTry) cout <<"Read (int) to Show Sioux window", cin >> kkk ;
     573                                if(NbTry)
     574                                  {
     575                                        cerr << " -- Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl;
     576                                        cerr << "   The mesh of the  Geometry is to fine: ";
     577                                        cerr << "     1)  a mesh edge  contening more than "<< mxe/2 << " geometrical edges." << endl;
     578                                        cerr << "     2)  code bug : be sure that we call   Triangles::SetVertexFieldOn() before " << endl;
     579                                        cerr << "   To solve the problem do a coarsening of the geometrical mesh " << endl;
     580                                        cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl;   
     581                                        MeshError(222);
     582                                  }
     583                                NbTry++;
     584                                goto retry;}
     585                                GeometricalEdge* tmpge = eg0;
     586                                if(NbTry)
     587                                 cout << "bug: --Edge @" <<  Number(tmpge)  << " = "<< Number(eg0) << ":" <<Number(eg0->Adj[0]) << "," << 
     588                                        Number(eg0->Adj[1]) <<"," ;
     589                                ge[--bge] =eg0 = eg0->Adj[sens0];
     590                                assert(bge>=0 && bge <= mxe);
     591                                sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]);
     592                                if(NbTry)
     593                                 cout << "bug: Edge "  <<  Number(eg0) << " "<< 1-sens0 <<  " S "
     594                                        << Number((*eg0)[1-sens0]) <<":" << Number(eg0->Adj[0]) << ","
     595                                        <<  Number(eg0->Adj[1]) <<"," << endl
     596                                        <<Number(eg0)<< (*eg0)[sens0].r << "v = " << Number((*eg1)(sens0)) << " e = " << eg0 <<  endl;
    565597                  }
    566             NbTry++;
    567             goto retry;}
    568         GeometricalEdge* tmpge = eg0;
    569          if(NbTry)
    570                 cout << "bug: --Edge @" <<  Number(tmpge)  << " = "<< Number(eg0) << ":" <<Number(eg0->Adj[0]) << "," << 
    571                          Number(eg0->Adj[1]) <<"," ;
    572         ge[--bge] =eg0 = eg0->Adj[sens0];
    573         assert(bge>=0 && bge <= mxe);
    574         sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]);
    575         if(NbTry)
    576                 cout << "bug: Edge "  <<  Number(eg0) << " "<< 1-sens0 <<  " S "
    577                      << Number((*eg0)[1-sens0]) <<":" << Number(eg0->Adj[0]) << ","
    578                      <<  Number(eg0->Adj[1]) <<"," << endl
    579                      <<Number(eg0)<< (*eg0)[sens0].r << "v = " << Number((*eg1)(sens0)) << " e = " << eg0 <<  endl;
    580      }
    581       if(NbTry) cout << Number((GeometricalEdge*) vg1) << " " << Number((GeometricalVertex*) vg1) << endl;
    582     while (eg1 != (GeometricalEdge*) vg1  &&  (*eg1)(sens1) != (GeometricalVertex*) vg1)
    583       {
    584         if(tge>=mxe ) {
    585           cerr << " --Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl;
    586           NbTry++;
    587           if (NbTry<2) goto retry;
    588           cerr << "   The mesh of the  Geometry is to fine:" ;
    589           cerr << "     1)  a mesh edge  contening more than "<< mxe/2 << " geometrical edges." << endl;
    590           cerr << "     2)  code bug : be sure that we call   Triangles::SetVertexFieldOn() before " << endl;
    591           cerr << "   To solve the problem do a coarsening of the geometrical mesh " << endl;
    592           cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl;       
    593           MeshError(223);
     598                if(NbTry) cout << Number((GeometricalEdge*) vg1) << " " << Number((GeometricalVertex*) vg1) << endl;
     599                while (eg1 != (GeometricalEdge*) vg1  &&  (*eg1)(sens1) != (GeometricalVertex*) vg1)
     600                  {
     601                        if(tge>=mxe ) {
     602                                cerr << " --Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl;
     603                                NbTry++;
     604                                if (NbTry<2) goto retry;
     605                                cerr << "   The mesh of the  Geometry is to fine:" ;
     606                                cerr << "     1)  a mesh edge  contening more than "<< mxe/2 << " geometrical edges." << endl;
     607                                cerr << "     2)  code bug : be sure that we call   Triangles::SetVertexFieldOn() before " << endl;
     608                                cerr << "   To solve the problem do a coarsening of the geometrical mesh " << endl;
     609                                cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl;
     610                                MeshError(223);
     611                        }
     612
     613                        GeometricalEdge* tmpge = eg1;
     614                        if(NbTry)
     615                         cout << "++Edge @" << tmpge << " = " <<  Number(eg1) <<"%" << Number(eg1->Adj[0]) << ","
     616                                <<  Number(eg1->Adj[1]) <<"," ;
     617                        ge[++tge] =eg1 = eg1->Adj[sens1];
     618                        sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1];
     619                        assert(tge>=0 && tge <= mxe);
     620                        if(NbTry)
     621                         cout << "  Edge "  <<  Number(eg1) << " " << sens1 << " S "
     622                                <<Number((*eg1)[sens1]) <<"%"<< Number(eg1->Adj[0]) << "," <<  Number(eg1->Adj[1]) <<","
     623                                <<Number(eg1)<< (*eg1)[sens1].r << "v = " << Number((*eg1)(sens1)) << " e = " << Number(eg1) <<  endl;
     624                  }
     625
     626                if(NbTry)    cout << endl;
     627
     628
     629                if ( (*eg0)(sens0) == (GeometricalVertex*) vg0 )
     630                 vg0 = VertexOnGeom( *(Vertex *) vg0,*eg0,sens0);
     631
     632                if ( (*eg1)(sens1) == (GeometricalVertex*) vg1)
     633                 vg1 = VertexOnGeom( *(Vertex *) vg1,*eg1,sens1);
     634
     635                Real8 sg;
     636                //   cout << "           " << Number(on) << " " <<  Number(eg0) << " " <<  Number(eg1) << " "  ;
     637                if (eg0 == eg1) {
     638                        register Real8 s0= vg0,s1=vg1;
     639                        sg =  s0 * (1.0-s) +  s * s1;
     640                        //    cout <<"                s0=" << s0 << " s1=" << s1
     641                        //             << " s = " << s << " sens= " << OppositeSens << "\t\t sg = " << sg << endl ;
     642                        on=eg0;}
     643                else {
     644                        R2 AA=V0,BB;
     645                        Real8 s0,s1;
     646
     647                        //cout << endl << "s= " << s << Number(eg0) << " " << (Real8) vg0 << " "
     648                        //          << Number(eg1) << " " << (Real8) vg1 << V0 << V1 << "  Interpol = "
     649                        // << V0*(1-s)+V1*s << ";;; " <<  endl;
     650                        int i=bge;
     651                        Real8 ll=0;
     652                        for(i=bge;i<tge;i++)
     653                          {
     654                                assert( i>=0 && i <= mxe);
     655                                BB =  (*ge[i])[sensge[i]];
     656                                lge[i]=ll += Norme2(AA-BB);
     657                                //   cout << " ll " << i << BB << ll << " " <<sensge[i] <<" on = " <<
     658                                // Number(ge[i]) << " sens= " << sensge[i] ;
     659                                AA=BB ;}
     660                                lge[tge]=ll+=Norme2(AA-V1);
     661                                // cout << " ll " << tge << " " << ll <<  sensge[tge]
     662                                //           <<" on = " << Number(ge[tge]) <<  " sens= " << sensge[tge] << endl;
     663                                // search the geometrical edge
     664                                assert(s <= 1.0);
     665                                Real8 ls= s*ll;
     666                                on =0;
     667                                s0 = vg0;
     668                                s1= sensge[bge];
     669                                Real8 l0=0,l1;
     670                                i=bge;
     671                                while (  (l1=lge[i]) < ls ) {
     672                                        assert(i >= 0 && i <= mxe);
     673                                        i++,s0=1-(s1=sensge[i]),l0=l1;}
     674                                        on=ge[i];
     675                                        if (i==tge)
     676                                         s1=vg1;
     677
     678                                        s=(ls-l0)/(l1-l0);
     679                                        //  cout << "on =" << Number(on) << sens0 << sens1 <<  "s0  " << s0 << " s1 ="
     680                                        //           << s1 << " l0 =" << l0 << " ls= " << ls << " l1= " << l1 << " s= " << s;
     681                                        sg =  s0 * (1.0-s) +  s * s1;   
     682                }
     683                assert(on);
     684                // assert(sg && sg-1);
     685                V.r= on->F(sg);
     686                //  if (eg0 != eg1)
     687                //        cout << "----- sg = "<< sg << " Sens =" << OppositeSens << " Edge = "
     688                //     << Number(on) <<"  V=" << V << endl;
     689                GV=VertexOnGeom(V,*on,sg);
     690                return on;
    594691        }
    595 
    596         GeometricalEdge* tmpge = eg1;
    597         if(NbTry)
    598                 cout << "++Edge @" << tmpge << " = " <<  Number(eg1) <<"%" << Number(eg1->Adj[0]) << ","
    599                  <<  Number(eg1->Adj[1]) <<"," ;
    600         ge[++tge] =eg1 = eg1->Adj[sens1];
    601         sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1];
    602         assert(tge>=0 && tge <= mxe);
    603          if(NbTry)
    604                 cout << "  Edge "  <<  Number(eg1) << " " << sens1 << " S "
    605                      <<Number((*eg1)[sens1]) <<"%"<< Number(eg1->Adj[0]) << "," <<  Number(eg1->Adj[1]) <<","
    606                      <<Number(eg1)<< (*eg1)[sens1].r << "v = " << Number((*eg1)(sens1)) << " e = " << Number(eg1) <<  endl;
    607       }
    608 
    609         if(NbTry)    cout << endl;
    610        
    611 
    612     if ( (*eg0)(sens0) == (GeometricalVertex*) vg0 )
    613       vg0 = VertexOnGeom( *(Vertex *) vg0,*eg0,sens0);
    614    
    615     if ( (*eg1)(sens1) == (GeometricalVertex*) vg1)
    616        vg1 = VertexOnGeom( *(Vertex *) vg1,*eg1,sens1);
    617 
    618     Real8 sg;
    619     //   cout << "           " << Number(on) << " " <<  Number(eg0) << " " <<  Number(eg1) << " "  ;
    620     if (eg0 == eg1) {
    621        register Real8 s0= vg0,s1=vg1;
    622        sg =  s0 * (1.0-s) +  s * s1;
    623        //    cout <<"                s0=" << s0 << " s1=" << s1
    624        //             << " s = " << s << " sens= " << OppositeSens << "\t\t sg = " << sg << endl ;
    625        on=eg0;}
    626     else {
    627        R2 AA=V0,BB;
    628        Real8 s0,s1;
    629        
    630        //cout << endl << "s= " << s << Number(eg0) << " " << (Real8) vg0 << " "
    631        //           << Number(eg1) << " " << (Real8) vg1 << V0 << V1 << "  Interpol = "
    632        // << V0*(1-s)+V1*s << ";;; " <<  endl;
    633        int i=bge;
    634        Real8 ll=0;
    635        for(i=bge;i<tge;i++)
    636          {
    637            assert( i>=0 && i <= mxe);
    638            BB =  (*ge[i])[sensge[i]];
    639            lge[i]=ll += Norme2(AA-BB);
    640            //   cout << " ll " << i << BB << ll << " " <<sensge[i] <<" on = " <<
    641            // Number(ge[i]) << " sens= " << sensge[i] ;
    642            AA=BB ;}
    643        lge[tge]=ll+=Norme2(AA-V1);
    644        // cout << " ll " << tge << " " << ll <<  sensge[tge]
    645        //            <<" on = " << Number(ge[tge]) <<  " sens= " << sensge[tge] << endl;
    646     // search the geometrical edge
    647       assert(s <= 1.0);
    648       Real8 ls= s*ll;
    649       on =0;
    650       s0 = vg0;
    651       s1= sensge[bge];
    652       Real8 l0=0,l1;
    653       i=bge;
    654       while (  (l1=lge[i]) < ls ) {
    655         assert(i >= 0 && i <= mxe);
    656         i++,s0=1-(s1=sensge[i]),l0=l1;}
    657       on=ge[i];
    658       if (i==tge)
    659         s1=vg1;
    660      
    661       s=(ls-l0)/(l1-l0);
    662       //  cout << "on =" << Number(on) << sens0 << sens1 <<  "s0  " << s0 << " s1 ="
    663       //             << s1 << " l0 =" << l0 << " ls= " << ls << " l1= " << l1 << " s= " << s;
    664        sg =  s0 * (1.0-s) +  s * s1;   
    665        }
    666     assert(on);
    667     // assert(sg && sg-1);
    668     V.r= on->F(sg);
    669     //  if (eg0 != eg1)
    670     //        cout << "----- sg = "<< sg << " Sens =" << OppositeSens << " Edge = "
    671     //     << Number(on) <<"  V=" << V << endl;
    672     GV=VertexOnGeom(V,*on,sg);
    673     return on;
    674  }
    675 
    676 void Geometry::AfterRead()
    677  {// -----------------
    678          long int verbosity=0;
    679 
    680     if (verbosity>20)
    681     cout << "Geometry::AfterRead()" <<  nbv << " " << nbe << endl;
    682     Int4 i,k=0;        ;
    683     int jj; // jj in [0,1]
    684     Int4 * hv = new Int4 [ nbv];
    685     Int4 * ev = new Int4 [ 2 * nbe ];
    686     float  * eangle = new float[ nbe ];
    687     {
    688       double eps = 1e-20;
    689       QuadTree quadtree; // to find same vertices
    690       Vertex * v0 = vertices;
    691       GeometricalVertex  * v0g = (GeometricalVertex  *) (void *) v0;   
    692       int k=0;
    693       for (i=0;i<nbv;i++)
    694         vertices[i].link = vertices +i;
    695       for (i=0;i<nbv;i++)
    696              {
    697               vertices[i].i = toI2(vertices[i].r); // set integer coordinate
    698               Vertex *v= quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
    699               if( v && Norme1(v->r - vertices[i]) < eps )
    700                { // link v & vertices[i]
    701                  // vieille ruse pour recuperer j
    702                  GeometricalVertex * vg = (GeometricalVertex  *) (void *) v;
    703                  int j = vg-v0g;
    704                  assert( v ==  & (Vertex &) vertices[j]);
    705                  vertices[i].link = vertices + j;
    706             k++;             
    707                }
    708               else  quadtree.Add(vertices[i]);
    709              }
    710       if (k) {
    711         cout << " Number of distinte vertices " << nbv - k << " Over " << nbv << endl;
    712         //if (verbosity>10)
    713         {
    714           cout << " The duplicate vertex " << endl;
    715           for (i=0;i<nbv;i++)
    716             if (!vertices[i].IsThe())
    717               cout << " " << i << " and " << Number(vertices[i].The()) << endl;
    718           MeshError(102);
    719           //throw(ErrorExec("exit",1));   
     692        /*}}}1*/
     693        /*FUNCTION  Geometry::AfterRead(){{{1*/
     694        void Geometry::AfterRead() {
     695                long int verbosity=0;
     696
     697                if (verbosity>20)
     698                 cout << "Geometry::AfterRead()" <<  nbv << " " << nbe << endl;
     699                Int4 i,k=0;        ;
     700                int jj; // jj in [0,1]
     701                Int4 * hv = new Int4 [ nbv];
     702                Int4 * ev = new Int4 [ 2 * nbe ];
     703                float  * eangle = new float[ nbe ];
     704                  {
     705                        double eps = 1e-20;
     706                        QuadTree quadtree; // to find same vertices
     707                        Vertex * v0 = vertices;
     708                        GeometricalVertex  * v0g = (GeometricalVertex  *) (void *) v0;   
     709                        int k=0;
     710                        for (i=0;i<nbv;i++)
     711                         vertices[i].link = vertices +i;
     712                        for (i=0;i<nbv;i++)
     713                          {
     714                                vertices[i].i = toI2(vertices[i].r); // set integer coordinate
     715                                Vertex *v= quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
     716                                if( v && Norme1(v->r - vertices[i]) < eps )
     717                                  { // link v & vertices[i]
     718                                        // vieille ruse pour recuperer j
     719                                        GeometricalVertex * vg = (GeometricalVertex  *) (void *) v;
     720                                        int j = vg-v0g;
     721                                        assert( v ==  & (Vertex &) vertices[j]);
     722                                        vertices[i].link = vertices + j;
     723                                        k++;         
     724                                  }
     725                                else  quadtree.Add(vertices[i]);
     726                          }
     727                        if (k) {
     728                                cout << " Number of distinte vertices " << nbv - k << " Over " << nbv << endl;
     729                                //if (verbosity>10)
     730                                  {
     731                                        cout << " The duplicate vertex " << endl;
     732                                        for (i=0;i<nbv;i++)
     733                                         if (!vertices[i].IsThe())
     734                                          cout << " " << i << " and " << Number(vertices[i].The()) << endl;
     735                                        MeshError(102);
     736                                        //throw(ErrorExec("exit",1));   
     737                                  }
     738                        }
     739
     740                        //  verification of cracked edge
     741                        int err =0;
     742                        for (i=0;i<nbe;i++)
     743                         if (edges[i].Cracked() )
     744                                {
     745                                 //    verification of crack
     746                                 GeometricalEdge & e1=edges[i];
     747                                 GeometricalEdge & e2=*e1.link;
     748                                 cerr << i << " " << e1[0].The() << " " << e2[0].The() << " " <<  e1[1].The() << " " << e2[1].The() << endl;
     749                                 if ( e1[0].The() == e2[0].The() && e1[1].The() == e2[1].The() )
     750                                        {
     751                                        }
     752                                 else
     753                                  if ( e1[0].The() == e2[1].The() && e1[1].The() == e2[0].The() )
     754                                         {
     755                                         }
     756                                  else
     757                                         {
     758                                          err++;
     759                                          cerr << " Cracked edges with no same vertex " << &e1-edges << " " << &e2 -edges << endl;
     760                                         }
     761                                }
     762                         else
     763                                {
     764                                 //  if (!edges[i][0].IsThe()) err++;
     765                                 // if (!edges[i][1].IsThe()) err++;
     766                                }
     767                        if (err)
     768                          {
     769                                cerr << " Some vertex was not distint and not on cracked edge " << err<< endl;
     770                                MeshError(222);
     771                          }
     772                  }
     773                if(verbosity>7)
     774                 for (i=0;i<nbv;i++)
     775                  if (vertices[i].Required())
     776                        cout << "     The geo vertices  " << i << " is required" << endl;
     777
     778                for (i=0;i<nbv;i++)
     779                 hv[i]=-1;// empty list
     780
     781                for (i=0;i<nbe;i++)
     782                  {
     783                        R2 v10  =  edges[i].v[1]->r -  edges[i].v[0]->r;
     784                        Real8 lv10 = Norme2(v10);
     785                        if(lv10 == 0) {
     786                                cerr << "The length  of " <<i<< "th Egde is 0 " << endl ;
     787                                MeshError(1);}
     788                                eangle[i] = atan2(v10.y,v10.x)  ; // angle in [ -Pi,Pi ]
     789                                if(verbosity>9)
     790                                 cout << "     angle edge " << i <<" " << eangle[i]*180/Pi<< v10<<endl;
     791                                for (jj=0;jj<2;jj++)
     792                                  { // generation of list
     793                                        Int4 v =  Number(edges[i].v[jj]);
     794                                        ev[k] = hv[v];
     795                                        hv[v] = k++;
     796                                  }
     797                  }
     798                // bulle sort on the angle of edge 
     799                for (i=0;i<nbv;i++) {
     800                        int exch = 1,ord =0;     
     801                        while (exch) {
     802                                exch = 0;
     803                                Int4  *p =  hv + i, *po = p;
     804                                Int4 n = *p;
     805                                register float angleold = -1000 ; // angle = - infini
     806                                ord = 0;
     807                                while (n >=0)
     808                                  {
     809                                        ord++;
     810                                        register Int4 i1= n /2;
     811                                        register Int4  j1 = n % 2;
     812                                        register Int4 *pn = ev + n;
     813                                        float angle = j1 ? OppositeAngle(eangle[i1]):  eangle[i1];
     814                                        n = *pn;
     815                                        if (angleold > angle) // exch to have : po -> pn -> p
     816                                         exch=1,*pn = *po,*po=*p,*p=n,po = pn;
     817                                        else //  to have : po -> p -> pn
     818                                         angleold =  angle, po = p,p  = pn;
     819                                  }
     820                        } // end while (exch)
     821
     822                        if (ord >= 1 )
     823                          { /*
     824                                         Int4 n = hv[i];
     825                                         while ( n >=0)
     826                                         { Int4 i1 = n/2,j1 = n%2;
     827                                //float a = 180*(j1 ? OppositeAngle(eangle[i1]): eangle[i1])/Pi;
     828                                n = ev[n];
     829                                }
     830                                */
     831                          }
     832                        if(ord == 2) { // angulare test to find a corner
     833                                Int4 n1 = hv[i];
     834                                Int4 n2 = ev[n1];
     835                                Int4 i1 = n1 /2, i2 = n2/2; // edge number
     836                                Int4  j1 = n1 %2, j2 = n2%2; // vertex in the edge
     837                                float angle1= j1 ? OppositeAngle(eangle[i1]) :  eangle[i1];
     838                                float angle2= !j2 ? OppositeAngle(eangle[i2]) :  eangle[i2];
     839                                float da12 = Abs(angle2-angle1);
     840                                if(verbosity>9)
     841                                 cout <<"     check angle " << i << " " << i1 << " " << i2  << " " << 180*(da12)/Pi
     842                                        << " " << 180*MaximalAngleOfCorner/Pi << vertices[i] << endl;
     843
     844                                if (( da12 >= MaximalAngleOfCorner )
     845                                                        && (da12 <= 2*Pi -MaximalAngleOfCorner)) {
     846                                        vertices[i].SetCorner() ;
     847                                        if(verbosity>7)
     848                                         cout << "     The vertex " << i << " is a corner (angle) "
     849                                                << 180*(da12)/ Pi<< " " << 180*MaximalAngleOfCorner/Pi << endl;}
     850                                        // if the ref a changing then is     SetRequired();
     851
     852                                        if (edges[i1].flag != edges[i2].flag || edges[i1].Required())
     853                                          {
     854                                                vertices[i].SetRequired();
     855                                                if(verbosity>7)
     856                                                 cout << "     The vertex " << i << " is Required the flag change (crack or equi, or require)" << endl;}
     857
     858                                                if (edges[i1].ref != edges[i2].ref) {
     859                                                        vertices[i].SetRequired();
     860                                                        if(verbosity>7)
     861                                                         cout << "     The vertex " << i << " is Required ref" << endl;}
     862                        } ;
     863
     864                        if(ord != 2) {
     865                                vertices[i].SetCorner();
     866                                if(verbosity>7)
     867                                 cout << "     the vertex " << i << " is a corner ordre = " << ord << endl;
     868                        }
     869                        // close the liste around the vertex
     870                          { Int4 no=-1, ne = hv[i];
     871                                while ( ne >=0)
     872                                 ne = ev[no=ne];       
     873                                if(no>=0)
     874                                 ev[no] = hv[i];
     875                          } // now the list around the vertex is circular
     876
     877                } // end for (i=0;i<nbv;i++)
     878
     879                k =0;
     880                for (i=0;i<nbe;i++)
     881                 for (jj=0;jj<2;jj++){
     882                         Int4 n1 = ev[k++];
     883                         Int4 i1 = n1/2 ,j1=n1%2;
     884                         if( edges[i1].v[j1] != edges[i].v[jj])
     885                                { cerr << " Bug Adj edge " << i << " " << jj <<
     886                                 " et " << i1 << " " << j1 << " k=" << k;
     887                                 cerr << Number(edges[i].v[jj]) <<" <> "
     888                                        << Number(edges[i1].v[j1])  <<endl;
     889                                 cerr << "edge " << Number(edges[i].v[0]) << " "
     890                                        << Number(edges[i].v[1]) << endl;
     891                                 //    cerr << "in file " <<filename <<endl;
     892                                 MeshError(1);
     893                                }
     894                         edges[i1].Adj[j1] = edges + i;
     895                         edges[i1].SensAdj[j1] = jj;
     896                         if (verbosity>10)
     897                          cout << " edges. Adj " << i1 << " " << j1 << " <--- " << i << " " << jj << endl;
     898                 }
     899
     900                // generation of  all the tangente
     901                for (i=0;i<nbe;i++) {
     902                        R2 AB = edges[i].v[1]->r -edges[i].v[0]->r;       
     903                        Real8 lAB = Norme2(AB); // length of current edge AB
     904                        Real8 ltg2[2];
     905                        ltg2[0]=0;ltg2[1]=0;
     906                        for (jj=0;jj<2;jj++) {
     907                                R2 tg =  edges[i].tg[jj];
     908                                Real8 ltg = Norme2(tg); // length of tg
     909                                if(ltg == 0) {// no tg
     910                                        if( ! edges[i].v[jj]->Corner())   { // not a Corner       
     911                                                tg =  edges[i].v[1-jj]->r
     912                                                  - edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r;
     913                                                ltg =  Norme2(tg);
     914                                                tg =  tg *(lAB/ltg),ltg=lAB;
     915                                                /*
     916                                                        if(edges[i].ref >=4)
     917                                                        cout << " tg " << tg.x << " "<< tg.y  << " " << edges[i].v[1-jj]->r << edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r << " y-y = "
     918                                                        << edges[i].v[1-jj]->r.y -edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r.y <<  endl;
     919                                                        */
     920                                        }
     921
     922                                        //else ;// a Corner with no tangent => nothing to do   
     923                                } // a tg
     924                                else
     925                                 tg = tg *(lAB/ltg),ltg=lAB;
     926                                ltg2[jj] = ltg;
     927                                if ( (tg,AB) < 0)
     928                                 tg = -tg;
     929                                //if(edges[i].ref >=4) cout << " tg = " << tg << endl;
     930                                edges[i].tg[jj] = tg;
     931                        }     // for (jj=0;jj<2;jj++)
     932
     933                        if (ltg2[0]!=0) edges[i].SetTgA();
     934                        if (ltg2[1]!=0) edges[i].SetTgB();
     935                } // for (i=0;i<nbe;i++)
     936
     937                if(verbosity>7)
     938                 for (i=0;i<nbv;i++)
     939                  if (vertices[i].Required())
     940                        cout << "     The  geo  vertices " << i << " is required " << endl;
     941
     942                for (int step=0;step<2;step++)
     943                  {
     944                        for (i=0;i<nbe;i++)
     945                         edges[i].SetUnMark();
     946
     947                        NbOfCurves = 0;
     948                        Int4  nbgem=0;
     949                        for (int level=0;level < 2 && nbgem != nbe;level++)
     950                         for (i=0;i<nbe;i++) {
     951                                 GeometricalEdge & ei = edges[i];   
     952                                 for(jj=0;jj<2;jj++)
     953                                  if (!ei.Mark() && (level || ei[jj].Required())) {
     954                                          // warning ei.Mark() can be change in loop for(jj=0;jj<2;jj++)
     955                                          int k0=jj,k1;
     956                                          GeometricalEdge *e = & ei;
     957                                          GeometricalVertex *a=(*e)(k0); // begin
     958                                          if(curves) {
     959                                                  curves[NbOfCurves].be=e;
     960                                                  curves[NbOfCurves].kb=k0;
     961                                          }
     962                                          int nee=0;
     963                                          for(;;) {
     964                                                  nee++;
     965                                                  k1 = 1-k0; // next vertex of the edge
     966                                                  e->SetMark();
     967                                                  nbgem++;
     968                                                  e->CurveNumber=NbOfCurves;
     969                                                  if(curves) {
     970                                                          curves[NbOfCurves].ee=e;
     971                                                          curves[NbOfCurves].ke=k1;
     972                                                  }
     973
     974                                                  GeometricalVertex *b=(*e)(k1);
     975                                                  if (a == b ||  b->Required() ) break;
     976                                                  k0 = e->SensAdj[k1];//  vertex in next edge
     977                                                  e = e->Adj[k1]; // next edge
     978
     979                                          }// for(;;)
     980                                          if(verbosity>10 && curves==0) cout << NbOfCurves <<" curve :  nb edges=  "<< nee<<  endl;
     981                                          NbOfCurves++;
     982                                          if(level) {
     983                                                  if(verbosity>4)
     984                                                        cout << "    Warning: Curve "<< NbOfCurves << " without required vertex "
     985                                                          << "so the vertex " << Number(a) << " become required " <<endl;
     986                                                  a->SetRequired();
     987                                          }
     988
     989                                  }}
     990                                 assert(nbgem && nbe);
     991
     992                                 if(step==0) {
     993                                         curves = new Curve[NbOfCurves];
     994                                 }
     995                  }
     996                for(int i=0;i<NbOfCurves ;i++)
     997                  {
     998                        GeometricalEdge * be=curves[i].be, *eqbe=be->link;
     999                        //GeometricalEdge * ee=curves[i].ee, *eqee=be->link;
     1000                        curves[i].master=true;
     1001                        if(be->Equi() || be->ReverseEqui() )
     1002                          {
     1003                                assert(eqbe);
     1004                                int nc = eqbe->CurveNumber;
     1005                                assert(i!=nc);
     1006                                curves[i].next=curves[nc].next;
     1007                                curves[i].master=false;
     1008                                curves[nc].next=curves+i;
     1009                                if(be->ReverseEqui())
     1010                                 curves[i].Reverse();           
     1011                          }
     1012                  }
     1013
     1014                if(verbosity>3)
     1015                 cout << "    End ReadGeometry: Number of curves in geometry is " << NbOfCurves <<endl;
     1016                if(verbosity>4)
     1017                 for(int i=0;i<NbOfCurves ;i++)
     1018                        {
     1019                         cout << " Curve " << i << " begin e=" << Number(curves[i].be) << " k=" << curves[i].kb
     1020                                << "  end e= " << Number(curves[i].ee) << " k=" << curves[i].ke << endl;
     1021                        }
     1022                delete []ev;
     1023                delete []hv;
     1024                delete []eangle;
     1025
    7201026        }
    721       }
    722      
    723       //  verification of cracked edge
    724       int err =0;
    725       for (i=0;i<nbe;i++)
    726         if (edges[i].Cracked() )
    727           {
    728             //    verification of crack
    729             GeometricalEdge & e1=edges[i];
    730             GeometricalEdge & e2=*e1.link;
    731             cerr << i << " " << e1[0].The() << " " << e2[0].The() << " " <<  e1[1].The() << " " << e2[1].The() << endl;
    732             if ( e1[0].The() == e2[0].The() && e1[1].The() == e2[1].The() )
    733               {
    734               }
    735             else
    736               if ( e1[0].The() == e2[1].The() && e1[1].The() == e2[0].The() )
    737                 {
    738                 }
    739               else
    740                 {
    741                   err++;
    742                   cerr << " Cracked edges with no same vertex " << &e1-edges << " " << &e2 -edges << endl;
    743                 }
    744           }
    745         else
    746           {
    747             //  if (!edges[i][0].IsThe()) err++;
    748             // if (!edges[i][1].IsThe()) err++;
    749           }
    750       if (err)
    751         {
    752           cerr << " Some vertex was not distint and not on cracked edge " << err<< endl;
    753           MeshError(222);
    754         }
    755     }
    756     if(verbosity>7)
    757       for (i=0;i<nbv;i++)
    758         if (vertices[i].Required())
    759           cout << "     The geo vertices  " << i << " is required" << endl;
    760 
    761     for (i=0;i<nbv;i++)
    762       hv[i]=-1;// empty list
    763 
    764     for (i=0;i<nbe;i++)
    765       {
    766         R2 v10  =  edges[i].v[1]->r -  edges[i].v[0]->r;
    767         Real8 lv10 = Norme2(v10);
    768         if(lv10 == 0) {
    769           cerr << "The length  of " <<i<< "th Egde is 0 " << endl ;
    770           MeshError(1);}
    771         eangle[i] = atan2(v10.y,v10.x)  ; // angle in [ -Pi,Pi ]
    772         if(verbosity>9)
    773           cout << "     angle edge " << i <<" " << eangle[i]*180/Pi<< v10<<endl;
    774         for (jj=0;jj<2;jj++)
    775           { // generation of list
    776             Int4 v =  Number(edges[i].v[jj]);
    777             ev[k] = hv[v];
    778             hv[v] = k++;
    779           }
    780       }
    781     // bulle sort on the angle of edge 
    782     for (i=0;i<nbv;i++) {
    783       int exch = 1,ord =0;     
    784       while (exch) {
    785         exch = 0;
    786         Int4  *p =  hv + i, *po = p;
    787         Int4 n = *p;
    788         register float angleold = -1000 ; // angle = - infini
    789         ord = 0;
    790         while (n >=0)
    791         {
    792           ord++;
    793           register Int4 i1= n /2;
    794           register Int4  j1 = n % 2;
    795           register Int4 *pn = ev + n;
    796           float angle = j1 ? OppositeAngle(eangle[i1]):  eangle[i1];
    797           n = *pn;
    798           if (angleold > angle) // exch to have : po -> pn -> p
    799             exch=1,*pn = *po,*po=*p,*p=n,po = pn;
    800           else //  to have : po -> p -> pn
    801             angleold =  angle, po = p,p  = pn;
    802         }
    803       } // end while (exch)
    804      
    805       if (ord >= 1 )
    806         { /*
    807           Int4 n = hv[i];
    808           while ( n >=0)
    809             { Int4 i1 = n/2,j1 = n%2;
    810             //float a = 180*(j1 ? OppositeAngle(eangle[i1]): eangle[i1])/Pi;
    811             n = ev[n];
    812             }
    813           */
    814         }
    815       if(ord == 2) { // angulare test to find a corner
    816         Int4 n1 = hv[i];
    817         Int4 n2 = ev[n1];
    818         Int4 i1 = n1 /2, i2 = n2/2; // edge number
    819         Int4  j1 = n1 %2, j2 = n2%2; // vertex in the edge
    820         float angle1= j1 ? OppositeAngle(eangle[i1]) :  eangle[i1];
    821         float angle2= !j2 ? OppositeAngle(eangle[i2]) :  eangle[i2];
    822         float da12 = Abs(angle2-angle1);
    823         if(verbosity>9)
    824           cout <<"     check angle " << i << " " << i1 << " " << i2  << " " << 180*(da12)/Pi
    825                << " " << 180*MaximalAngleOfCorner/Pi << vertices[i] << endl;
    826 
    827         if (( da12 >= MaximalAngleOfCorner )
    828             && (da12 <= 2*Pi -MaximalAngleOfCorner)) {
    829           vertices[i].SetCorner() ;
    830           if(verbosity>7)
    831             cout << "     The vertex " << i << " is a corner (angle) "
    832                  << 180*(da12)/ Pi<< " " << 180*MaximalAngleOfCorner/Pi << endl;}
    833         // if the ref a changing then is     SetRequired();
    834        
    835         if (edges[i1].flag != edges[i2].flag || edges[i1].Required())
    836          {
    837           vertices[i].SetRequired();
    838           if(verbosity>7)
    839             cout << "     The vertex " << i << " is Required the flag change (crack or equi, or require)" << endl;}
    840          
    841         if (edges[i1].ref != edges[i2].ref) {
    842           vertices[i].SetRequired();
    843           if(verbosity>7)
    844             cout << "     The vertex " << i << " is Required ref" << endl;}
    845       } ;
    846      
    847       if(ord != 2) {
    848         vertices[i].SetCorner();
    849         if(verbosity>7)
    850           cout << "     the vertex " << i << " is a corner ordre = " << ord << endl;
    851       }
    852       // close the liste around the vertex
    853       { Int4 no=-1, ne = hv[i];
    854         while ( ne >=0)
    855                  ne = ev[no=ne];       
    856             if(no>=0)
    857               ev[no] = hv[i];
    858           } // now the list around the vertex is circular
    859      
    860     } // end for (i=0;i<nbv;i++)
    861  
    862     k =0;
    863     for (i=0;i<nbe;i++)
    864       for (jj=0;jj<2;jj++){
    865             Int4 n1 = ev[k++];
    866             Int4 i1 = n1/2 ,j1=n1%2;
    867             if( edges[i1].v[j1] != edges[i].v[jj])
    868               { cerr << " Bug Adj edge " << i << " " << jj <<
    869                   " et " << i1 << " " << j1 << " k=" << k;
    870                 cerr << Number(edges[i].v[jj]) <<" <> "
    871                      << Number(edges[i1].v[j1])  <<endl;
    872                 cerr << "edge " << Number(edges[i].v[0]) << " "
    873                      << Number(edges[i].v[1]) << endl;
    874             //    cerr << "in file " <<filename <<endl;
    875                 MeshError(1);
    876               }
    877             edges[i1].Adj[j1] = edges + i;
    878             edges[i1].SensAdj[j1] = jj;
    879             if (verbosity>10)
    880               cout << " edges. Adj " << i1 << " " << j1 << " <--- " << i << " " << jj << endl;
    881       }
    882    
    883     // generation of  all the tangente
    884     for (i=0;i<nbe;i++) {
    885         R2 AB = edges[i].v[1]->r -edges[i].v[0]->r;       
    886         Real8 lAB = Norme2(AB); // length of current edge AB
    887         Real8 ltg2[2];
    888         ltg2[0]=0;ltg2[1]=0;
    889         for (jj=0;jj<2;jj++) {
    890              R2 tg =  edges[i].tg[jj];
    891              Real8 ltg = Norme2(tg); // length of tg
    892              if(ltg == 0) {// no tg
    893                if( ! edges[i].v[jj]->Corner())   { // not a Corner       
    894                  tg =  edges[i].v[1-jj]->r
    895                    - edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r;
    896                  ltg =  Norme2(tg);
    897                  tg =  tg *(lAB/ltg),ltg=lAB;
    898                 /*
    899                    if(edges[i].ref >=4)
    900                    cout << " tg " << tg.x << " "<< tg.y  << " " << edges[i].v[1-jj]->r << edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r << " y-y = "
    901                      << edges[i].v[1-jj]->r.y -edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r.y <<  endl;
    902                 */
    903                }
    904                
    905                //else ;// a Corner with no tangent => nothing to do   
    906              } // a tg
    907              else
    908                tg = tg *(lAB/ltg),ltg=lAB;
    909              ltg2[jj] = ltg;
    910              if ( (tg,AB) < 0)
    911                tg = -tg;
    912             //if(edges[i].ref >=4) cout << " tg = " << tg << endl;
    913              edges[i].tg[jj] = tg;
    914         }     // for (jj=0;jj<2;jj++)
    915        
    916         if (ltg2[0]!=0) edges[i].SetTgA();
    917         if (ltg2[1]!=0) edges[i].SetTgB();
    918     } // for (i=0;i<nbe;i++)
    919 
    920     if(verbosity>7)
    921       for (i=0;i<nbv;i++)
    922         if (vertices[i].Required())
    923           cout << "     The  geo  vertices " << i << " is required " << endl;
    924  
    925    for (int step=0;step<2;step++)
    926    {
    927     for (i=0;i<nbe;i++)
    928       edges[i].SetUnMark();
    929    
    930     NbOfCurves = 0;
    931     Int4  nbgem=0;
    932     for (int level=0;level < 2 && nbgem != nbe;level++)
    933       for (i=0;i<nbe;i++) {
    934         GeometricalEdge & ei = edges[i];   
    935         for(jj=0;jj<2;jj++)
    936           if (!ei.Mark() && (level || ei[jj].Required())) {
    937             // warning ei.Mark() can be change in loop for(jj=0;jj<2;jj++)
    938             int k0=jj,k1;
    939             GeometricalEdge *e = & ei;
    940             GeometricalVertex *a=(*e)(k0); // begin
    941             if(curves) {
    942               curves[NbOfCurves].be=e;
    943               curves[NbOfCurves].kb=k0;
    944             }
    945             int nee=0;
    946             for(;;) {
    947                 nee++;
    948               k1 = 1-k0; // next vertex of the edge
    949               e->SetMark();
    950               nbgem++;
    951               e->CurveNumber=NbOfCurves;
    952               if(curves) {
    953               curves[NbOfCurves].ee=e;
    954               curves[NbOfCurves].ke=k1;
    955               }
    956              
    957               GeometricalVertex *b=(*e)(k1);
    958               if (a == b ||  b->Required() ) break;
    959               k0 = e->SensAdj[k1];//  vertex in next edge
    960               e = e->Adj[k1]; // next edge
    961              
    962             }// for(;;)
    963               if(verbosity>10 && curves==0) cout << NbOfCurves <<" curve :  nb edges=  "<< nee<<  endl;
    964             NbOfCurves++;
    965             if(level) {
    966               if(verbosity>4)
    967                 cout << "    Warning: Curve "<< NbOfCurves << " without required vertex "
    968                      << "so the vertex " << Number(a) << " become required " <<endl;
    969               a->SetRequired();
    970             }
    971        
    972           }}
    973           assert(nbgem && nbe);
    974 
    975           if(step==0) {
    976             curves = new Curve[NbOfCurves];
    977           }
    978         }
    979     for(int i=0;i<NbOfCurves ;i++)
    980      {
    981        GeometricalEdge * be=curves[i].be, *eqbe=be->link;
    982        //GeometricalEdge * ee=curves[i].ee, *eqee=be->link;
    983        curves[i].master=true;
    984        if(be->Equi() || be->ReverseEqui() )
    985         {
    986           assert(eqbe);
    987           int nc = eqbe->CurveNumber;
    988           assert(i!=nc);
    989           curves[i].next=curves[nc].next;
    990           curves[i].master=false;
    991           curves[nc].next=curves+i;
    992           if(be->ReverseEqui())
    993            curves[i].Reverse();           
    994         }
    995      }
    996          
    997     if(verbosity>3)
    998       cout << "    End ReadGeometry: Number of curves in geometry is " << NbOfCurves <<endl;
    999     if(verbosity>4)
    1000     for(int i=0;i<NbOfCurves ;i++)
    1001      {
    1002         cout << " Curve " << i << " begin e=" << Number(curves[i].be) << " k=" << curves[i].kb
    1003              << "  end e= " << Number(curves[i].ee) << " k=" << curves[i].ke << endl;
    1004      }
    1005     delete []ev;
    1006     delete []hv;
    1007     delete []eangle;
    1008    
    1009 }
    1010 Geometry::~Geometry()
    1011 {
    1012         long int verbosity=0;
    1013 
    1014   assert(NbRef<=0);
    1015   if(verbosity>9)
    1016     cout << "DELETE      ~Geometry "<< this  << endl;
    1017   if(vertices)  delete [] vertices;vertices=0;
    1018   if(edges)     delete [] edges;edges=0;
    1019  // if(edgescomponante) delete [] edgescomponante; edgescomponante=0;
    1020   if(triangles) delete [] triangles;triangles=0;
    1021   if(quadtree)  delete  quadtree;quadtree=0;
    1022   if(curves)  delete  []curves;curves=0;NbOfCurves=0;
    1023   if(name) delete [] name;name=0;
    1024   if(subdomains) delete [] subdomains;subdomains=0;
    1025 //  if(ordre)     delete [] ordre;
    1026   EmptyGeometry();
    1027 
    1028 }
    1029 
     1027        /*}}}1*/
    10301028
    10311029}

  • issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp

    r2838 r2841  
    1414namespace bamg {
    1515
    16         /*Constructor{{{1*/
     16        /*Constructor*/
     17        /*FUNCTION MatVVP2x2::MatVVP2x2(const MetricAnIso M){{{1*/
    1718        MatVVP2x2::MatVVP2x2(const MetricAnIso M){
    1819                double a11=M.a11,a21=M.a21,a22=M.a22;
    … …  
    4243        }
    4344        /*}}}1*/
     45
    4446}

  • issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp

    r2838 r2841  
    5252        /*}}}1*/
    5353
    54         /*Others*/
     54        /*Methods*/
     55        /*FUNCTION MetricAnIso::IntersectWith{{{1*/
     56        int MetricAnIso::IntersectWith(const MetricAnIso M2) {
     57                int r=0;
     58                MetricAnIso & M1 = *this;
     59                D2xD2 M;
     60                double l1,l2;
     61
     62                ReductionSimultanee(*this,M2,l1,l2,M);
     63                R2 v1(M.x.x,M.y.x);
     64                R2 v2(M.x.y,M.y.y);
     65                double l11=M1(v1,v1);
     66                double l12=M1(v2,v2);
     67                double l21=M2(v1,v1);
     68                double l22=M2(v2,v2);
     69                if ( l11 < l21 )  r=1,l11=l21;
     70                if ( l12 < l22 )  r=1,l12=l22;
     71                if (r) { // change
     72                        D2xD2 M_1(M.inv());
     73                        D2xD2 D(l11,0,0,l12);
     74                        D2xD2 Mi(M_1.t()*D*M_1);
     75                        a11=Mi.x.x;
     76                        a21=0.5*(Mi.x.y+Mi.y.x);
     77                        a22=Mi.y.y; }
     78                        return r;
     79        }
     80        /*}}}1*/
     81
     82        /*Intermediary*/
    5583        /*FUNCTION ReductionSimultanee{{{1*/
    5684        void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
    … …  
    118146        }
    119147        /*}}}1*/
    120         /*FUNCTION MetricAnIso::IntersectWith{{{1*/
    121         int MetricAnIso::IntersectWith(const MetricAnIso M2) {
    122                 int r=0;
    123                 MetricAnIso & M1 = *this;
    124                 D2xD2 M;
    125                 double l1,l2;
    126 
    127                 ReductionSimultanee(*this,M2,l1,l2,M);
    128                 R2 v1(M.x.x,M.y.x);
    129                 R2 v2(M.x.y,M.y.y);
    130                 double l11=M1(v1,v1);
    131                 double l12=M1(v2,v2);
    132                 double l21=M2(v1,v1);
    133                 double l22=M2(v2,v2);
    134                 if ( l11 < l21 )  r=1,l11=l21;
    135                 if ( l12 < l22 )  r=1,l12=l22;
    136                 if (r) { // change
    137                         D2xD2 M_1(M.inv());
    138                         D2xD2 D(l11,0,0,l12);
    139                         D2xD2 Mi(M_1.t()*D*M_1);
    140                         a11=Mi.x.x;
    141                         a21=0.5*(Mi.x.y+Mi.y.x);
    142                         a22=Mi.y.y; }
    143                         return r;
    144         }
    145         /*}}}1*/
    146148        /*FUNCTION LengthInterpole{{{1*/
    147149        Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB) {

  • issm/trunk/src/c/Bamgx/objects/Triangle.cpp

    r2820 r2841  
    1414namespace bamg {
    1515
    16         /*Others*/
     16        /*Methods*/
    1717        /*FUNCTION Triangle::FindBoundaryEdge{{{1*/
    1818        TriangleAdjacent Triangle::FindBoundaryEdge(int i) const

  • issm/trunk/src/c/Bamgx/objects/Triangles.cpp

    r2820 r2841  
    720720        /*}}}1*/
    721721
    722         /*Others*/
     722        /*Methods*/
    723723        /*FUNCTION Triangles::ConsGeometry{{{1*/
    724724        void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges) // construct a geometry if no geo
    … …  
    57625762        /*}}}1*/
    57635763
    5764 } // end of namespace bamg
     5764}

  • issm/trunk/src/c/Makefile.am

    r2840 r2841  
    329329                                        ./Bamgx/objects/Triangle.cpp    \
    330330                                        ./Bamgx/objects/Geometry.cpp    \
     331                                        ./Bamgx/objects/QuadTree.cpp \
    331332                                        ./Bamgx/Metric.h \
    332                                         ./Bamgx/QuadTree.cpp \
    333333                                        ./Bamgx/QuadTree.h \
    334334                                        ./Bamgx/R2.h \
    … …  
    669669                                        ./Bamgx/objects/Triangle.cpp    \
    670670                                        ./Bamgx/objects/Geometry.cpp    \
     671                                        ./Bamgx/objects/QuadTree.cpp \
    671672                                        ./Bamgx/Metric.h \
    672                                         ./Bamgx/QuadTree.cpp \
    673673                                        ./Bamgx/QuadTree.h \
    674674                                        ./Bamgx/R2.h \
Note: See TracChangeset for help on using the changeset viewer.