Changeset 5120


Ignore:
Timestamp:
08/10/10 13:39:21 (15 years ago)
Author:
Mathieu Morlighem
Message:

improved Bamg objects

Location:
issm/trunk/src/c
Files:
21 edited
2 moved

Legend:

Unmodified
Added
Removed

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

    r5114 r5120  
    6969                                        ./objects/Bamg/Mesh.cpp\
    7070                                        ./objects/Bamg/Mesh.h\
    71                                         ./objects/Bamg/MeshVertex.cpp\
    72                                         ./objects/Bamg/MeshVertex.h\
     71                                        ./objects/Bamg/BamgVertex.cpp\
     72                                        ./objects/Bamg/BamgVertex.h\
    7373                                        ./objects/Bamg/VertexOnEdge.h\
    7474                                        ./objects/Bamg/VertexOnEdge.cpp\
    … …  
    605605                                        ./objects/Bamg/Mesh.h\
    606606                                        ./objects/Bamg/Mesh.cpp\
    607                                         ./objects/Bamg/MeshVertex.cpp\
    608                                         ./objects/Bamg/MeshVertex.h\
     607                                        ./objects/Bamg/BamgVertex.cpp\
     608                                        ./objects/Bamg/BamgVertex.h\
    609609                                        ./objects/Bamg/VertexOnEdge.h\
    610610                                        ./objects/Bamg/VertexOnEdge.cpp\

  • issm/trunk/src/c/objects/Bamg/BamgVertex.cpp

    r5100 r5120  
    99
    1010        /*Methods*/
    11         /*FUNCTION MeshVertex::Smoothing{{{1*/
    12         double  MeshVertex::Smoothing(Mesh &Th,const Mesh &BTh,Triangle* &tstart ,double omega){
     11        /*FUNCTION BamgVertex::Smoothing{{{1*/
     12        double  BamgVertex::Smoothing(Mesh &Th,const Mesh &BTh,Triangle* &tstart ,double omega){
    1313                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Smoothing)*/
    1414
    15                 register MeshVertex* s=this;
    16                 MeshVertex &vP = *s,vPsave=vP;
     15                register BamgVertex* s=this;
     16                BamgVertex &vP = *s,vPsave=vP;
    1717
    1818                register Triangle* tbegin= t , *tria = t , *ttc;
    … …  
    2424
    2525                        if (!tria->Hidden(j)){
    26                                 MeshVertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]];
     26                                BamgVertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]];
    2727
    2828                                R2 Q = vQ,QP(P-Q);
    … …  
    7070                vP.i = Th.toI2(PNew);
    7171
    72                 MeshVertex vPnew = vP;
     72                BamgVertex vPnew = vP;
    7373
    7474                int ok=1;
    … …  
    8282                                tria->det =  bamg::det( (*tria)[0],(*tria)[1]  ,(*tria)[2]);
    8383                                if (loop) {
    84                                         MeshVertex *v0,*v1,*v2,*v3;
     84                                        BamgVertex *v0,*v1,*v2,*v3;
    8585                                        if (tria->det<0) ok =1;                       
    8686                                        else if (tria->Quadrangle(v0,v1,v2,v3))
    … …  
    113113        }
    114114        /*}}}1*/
    115         /*FUNCTION MeshVertex::MetricFromHessian{{{1*/
    116         void MeshVertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,double err,BamgOpts* bamgopts){
     115        /*FUNCTION BamgVertex::MetricFromHessian{{{1*/
     116        void BamgVertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,double err,BamgOpts* bamgopts){
    117117                /*Compute Metric from Hessian*/
    118118
    … …  
    192192        }
    193193        /*}}}1*/
    194         /*FUNCTION MeshVertex::Echo {{{1*/
    195 
    196         void MeshVertex::Echo(void){
     194        /*FUNCTION BamgVertex::Echo {{{1*/
     195
     196        void BamgVertex::Echo(void){
    197197
    198198                printf("Vertex:\n");
    … …  
    205205        }
    206206        /*}}}*/
    207         /*FUNCTION MeshVertex::Optim {{{1*/
    208         long MeshVertex::Optim(int i,int koption){
     207        /*FUNCTION BamgVertex::Optim {{{1*/
     208        long BamgVertex::Optim(int i,int koption){
    209209                long ret=0;
    210210                if ( t && (vint >= 0 ) && (vint <3) ){
    … …  
    221221        /*Intermediary*/
    222222        /*FUNCTION QuadQuality{{{1*/
    223         double QuadQuality(const MeshVertex & a,const MeshVertex &b,const MeshVertex &c,const MeshVertex &d) {
     223        double QuadQuality(const BamgVertex & a,const BamgVertex &b,const BamgVertex &c,const BamgVertex &d) {
    224224                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/QuadQuality)*/
    225225

  • issm/trunk/src/c/objects/Bamg/BamgVertex.h

    r5100 r5120  
    1515        class VertexOnEdge;
    1616
    17         class MeshVertex {
     17        class BamgVertex {
    1818
    1919                public:
    20                         I2 i;  // integer coordinates
    21                         R2 r;  // real coordinates
    22                         Metric m;
    23                         long ReferenceNumber;
     20
     21                        I2        i;                 // integer coordinates
     22                        R2        r;                 // real coordinates
     23                        Metric    m;
     24                        long      ReferenceNumber;
    2425                        Direction DirOfSearch;
    25                         short vint;  // the vertex number in triangle; varies between 0 and 2 in t
     26                        short     vint;              // the vertex number in triangle; varies between 0 and 2 in t
     27
    2628                        union {
    27                                 Triangle* t;   // one triangle which is containing the vertex
    28                                 long      color;
    29                                 MeshVertex*   to;  // used in geometry MeshVertex to know the Mesh MeshVertex associated
    30                                 VertexOnGeom* onGeometry;        // if vint == 8; // set with Mesh::SetVertexFieldOn()
    31                                 MeshVertex*       onBackgroundVertex;// if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh()
    32                                 VertexOnEdge* onBackgroundEdge;  // if vint == 32 on Background edge
     29                                Triangle     *t;                    // one triangle which is containing the vertex
     30                                long          color;
     31                                BamgVertex   *to;                   // used in geometry BamgVertex to know the Mesh BamgVertex associated
     32                                VertexOnGeom *onGeometry;           // if vint == 8; // set with Mesh::SetVertexFieldOn()
     33                                BamgVertex   *onBackgroundVertex;   // if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh()
     34                                VertexOnEdge *onBackgroundEdge;     // if vint == 32 on Background edge
    3335                        };
    3436
    … …  
    4749
    4850                        //inline functions
    49                         inline void Set(const MeshVertex &rec,const Mesh & ,Mesh & ){*this=rec;}
     51                        inline void Set(const BamgVertex &rec,const Mesh & ,Mesh & ){*this=rec;}
    5052        };
    5153
    5254        //Intermediary
    53         double QuadQuality(const MeshVertex &,const MeshVertex &,const MeshVertex &,const MeshVertex &);
     55        double QuadQuality(const BamgVertex &,const BamgVertex &,const BamgVertex &,const BamgVertex &);
    5456}
    5557#endif

  • issm/trunk/src/c/objects/Bamg/Edge.h

    r5095 r5120  
    33
    44#include "./include.h"
    5 #include "./MeshVertex.h"
     5#include "./BamgVertex.h"
    66#include "../../include/include.h"
    77#include "../../shared/Exceptions/exceptions.h"
    … …  
    1616
    1717                public:
    18                         MeshVertex* v[2];
     18                        BamgVertex* v[2];
    1919                        long ref;
    2020                        GeometricalEdge* onGeometry;
    … …  
    2222
    2323                        //Operators
    24                         MeshVertex & operator[](int i){return *v[i];};
    25                         MeshVertex * operator()(int i){return v[i];};
     24                        BamgVertex & operator[](int i){return *v[i];};
     25                        BamgVertex * operator()(int i){return v[i];};
    2626                        R2       operator()(double t) const; // return the point
    27                         const MeshVertex & operator[](int i) const { return *v[i];};
     27                        const BamgVertex & operator[](int i) const { return *v[i];};
    2828
    2929                        //Methods
    30                         void ReNumbering(MeshVertex *vb,MeshVertex *ve, long *renu){
     30                        void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
    3131                                if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
    3232                                if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];

  • issm/trunk/src/c/objects/Bamg/GeometricalVertex.h

    r3913 r5120  
    33
    44#include "./include.h"
    5 #include "MeshVertex.h"
     5#include "BamgVertex.h"
    66
    77namespace bamg {
    … …  
    99        class Geometry;
    1010
    11         class GeometricalVertex : public MeshVertex {
     11        class GeometricalVertex : public BamgVertex {
    1212
    1313                public:

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

    r5095 r5120  
    1313
    1414        /*Constructors/Destructors*/
    15         /*FUNCTION  Geometry::Geometry(const Geometry & Gh){{{1*/
     15        Geometry::Geometry(){
     16                Init();
     17        }
     18        Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
     19                Init();
     20                ReadGeometry(bamggeom,bamgopts);
     21                AfterRead();
     22        }
     23        /*FUNCTION  Geometry::Geometry(const Geometry & Gh) (COPY operator){{{1*/
    1624        Geometry::Geometry(const Geometry & Gh) {
    1725                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)*/
    … …  
    2230                quadtree=0;
    2331                vertices = nbv ? new GeometricalVertex[nbv] : NULL;
    24                 triangles = nbt ? new  Triangle[nbt]:NULL;
    2532                edges = nbe ? new GeometricalEdge[nbe]:NULL;
    2633                curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;
    … …  
    3441                for (i=0;i<NbSubDomains;i++)
    3542                 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
    36                 if (nbt);  {
    37                         ISSMERROR("nbt");
    38                 }
    3943        }
    4044        /*}}}1*/
    … …  
    4246        Geometry::~Geometry() {
    4347                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)*/
    44 
    45                 long int verbose=0;
    46 
    47                 if (NbRef>0){
    48                         ISSMERROR("NbRef>0");
    49                 }
    50                 if(vertices)  delete [] vertices;vertices=0;
    51                 if(edges)     delete [] edges;edges=0;
    52                 // if(edgescomponante) delete [] edgescomponante; edgescomponante=0;
    53                 if(triangles) delete [] triangles;triangles=0;
    54                 if(quadtree)  delete  quadtree;quadtree=0;
    55                 if(curves)  delete  []curves;curves=0;NbOfCurves=0;
     48                ISSMASSERT(NbRef<=0);
     49                if(vertices)  delete [] vertices;   vertices=0;
     50                if(edges)     delete [] edges;      edges=0;
     51                if(triangles) delete [] triangles;  triangles=0;
     52                if(quadtree)  delete  quadtree;     quadtree=0;
     53                if(curves)    delete  []curves;     curves=0;NbOfCurves=0;
    5654                if(subdomains) delete [] subdomains;subdomains=0;
    57                 //  if(ordre)     delete [] ordre;
    58                 EmptyGeometry();
     55                Init();
    5956        }
    6057        /*}}}1*/
    … …  
    6562
    6663                int verbose;
    67                 nbiv=nbv=nbvx=0;
    68                 nbe=nbt=nbtx=0;
     64                nbv=0;
     65                nbe=0;
    6966                NbOfCurves=0;
    7067
    … …  
    7269                int i,j,k,n,i1,i2;
    7370
    74                 //initialize some variables
    75                 int Version=1,dim=2;
    76                 nbv=bamggeom->VerticesSize[0];
    77                 nbe=bamggeom->EdgesSize[0];
    78                 nbvx = nbv;
    79                 nbiv = nbv;
     71                /*initialize some variables*/
    8072                verbose=bamgopts->verbose;
     73                nbv  = bamggeom->VerticesSize[0];
     74                nbe  = bamggeom->EdgesSize[0];
    8175
    8276                //some checks
    83                 if (bamggeom->Vertices==NULL){
    84                         ISSMERROR("the domain provided does not contain any vertex");
    85                 }
    86                 if (bamggeom->Edges==NULL){
    87                         ISSMERROR("the domain provided does not contain any edge");
    88                 }
     77                if (bamggeom->Vertices==NULL) ISSMERROR("the domain provided does not contain any vertex");
     78                if (bamggeom->Edges==NULL) ISSMERROR("the domain provided does not contain any edge");
    8979
    9080                //Vertices
    … …  
    9282                        if(verbose>5) printf("      processing Vertices\n");
    9383                        if (bamggeom->VerticesSize[1]!=3) ISSMERROR("Vertices should have 3 columns");
    94                         vertices = new GeometricalVertex[nbvx];
     84                        vertices = new GeometricalVertex[nbv];
    9585                        for (i=0;i<nbv;i++) {
    9686                                vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
    … …  
    10191                                vertices[i].Set();
    10292                        }
    103                         //find domain extrema for pmin,pmax that will define the square box
    104                         //used for by the quadtree
     93                        /*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/
    10594                        pmin =  vertices[0].r;
    10695                        pmax =  vertices[0].r;
    … …  
    111100                                pmax.y = Max(pmax.y,vertices[i].r.y);
    112101                        }
    113                         R2 DD05 = (pmax-pmin)*0.05;
    114                         pmin -=  DD05;
    115                         pmax +=  DD05;
    116                         //coefIcoor is the coefficient used to have coordinates
    117                         //in ]0 1[^2:  x'=coefIcoor*(x-pmin.x)
    118                         coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
    119                         if(coefIcoor <=0){
    120                                 ISSMERROR("coefIcoor should be positive");
    121                         }
     102                        /*Offset pmin and pmax to avoid round-off errors*/
     103                        R2 offset = (pmax-pmin)*0.05;
     104                        pmin -= offset;
     105                        pmax += offset;
     106                        /*coefIcoor is the coefficient used for integer coordinates:
     107                         *                       (x-pmin.x)
     108                         * Icoor x = (2^30 -1) ------------
     109                         *                          D
     110                         * where D is the longest side of the domain (direction x or y)
     111                         * so that (x-pmin.x)/D is in ]0 1[
     112                         */
     113                        coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
     114                        if(coefIcoor<=0) ISSMERROR("coefIcoor should be positive");
    122115                }
    123116                else{
    124                         ISSMERROR("No MeshVertex provided");
     117                        ISSMERROR("No BamgVertex provided");
    125118                }
    126119
    … …  
    427420                double eps=1e-20;
    428421                QuadTree quadtree; // build quadtree to find duplicates
    429                 MeshVertex* v0=vertices;
     422                BamgVertex* v0=vertices;
    430423                GeometricalVertex* v0g=(GeometricalVertex*) (void*)v0;   
    431424
    … …  
    446439
    447440                        //find nearest vertex already present in the quadtree (NULL if empty)
    448                         MeshVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
     441                        BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
    449442
    450443                        //if there is a vertex found that is to close to vertices[i] -> error
    … …  
    454447                                j=vg-v0g;
    455448                                //check that the clostest vertex is not itself...
    456                                 if ( v !=  &(MeshVertex &) vertices[j]){
    457                                         ISSMERROR(" v !=  &(MeshVertex &) vertices[j]");
     449                                if ( v !=  &(BamgVertex &) vertices[j]){
     450                                        ISSMERROR(" v !=  &(BamgVertex &) vertices[j]");
    458451                                }
    459452                                vertices[i].link = vertices + j;
    … …  
    748741                GeometricalEdge* on =start,* pon=0;
    749742                // walk with the cos on geometry
    750                 int k=0;
     743                int counter=0;
    751744                while(pon != on){ 
    752                         k++;
     745                        counter++;
     746                        ISSMASSERT(counter<100);
    753747                        pon = on;
    754                         if (k>=100){
    755                                 ISSMERROR("k>=100");
    756                         }
    757748                        R2 A= (*on)[0];
    758749                        R2 B= (*on)[1];
    … …  
    775766
    776767                printf("Geometry:\n");
    777                 printf("   nbvx (maximum number of vertices) : %i\n",nbvx);
    778                 printf("   nbtx (maximum number of triangles): %i\n",nbtx);
    779768                printf("   nbv  (number of vertices) : %i\n",nbv);
    780                 printf("   nbt  (number of triangles): %i\n",nbt);
    781769                printf("   nbe  (number of edges)    : %i\n",nbe);
    782                 printf("   nbv  (number of initial vertices) : %i\n",nbiv);
    783770                printf("   NbSubDomains: %i\n",NbSubDomains);
    784771                printf("   NbOfCurves: %i\n",NbOfCurves);
    … …  
    797784        }
    798785        /*}}}*/
    799         /*FUNCTION  Geometry::EmptyGeometry(){{{1*/
    800         void Geometry::EmptyGeometry() {
     786        /*FUNCTION  Geometry::Init{{{1*/
     787        void Geometry::Init(void){
    801788                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/
    802789
    803790                NbRef=0;
     791                nbv=0;
     792                nbe=0;
    804793                quadtree=0;
    805794                curves=0;
    806                 // edgescomponante=0;
    807795                triangles=0;
    808796                edges=0;
    809797                vertices=0;
    810798                NbSubDomains=0;
    811                 //  nbtf=0;
    812                 //  BeginOfCurve=0; 
    813                 nbiv=nbv=nbvx=0;
    814                 nbe=nbt=nbtx=0;
    815799                NbOfCurves=0;
    816                 //  BeginOfCurve=0;
    817800                subdomains=0;
    818                 MaxCornerAngle = 10*Pi/180;
     801                MaxCornerAngle = 10*Pi/180; //default is 10 degres
    819802        }
    820803        /*}}}1*/
    821804        /*FUNCTION  Geometry::ProjectOnCurve {{{1*/
    822         GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,MeshVertex &V,VertexOnGeom &GV) const {
     805        GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
    823806                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
    824807                /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/
    … …  
    839822
    840823                //Get the two vertices of the edge
    841                 const MeshVertex &v0=e[0];
    842                 const MeshVertex &v1=e[1];
     824                const BamgVertex &v0=e[0];
     825                const BamgVertex &v1=e[1];
    843826
    844827                //Get position of V0, V1 and vector v0->v1
    … …  
    906889
    907890                if ((*eg0)(sens0)==(GeometricalVertex*)vg0)
    908                  vg0=VertexOnGeom(*(MeshVertex*) vg0,*eg0,sens0); //vg0 = absisce
     891                 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,sens0); //vg0 = absisce
    909892
    910893                if ((*eg1)(sens1)==(GeometricalVertex*)vg1)
    911                  vg1=VertexOnGeom(*(MeshVertex*) vg1,*eg1,sens1);
     894                 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,sens1);
    912895
    913896                double sg;

  • issm/trunk/src/c/objects/Bamg/Geometry.h

    r5091 r5120  
    1919
    2020                public:
    21                         long NbRef;     // counter of ref on the this class if 0 we can delete
    22                         long nbvx,nbtx; // maximum number of vertices
    23                         long nbv,nbt,nbiv,nbe; // number of vertices
    24                         long NbSubDomains; //
    25                         long NbOfCurves;
    26                         GeometricalVertex* vertices;
    27                         Triangle* triangles;
    28                         GeometricalEdge* edges;
    29                         QuadTree* quadtree;
    30                         GeometricalSubDomain* subdomains;
    31                         Curve* curves;
    32                         R2 pmin,pmax; // extrema
    33                         double coefIcoor;  // coef to integer Icoor1;
    34                         double MaxCornerAngle;
     21
     22                        long                  NbRef;                         // counter of ref on the this class if 0 we can delete
     23                        long                  nbv;                           // number of vertices
     24                        long                  nbe;                           // number of edges
     25                        long                  NbSubDomains;
     26                        long                  NbOfCurves;
     27                        GeometricalVertex    *vertices;
     28                        Triangle             *triangles;
     29                        GeometricalEdge      *edges;
     30                        QuadTree             *quadtree;
     31                        GeometricalSubDomain *subdomains;
     32                        Curve                *curves;
     33                        R2                    pmin,pmax;                     // domain extrema coordinates
     34                        double                coefIcoor;                     // coef to integer Icoor1;
     35                        double                MaxCornerAngle;
    3536
    3637                        //Constructor/Destructors
    3738                        ~Geometry();
    38                         Geometry(const Geometry & Gh); //Copy  Operator
    39                         Geometry(int nbg,const Geometry** ag); // intersection operator
     39                        Geometry();
     40                        Geometry(const Geometry & Gh);
     41                        Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
    4042
    4143                        //Operators
    … …  
    4648
    4749                        //Methods
    48                         int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }
    4950                        void Echo();
    5051                        I2 toI2(const R2 & P) const {
    5152                                return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
    52                                                         ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
     53                                                        ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
     54                        }
    5355                        double MinimalHmin() {return 2.0/coefIcoor;}
    5456                        double MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
    5557                        void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
    56                         void EmptyGeometry();
    57                         Geometry() {EmptyGeometry();}// empty Geometry
     58                        void Init(void);
    5859                        void AfterRead();
    59                         Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();}
    6060                        long Number(const GeometricalVertex & t) const  { return &t - vertices;}
    6161                        long Number(const GeometricalVertex * t) const  { return t - vertices;}
    … …  
    6464                        long Number(const Curve * c) const  { return c - curves;}
    6565                        void UnMarkEdges() {for (int i=0;i<nbe;i++) edges[i].SetUnMark();}
    66                         GeometricalEdge *  ProjectOnCurve(const Edge & ,double,MeshVertex &,VertexOnGeom &) const ;
     66                        GeometricalEdge *  ProjectOnCurve(const Edge & ,double,BamgVertex &,VertexOnGeom &) const ;
    6767                        GeometricalEdge *  Containing(const R2 P,  GeometricalEdge * start) const;
    6868                        void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);

  • issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp

    r5095 r5120  
    4444                                double ba,bb;
    4545                                TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
    46                                 MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
     46                                BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
    4747                                NewItem(A,Metric(ba,v0,bb,v1));
    4848                                t=edge;
    … …  
    5050                                if (det(v0.i,v1.i,b)>=0) {
    5151                                        TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
    52                                         MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
     52                                        BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
    5353                                        NewItem(A,Metric(ba,v0,bb,v1));
    5454                                        return;
    … …  
    8080                                        long int verbose=2;
    8181                                        TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
    82                                         MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
     82                                        BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
    8383                                        NewItem(A,Metric(ba,v0,bb,v1));
    8484                                        return;
    … …  
    232232                        lIntTria[Size].x = x;
    233233                        Metric m0,m1,m2;
    234                         register MeshVertex * v;
     234                        register BamgVertex * v;
    235235                        if ((v=(*tt)(0))) m0    = v->m;
    236236                        if ((v=(*tt)(1))) m1    = v->m;
    … …  
    307307        /*}}}1*/
    308308        /*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/
    309         long ListofIntersectionTriangles::NewPoints(MeshVertex* vertices,long &nbv,long nbvx){
     309        long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){
    310310                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
    311311
    … …  
    354354
    355355                        si += sint;
    356                         if ( nbv<nbvx) {
     356                        if ( nbv<maxnbv) {
    357357                                vertices[nbv].r = C;
    358358                                vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);

  • issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h

    r5095 r5120  
    7171                        void  SplitEdge(const Mesh & ,const R2 &,const R2  &,int nbegin=0);
    7272                        double Length();
    73                         long  NewPoints(MeshVertex *,long & nbv,long nbvx);
     73                        long  NewPoints(BamgVertex *,long & nbv,long maxnbv);
    7474                        void  NewSubSeg(GeometricalEdge *e,double s0,double s1){
    7575                                long int verbosity=0;

  • issm/trunk/src/c/objects/Bamg/Mesh.cpp

    r5095 r5120  
    1616
    1717                /*Initialize fields*/
    18                 PreInit(0);
     18                Init(0);
    1919
    2020                /*Read Geometry if provided*/
    … …  
    4545        Mesh::Mesh(double* index,double* x,double* y,int nods,int nels):Gh(*(new Geometry())),BTh(*this){
    4646
    47                 PreInit(0);
     47                Init(0);
    4848                ReadMesh(index,x,y,nods,nels);
    4949                SetIntCoor();
    … …  
    100100                  printf("   number of triangles %i, remove = %i\n",kt,nbInT-kt);
    101101                  printf("   number of New boundary edge %i\n",nbNewBedge);
    102                   long inbvx =k;
    103                   PreInit(inbvx);
     102                  long imaxnbv =k;
     103                  Init(imaxnbv);
    104104                  for (i=0;i<Tho.nbv;i++)
    105105                        if (kk[i]>=0)
    … …  
    110110                                nbv++;
    111111                          }
    112                   if (inbvx != nbv){
    113                           ISSMERROR("inbvx != nbv");
     112                  if (imaxnbv != nbv){
     113                          ISSMERROR("imaxnbv != nbv");
    114114                  }
    115115                  for (i=0;i<Tho.nbt;i++)
    … …  
    154154          }
    155155        /*}}}1*/
    156         /*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long nbvxx) COPY{{{1*/
    157         Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long nbvxx)
     156        /*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in) COPY{{{1*/
     157        Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in)
    158158          : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this)) {
    159159                  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)*/
    160160                  Gh.NbRef++;
    161                   nbvxx = Max(nbvxx,Th.nbv);
     161                  maxnbv_in = Max(maxnbv_in,Th.nbv);
    162162                  long i;
    163163                  // do all the allocation to be sure all the pointer existe
    164164
    165                   PreInit(nbvxx);// to make the allocation
     165                  Init(maxnbv_in);// to make the allocation
    166166                  // copy of triangles
    167167                  nt=Th.nt;
    168168                  nbv = Th.nbv;
    169169                  nbt = Th.nbt;
    170                   nbiv = Th.nbiv;
    171170                  nbe = Th.nbe;
    172171                  NbSubDomains = Th.NbSubDomains;
    … …  
    245244                        else if (BTh.NbRef==0) delete &BTh;
    246245                }
    247                 PreInit(0); // set all to zero
     246                Init(0); // set all to zero
    248247        }
    249248        /*}}}1*/
    … …  
    260259
    261260                nbv=nods;
    262                 nbvx=nbv;
     261                maxnbv=nbv;
    263262                nbt=nels;
    264                 nbiv=nbvx;
    265263
    266264                //Vertices
    267265                if (verbose) printf("Reading vertices (%i)\n",nbv);
    268                 vertices=(MeshVertex*)xmalloc(nbv*sizeof(MeshVertex));
    269                 ordre=(MeshVertex**)xmalloc(nbv*sizeof(MeshVertex*));
     266                vertices=(BamgVertex*)xmalloc(nbv*sizeof(BamgVertex));
     267                ordre=(BamgVertex**)xmalloc(nbv*sizeof(BamgVertex*));
    270268                for (i=0;i<nbv;i++){
    271269                        vertices[i].r.x=x[i];
    … …  
    276274                        vertices[i].color=0;
    277275                }
    278                 nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals
     276                nbtx=2*maxnbv-2; // for filling The Holes and quadrilaterals
    279277
    280278                //Triangles
    … …  
    311309
    312310                nbv=bamgmesh->VerticesSize[0];
    313                 nbvx=nbv;
     311                maxnbv=nbv;
    314312                nbt=bamgmesh->TrianglesSize[0];
    315                 nbiv=nbvx;
    316313
    317314                //Vertices
    … …  
    319316                        if(verbose>5) printf("      processing Vertices\n");
    320317
    321                         vertices=(MeshVertex*)xmalloc(nbv*sizeof(MeshVertex));
    322                         ordre=(MeshVertex**)xmalloc(nbv*sizeof(MeshVertex*));
     318                        vertices=(BamgVertex*)xmalloc(nbv*sizeof(BamgVertex));
     319                        ordre=(BamgVertex**)xmalloc(nbv*sizeof(BamgVertex*));
    323320
    324321                        for (i=0;i<nbv;i++){
    … …  
    330327                                vertices[i].color=0;
    331328                        }
    332                         nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals
     329                        nbtx=2*maxnbv-2; // for filling The Holes and quadrilaterals
    333330                }
    334331                else{
    … …  
    450447                        for (i=0;i<nbe;i++){
    451448                                for (j=0;j<2;j++) {
    452                                         MeshVertex *v=edges[i].v[j];
     449                                        BamgVertex *v=edges[i].v[j];
    453450                                        long i0=v->color,j0;
    454451                                        if(i0==-1){
    … …  
    726723                                Triangle &t =triangles[i];
    727724                                Triangle* ta;
    728                                 MeshVertex *v0,*v1,*v2,*v3;
     725                                BamgVertex *v0,*v1,*v2,*v3;
    729726                                if (reft[i]<0) continue;
    730727                                if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) {
    … …  
    775772                                VertexOnGeom &v=VerticesOnGeomVertex[i];
    776773                                ISSMASSERT(v.OnGeomVertex());
    777                                 bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((MeshVertex*)v)+1; //back to Matlab indexing
     774                                bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((BamgVertex*)v)+1; //back to Matlab indexing
    778775                                bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number((GeometricalVertex*)v)+1; //back to Matlab indexing
    779776                        }
    … …  
    791788                                        ISSMERROR("A vertices supposed to be OnGeometricEdge is actually not");
    792789                                }
    793                                 bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((MeshVertex*)v)+1; //back to Matlab indexing
     790                                bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((BamgVertex*)v)+1; //back to Matlab indexing
    794791                                bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing
    795792                                bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v; //absisce
    … …  
    10201017                        for (int j=0;j<2;j++){
    10211018
    1022                                 MeshVertex V;
     1019                                BamgVertex V;
    10231020                                VertexOnGeom GV;
    10241021                                Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
    … …  
    10661063        /*}}}1*/
    10671064        /*FUNCTION Mesh::AddVertex{{{1*/
    1068         void Mesh::AddVertex( MeshVertex &s,Triangle* t, Icoor2* det3) {
     1065        void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3) {
    10691066                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)*/
    10701067                // -------------------------------------------
    … …  
    10861083                Triangle* tt[3];
    10871084                //three vertices of t
    1088                 MeshVertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2];
     1085                BamgVertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2];
    10891086                //three determinants
    10901087                Icoor2 det3local[3];
    … …  
    14351432                        for (j=0;j<2;j++){
    14361433                                //get current vertex
    1437                                 MeshVertex* v=edges[i].v[j];
     1434                                BamgVertex* v=edges[i].v[j];
    14381435                                //get vertex color (i0)
    14391436                                long i0=v->color;
    … …  
    15751572                for (i=0;i<nbv;i++){
    15761573                        if((j=colorV[i])>=0){
    1577                                 MeshVertex & v = Gh.vertices[j];
     1574                                BamgVertex & v = Gh.vertices[j];
    15781575                                v = vertices[i];
    15791576                                v.color =0;
    … …  
    23262323                delete [] Edgeflags;
    23272324
    2328                 //Reset MeshVertex to On
     2325                //Reset BamgVertex to On
    23292326                SetVertexFieldOn();
    23302327
    … …  
    25572554                                                ISSMERROR("&e==NULL");
    25582555                                        }
    2559                                         MeshVertex * v0 =  e(0),*v1 = e(1);
     2556                                        BamgVertex * v0 =  e(0),*v1 = e(1);
    25602557                                        Triangle *t  = v0->t;
    25612558                                        int sens = Gh.subdomains[i].sens;
    … …  
    26452642
    26462643                        /*Call NearestVertex*/
    2647                         MeshVertex *a = quadtree->NearestVertex(B.x,B.y) ;
     2644                        BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
    26482645
    26492646                        /*Check output (Vertex a)*/
    … …  
    27222719        }
    27232720        /*}}}1*/
    2724         /*FUNCTION Mesh::GeomToTriangles0{{{1*/
    2725         void Mesh::GeomToTriangles0(long inbvx,BamgOpts* bamgopts){
     2721        /*FUNCTION Mesh::Init{{{1*/
     2722        void Mesh::Init(long maxnbv_in) {
     2723                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/
     2724
     2725                long int verbose=0;
     2726
     2727                srand(19999999);
     2728                NbRef=0;
     2729                NbOfTriangleSearchFind =0;
     2730                NbOfSwapTriangle =0;
     2731                nbv=0;
     2732                maxnbv=maxnbv_in;
     2733                nbt=0;
     2734                NbOfQuad = 0;
     2735                nbtx=2*maxnbv_in-2;
     2736                NbSubDomains=0;
     2737                NbVertexOnBThVertex=0;
     2738                NbVertexOnBThEdge=0;
     2739                VertexOnBThVertex=NULL;
     2740                VertexOnBThEdge=NULL;
     2741
     2742                NbCrackedVertices=0;
     2743                CrackedVertices  =NULL; 
     2744                NbCrackedEdges   =0;
     2745                CrackedEdges     =NULL; 
     2746                nbe = 0;
     2747
     2748                if (maxnbv_in) {
     2749                        vertices=new BamgVertex[maxnbv];
     2750                        ISSMASSERT(vertices);
     2751                        ordre=new (BamgVertex* [maxnbv]);
     2752                        ISSMASSERT(ordre);
     2753                        triangles=new Triangle[nbtx];
     2754                        ISSMASSERT(triangles);
     2755                }
     2756                else {
     2757                        vertices=NULL;
     2758                        ordre=NULL;
     2759                        triangles=NULL;
     2760                        nbtx=0;
     2761                }
     2762
     2763                quadtree=NULL;
     2764                edges=NULL;
     2765                VerticesOnGeomVertex=NULL;
     2766                VerticesOnGeomEdge=NULL;
     2767                NbVerticesOnGeomVertex=0;
     2768                NbVerticesOnGeomEdge=0;
     2769                subdomains=NULL;
     2770                NbSubDomains=0;
     2771        }
     2772        /*}}}1*/
     2773        /*FUNCTION Mesh::Insert{{{1*/
     2774        void Mesh::Insert() {
     2775                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/
     2776
     2777                /*Insert points in the existing Geometry*/
     2778
     2779                //Intermediary
     2780                int i;
     2781
     2782                /*Get options*/
     2783                long int verbose=2;
     2784
     2785                //Display info
     2786                if (verbose>2) printf("   Insert initial %i vertices\n",nbv);
     2787
     2788                //Compute integer coordinates and determinants for the existing vertices (from Geometry)
     2789                SetIntCoor();
     2790
     2791                /*Now we want to build a list (ordre) of the vertices in a random
     2792                 * order. To do so, we use the following method:
     2793                 *
     2794                 * From an initial k0 in [0 nbv[ random (vertex number)
     2795                 * the next k (vertex number) is computed using a big
     2796                 * prime number (PN>>nbv) following:
     2797                 *
     2798                 * k_{i+1} = k_i + PN  [nbv]
     2799                 *
     2800                 * let's show that:
     2801                 *
     2802                 *   for all j in [0 nbv[, ∃! i in [0 nbv[ such that k_i=j
     2803                 *
     2804                 * Let's assume that there are 2 distinct j1 and j2 such that
     2805                 * k_j1 = k_j2
     2806                 *
     2807                 * This means that
     2808                 * 
     2809                 *  k0+j1*PN = k0+j2*PN [nbv]
     2810                 *  (j1-j2)*PN =0       [nbv]
     2811                 * since PN is a prime number larger than nbv, and nbv!=1
     2812                 *  j1-j2=0             [nbv]
     2813                 * BUT
     2814                 *  j1 and j2 are in [0 nbv[ which is impossible.
     2815                 *
     2816                 *  We hence have built a random list of nbv elements of
     2817                 *  [0 nbv[ all distincts*/
     2818                for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ;
     2819                const long PrimeNumber= BigPrimeNumber(nbv) ;
     2820                int   k0=rand()%nbv;
     2821                for (int is3=0; is3<nbv; is3++){
     2822                        ordre[is3]= &vertices[k0=(k0+PrimeNumber)%nbv];
     2823                }
     2824
     2825                /*Modify ordre such that the first 3 vertices form a triangle*/
     2826
     2827                //get first vertex i such that [0,1,i] are not aligned
     2828                for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){
     2829                        //if i is higher than nbv, it means that all the determinants are 0,
     2830                        //all vertices are aligned!
     2831                        if  ( ++i >= nbv) {
     2832                                ISSMERROR("all the vertices are aligned");
     2833                        }
     2834                }
     2835                // exchange i et 2 in "ordre" so that
     2836                // the first 3 vertices are not aligned (real triangle)
     2837                Exchange(ordre[2], ordre[i]);
     2838
     2839                /*Take the first edge formed by the first two vertices and build
     2840                 * the initial simple mesh from this edge and 2 boundary triangles*/
     2841
     2842                BamgVertex *  v0=ordre[0], *v1=ordre[1];
     2843
     2844                nbt = 2;
     2845                triangles[0](0) = NULL; //infinite vertex
     2846                triangles[0](1) = v0;
     2847                triangles[0](2) = v1;
     2848                triangles[1](0) = NULL;//infinite vertex
     2849                triangles[1](2) = v0;
     2850                triangles[1](1) = v1;
     2851
     2852                //Build adjacence
     2853                const int e0 = OppositeEdge[0];
     2854                const int e1 = NextEdge[e0];
     2855                const int e2 = PreviousEdge[e0];
     2856                triangles[0].SetAdj2(e0, &triangles[1] ,e0);
     2857                triangles[0].SetAdj2(e1, &triangles[1] ,e2);
     2858                triangles[0].SetAdj2(e2, &triangles[1] ,e1);
     2859
     2860                triangles[0].det = -1;  //boundary triangle: det = -1
     2861                triangles[1].det = -1;  //boundary triangle: det = -1
     2862
     2863                triangles[0].SetTriangleContainingTheVertex();
     2864                triangles[1].SetTriangleContainingTheVertex();
     2865
     2866                triangles[0].link=&triangles[1];
     2867                triangles[1].link=&triangles[0];
     2868
     2869                //build quadtree
     2870                if (!quadtree)  quadtree = new QuadTree(this,0);
     2871                quadtree->Add(*v0);
     2872                quadtree->Add(*v1);
     2873
     2874                /*Now, add the vertices One by One*/
     2875
     2876                long NbSwap=0;
     2877                if (verbose>3) printf("   Begining of insertion process...\n");
     2878
     2879                for (int icount=2; icount<nbv; icount++) {
     2880
     2881                        //Get new vertex
     2882                        BamgVertex *newvertex=ordre[icount];
     2883
     2884                        //Find the triangle in which newvertex is located
     2885                        Icoor2 dete[3];
     2886                        Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates)
     2887
     2888                        //Add newvertex to the quadtree
     2889                        quadtree->Add(*newvertex);
     2890
     2891                        //Add newvertex to the existing mesh
     2892                        AddVertex(*newvertex,tcvi,dete);
     2893
     2894                        //Make the mesh Delaunay around newvertex by swaping the edges
     2895                        NbSwap += newvertex->Optim(1,0);
     2896                }
     2897
     2898                //Display info
     2899                if (verbose>3) {
     2900                        printf("      NbSwap of insertion: %i\n",NbSwap);
     2901                        printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
     2902                }
     2903
     2904#ifdef NBLOOPOPTIM
     2905
     2906                k0 = rand()%nbv ;
     2907                for (int is4=0; is4<nbv; is4++)
     2908                 ordre[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv];
     2909
     2910                for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){
     2911                        long  NbSwap = 0;
     2912                        for (int is1=0; is1<nbv; is1++)
     2913                         NbSwap += ordre[is1]->Optim(0,0);
     2914                        if (verbose>3) {
     2915                                printf("      Optim Loop: %i\n",Nbloop);
     2916                                printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
     2917                        }
     2918                        if(!NbSwap) break;
     2919                }
     2920                ReMakeTriangleContainingTheVertex();
     2921                // because we break the TriangleContainingTheVertex
     2922#endif
     2923        }
     2924        /*}}}1*/
     2925        /*FUNCTION Mesh::InsertNewPoints{{{1*/
     2926        long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {
     2927                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/
     2928
     2929                long int verbose=0;
     2930                double seuil= 1.414/2 ;// for two close point
     2931                long i;
     2932                long NbSwap=0;
     2933                Icoor2 dete[3];
     2934
     2935                //number of new points
     2936                const long nbvnew=nbv-nbvold;
     2937
     2938                //display info if required
     2939                if (verbose>5) printf("      Try to Insert %i new points\n",nbvnew);
     2940
     2941                //return if no new points
     2942                if (!nbvnew) return 0;
     2943
     2944                /*construction of a random order*/
     2945                const long PrimeNumber= BigPrimeNumber(nbv)  ;
     2946                //remainder of the division of rand() by nbvnew
     2947                long k3 = rand()%nbvnew;
     2948                //loop over the new points
     2949                for (int is3=0; is3<nbvnew; is3++){
     2950                        register long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
     2951                        register long i=nbvold+is3;
     2952                        ordre[i]= vertices + j;
     2953                        ordre[i]->ReferenceNumber=i;
     2954                }
     2955
     2956                // for all the new point
     2957                long iv=nbvold;
     2958                for (i=nbvold;i<nbv;i++){
     2959                        BamgVertex &vi=*ordre[i];
     2960                        vi.i=toI2(vi.r);
     2961                        vi.r=toR2(vi.i);
     2962                        double hx,hy;
     2963                        vi.m.Box(hx,hy);
     2964                        Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);
     2965                        if (!quadtree->ToClose(vi,seuil,hi,hj)){
     2966                                // a good new point
     2967                                BamgVertex &vj = vertices[iv];
     2968                                long  j=vj.ReferenceNumber;
     2969                                if (&vj!=ordre[j]){
     2970                                        ISSMERROR("&vj!= ordre[j]");
     2971                                }
     2972                                if(i!=j){
     2973                                        Exchange(vi,vj);
     2974                                        Exchange(ordre[j],ordre[i]);
     2975                                }
     2976                                vj.ReferenceNumber=0;
     2977                                Triangle *tcvj=FindTriangleContaining(vj.i,dete);
     2978                                if (tcvj && !tcvj->link){
     2979                                        tcvj->Echo();
     2980                                        ISSMERROR("problem inserting point in InsertNewPoints (tcvj=%p and tcvj->link=%i)",tcvj,tcvj->link);
     2981                                }
     2982                                quadtree->Add(vj);
     2983                                AddVertex(vj,tcvj,dete);
     2984                                NbSwap += vj.Optim(1);         
     2985                                iv++;
     2986                        }
     2987                }
     2988                if (verbose>3) {
     2989                        printf("         number of new points: %i\n",iv);
     2990                        printf("         number of to close (?) points: %i\n",nbv-iv);
     2991                        printf("         number of swap after: %i\n",NbSwap);
     2992                }
     2993                nbv = iv;
     2994
     2995                for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1); 
     2996                if (verbose>3) printf("   NbSwap=%i\n",NbSwap);
     2997
     2998                NbTSwap +=  NbSwap ;
     2999                return nbv-nbvold;
     3000        }
     3001        /*}}}1*/
     3002        /*FUNCTION Mesh::MakeGeometricalEdgeToEdge{{{1*/
     3003        Edge** Mesh::MakeGeometricalEdgeToEdge() {
     3004                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeometricalEdgeToEdge)*/
     3005
     3006                if (!Gh.nbe){
     3007                        ISSMERROR("!Gh.nbe");
     3008                }
     3009                Edge **e= new (Edge* [Gh.nbe]);
     3010
     3011                long i;
     3012                for ( i=0;i<Gh.nbe ; i++)
     3013                 e[i]=NULL;
     3014                for ( i=0;i<nbe ; i++)
     3015                  {
     3016                        Edge * ei = edges+i;
     3017                        GeometricalEdge *onGeometry = ei->onGeometry;
     3018                        e[Gh.Number(onGeometry)] = ei;   
     3019                  }
     3020                for ( i=0;i<nbe ; i++)
     3021                 for (int ii=0;ii<2;ii++) {
     3022                         Edge * ei = edges+i;
     3023                         GeometricalEdge *onGeometry = ei->onGeometry;
     3024                         int j= ii;
     3025                         while (!(*onGeometry)[j].Required()) {
     3026                                 Adj(onGeometry,j); // next geom edge
     3027                                 j=1-j;
     3028                                 if (e[Gh.Number(onGeometry)])  break; // optimisation
     3029                                 e[Gh.Number(onGeometry)] = ei;
     3030                         }
     3031                 }
     3032
     3033                int kk=0;
     3034                for ( i=0;i<Gh.nbe ; i++){
     3035                        if (!e[i]){
     3036                                kk++;
     3037                                if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i);
     3038                        }
     3039                }
     3040                if(kk) ISSMERROR("See above");
     3041
     3042                return e;
     3043        }
     3044        /*}}}1*/
     3045        /*FUNCTION Mesh::MakeQuadrangles{{{1*/
     3046        void Mesh::MakeQuadrangles(double costheta){
     3047                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadrangles)*/
     3048
     3049                long int verbose=0;
     3050
     3051                if (verbose>2) printf("MakeQuadrangles costheta = %g\n",costheta);
     3052
     3053                if (costheta >1) {
     3054                        if (verbose>5) printf("   do nothing: costheta > 1\n");
     3055                }
     3056
     3057
     3058                        long nbqq = (nbt*3)/2;
     3059                        DoubleAndInt *qq = new DoubleAndInt[nbqq];
     3060
     3061                        long i,ij;
     3062                        int j;
     3063                        long k=0;
     3064                        for (i=0;i<nbt;i++)
     3065                         for (j=0;j<3;j++)
     3066                          if ((qq[k].q=triangles[i].QualityQuad(j))>=costheta)
     3067                                qq[k++].i3j=i*3+j;
     3068                        //  sort  qq
     3069                        HeapSort(qq,k);
     3070
     3071                        long kk=0;
     3072                        for (ij=0;ij<k;ij++) {
     3073                                i=qq[ij].i3j/3;
     3074                                j=(int) (qq[ij].i3j%3);
     3075                                // optisamition no float computation 
     3076                                if (triangles[i].QualityQuad(j,0) >=costheta)
     3077                                 triangles[i].SetHidden(j),kk++;
     3078                          }
     3079                        NbOfQuad = kk;
     3080                        if (verbose>2){
     3081                                printf("   number of quadrilaterals    = %i\n",NbOfQuad);
     3082                                printf("   number of triangles         = %i\n",nbt-NbOutT- NbOfQuad*2);
     3083                                printf("   number of outside triangles = %i\n",NbOutT);
     3084                        }
     3085                        delete [] qq;
     3086        }
     3087        /*}}}1*/
     3088        /*FUNCTION Mesh::MakeQuadTree{{{1*/
     3089        void Mesh::MakeQuadTree() { 
     3090                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadTree)*/
     3091
     3092                long int verbose=0;
     3093                if (  !quadtree )  quadtree = new QuadTree(this);
     3094
     3095        }
     3096        /*}}}1*/
     3097        /*FUNCTION Mesh::MaxSubDivision{{{1*/
     3098        void  Mesh::MaxSubDivision(double maxsubdiv) {
     3099                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/
     3100
     3101                long int verbose=0;
     3102
     3103                const  double maxsubdiv2 = maxsubdiv*maxsubdiv;
     3104                if(verbose>1) printf("   Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv);
     3105                // for all the edges
     3106                // if the len of the edge is to long
     3107                long it,nbchange=0;   
     3108                double lmax=0;
     3109                for (it=0;it<nbt;it++){
     3110                        Triangle &t=triangles[it];
     3111                        for (int j=0;j<3;j++){
     3112                                Triangle &tt = *t.TriangleAdj(j);
     3113                                if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)){
     3114                                        BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
     3115                                        BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
     3116                                        R2 AB= (R2) v1-(R2) v0;
     3117                                        Metric M = v0;
     3118                                        double l = M(AB,AB);
     3119                                        lmax = Max(lmax,l);
     3120                                        if(l> maxsubdiv2){
     3121                                          R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
     3122                                                double lc = M(AC,AC);
     3123                                                D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
     3124                                                D2xD2 Rt1(Rt.inv());
     3125                                                D2xD2 D(maxsubdiv2,0,0,lc);
     3126                                                D2xD2 MM = Rt1*D*Rt1.t();
     3127                                                v0.m =  M = Metric(MM.x.x,MM.y.x,MM.y.y);
     3128                                                nbchange++;
     3129                                        }
     3130                                        M = v1;
     3131                                        l = M(AB,AB);
     3132                                        lmax = Max(lmax,l);
     3133                                        if(l> maxsubdiv2){
     3134                                          R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
     3135                                                double lc = M(AC,AC);
     3136                                                D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
     3137                                                D2xD2 Rt1(Rt.inv());
     3138                                                D2xD2 D(maxsubdiv2,0,0,lc);
     3139                                                D2xD2  MM = Rt1*D*Rt1.t();
     3140                                                v1.m =  M = Metric(MM.x.x,MM.y.x,MM.y.y);
     3141                                                nbchange++;
     3142                                        }
     3143                                }
     3144                        }
     3145                }
     3146                if(verbose>3){
     3147                        printf("      number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5));
     3148                }
     3149        }
     3150        /*}}}1*/
     3151        /*FUNCTION Mesh::MetricAt{{{1*/
     3152        Metric Mesh::MetricAt(const R2 & A) const {
     3153                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/
     3154
     3155                I2 a = toI2(A);
     3156                Icoor2 deta[3];
     3157                Triangle * t =FindTriangleContaining(a,deta);
     3158                if (t->det <0) { // outside
     3159                        double ba,bb;
     3160                        TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ;
     3161                        return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}
     3162                else { // inside
     3163                        double   aa[3];
     3164                        double s = deta[0]+deta[1]+deta[2];
     3165                        aa[0]=deta[0]/s;
     3166                        aa[1]=deta[1]/s;
     3167                        aa[2]=deta[2]/s;
     3168                        return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);
     3169                }
     3170        }
     3171        /*}}}1*/
     3172/*FUNCTION Mesh::NearestVertex{{{1*/
     3173BamgVertex* Mesh::NearestVertex(Icoor1 i,Icoor1 j) {
     3174        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NearestVertex)*/
     3175        return  quadtree->NearestVertex(i,j);
     3176}
     3177/*}}}1*/
     3178        /*FUNCTION Mesh::NewPoints{{{1*/
     3179        void  Mesh::NewPoints(Mesh & Bh,BamgOpts* bamgopts,int KeepVertices){
     3180                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
     3181
     3182                int i,j,k;
     3183                long NbTSwap=0;
     3184                long nbtold=nbt;
     3185                long nbvold=nbv;
     3186                long Headt=0;
     3187                long next_t;
     3188                long* first_np_or_next_t=new long[nbtx];
     3189                Triangle* t=NULL;
     3190
     3191                /*Recover options*/
     3192                int verbose=bamgopts->verbose;
     3193
     3194                /*First, insert old points if requested*/
     3195                if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< maxnbv)){
     3196                        if (verbose>5) printf("         Inserting initial mesh points\n");
     3197                        for (i=0;i<Bh.nbv;i++){
     3198                                BamgVertex &bv=Bh[i];
     3199                                if (!bv.onGeometry){
     3200                                        vertices[nbv].r   = bv.r;
     3201                                        vertices[nbv++].m = bv.m;
     3202                                }
     3203                        }
     3204                        Bh.ReMakeTriangleContainingTheVertex();     
     3205                        InsertNewPoints(nbvold,NbTSwap)   ;           
     3206                } 
     3207                else Bh.ReMakeTriangleContainingTheVertex();     
     3208
     3209                // generation of the list of next Triangle
     3210                for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);
     3211                // the next traingle of i is -first_np_or_next_t[i]
     3212
     3213                // Big loop (most time consuming)
     3214                int iter=0;
     3215                if (verbose>5) printf("         Big loop\n");
     3216                do {
     3217                        /*Update variables*/
     3218                        iter++;
     3219                        nbtold=nbt;
     3220                        nbvold=nbv;
     3221
     3222                        /*We test all triangles*/
     3223                        i=Headt;
     3224                        next_t=-first_np_or_next_t[i];
     3225                        for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){
     3226
     3227                                //check i
     3228                                if (i<0 || i>=nbt ){
     3229                                        ISSMERROR("Index problem in NewPoints (i=%i not in [0 %i])",i,nbt-1);
     3230                                }
     3231                                //change first_np_or_next_t[i]
     3232                                first_np_or_next_t[i] = iter;
     3233
     3234                                //Loop over the edges of t
     3235                                for(j=0;j<3;j++){
     3236                                        TriangleAdjacent tj(t,j);
     3237                                        BamgVertex &vA = *tj.EdgeVertex(0);
     3238                                        BamgVertex &vB = *tj.EdgeVertex(1);
     3239
     3240                                        //if t is a boundary triangle, or tj locked, continue
     3241                                        if (!t->link)     continue;
     3242                                        if (t->det <0)    continue;
     3243                                        if (t->Locked(j)) continue;
     3244
     3245                                        TriangleAdjacent tadjj = t->Adj(j);       
     3246                                        Triangle* ta=tadjj;
     3247
     3248                                        //if the adjacent triangle is a boundary triangle, continur
     3249                                        if (ta->det<0) continue;         
     3250
     3251                                        R2 A=vA;
     3252                                        R2 B=vB;
     3253                                        k=Number(ta);
     3254
     3255                                        //if this edge has already been done, go to next edge of triangle
     3256                                        if(first_np_or_next_t[k]==iter) continue;
     3257
     3258                                        lIntTria.SplitEdge(Bh,A,B);
     3259                                        lIntTria.NewPoints(vertices,nbv,maxnbv);
     3260                                  } // end loop for each edge
     3261                          }// for triangle   
     3262
     3263                        if (!InsertNewPoints(nbvold,NbTSwap)) break;
     3264                        for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
     3265                        Headt = nbt; // empty list
     3266
     3267                        // for all the triangle containing the vertex i
     3268                        for (i=nbvold;i<nbv;i++){
     3269                                BamgVertex*          s  = vertices + i;
     3270                                TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
     3271                                Triangle*        tbegin= (Triangle*) ta;
     3272                                long kt;
     3273                                do {
     3274                                        kt = Number((Triangle*) ta);
     3275                                        if (first_np_or_next_t[kt]>0){
     3276                                                first_np_or_next_t[kt]=-Headt;
     3277                                                Headt=kt;
     3278                                        }
     3279                                        if (ta.EdgeVertex(0)!=s){
     3280                                                ISSMERROR("ta.EdgeVertex(0)!=s");
     3281                                        }
     3282                                        ta = Next(Adj(ta));
     3283                                } while ( (tbegin != (Triangle*) ta));
     3284                        }
     3285
     3286                } while (nbv!=nbvold);
     3287
     3288                delete []  first_np_or_next_t;
     3289
     3290                long NbSwapf =0,NbSwp;
     3291
     3292                NbSwp = NbSwapf;
     3293                for (i=0;i<nbv;i++)
     3294                 NbSwapf += vertices[i].Optim(0);
     3295                NbTSwap +=  NbSwapf ;
     3296        }
     3297        /*}}}1*/
     3298        /*FUNCTION Mesh::ProjectOnCurve{{{1*/
     3299        GeometricalEdge*   Mesh::ProjectOnCurve( Edge & BhAB, BamgVertex &  vA, BamgVertex & vB,
     3300                                double theta,BamgVertex & R,VertexOnEdge &  BR,VertexOnGeom & GR) {
     3301                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/ProjectOnCurve)*/
     3302
     3303                void *pA=0,*pB=0;
     3304                double tA=0,tB=0;
     3305                R2 A=vA,B=vB;
     3306                BamgVertex * pvA=&vA, * pvB=&vB;
     3307                if (vA.vint == IsVertexOnVertex){
     3308                        pA=vA.onBackgroundVertex;
     3309                }
     3310                else if (vA.vint == IsVertexOnEdge){
     3311                        pA=vA.onBackgroundEdge->be;
     3312                        tA=vA.onBackgroundEdge->abcisse;
     3313                }
     3314                else {
     3315                        ISSMERROR("ProjectOnCurve On BamgVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA));
     3316                }
     3317
     3318                if (vB.vint == IsVertexOnVertex){
     3319                        pB=vB.onBackgroundVertex;
     3320                }
     3321                else if(vB.vint == IsVertexOnEdge){
     3322                        pB=vB.onBackgroundEdge->be;
     3323                        tB=vB.onBackgroundEdge->abcisse;
     3324                }
     3325                else {
     3326                        ISSMERROR("ProjectOnCurve On BamgVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB));
     3327                }
     3328                Edge * e = &BhAB;
     3329                if (!pA || !pB || !e){
     3330                        ISSMERROR("!pA || !pB || !e");
     3331                }
     3332                // be carefull the back ground edge e is on same geom edge
     3333                // of the initiale edge def by the 2 vertex A B;
     3334                //check Is a background Mesh;   
     3335                if (e<BTh.edges || e>=BTh.edges+BTh.nbe){
     3336                        ISSMERROR("e<BTh.edges || e>=BTh.edges+BTh.nbe");
     3337                }
     3338                // walk on BTh edge
     3339                //not finish ProjectOnCurve with BackGround Mesh);
     3340                // 1 first find a back ground edge contening the vertex A
     3341                // 2 walk n back gound boundary to find the final vertex B
     3342
     3343                if( vA.vint == IsVertexOnEdge)
     3344                  { // find the start edge
     3345                        e = vA.onBackgroundEdge->be;     
     3346
     3347                  }
     3348                else if (vB.vint == IsVertexOnEdge)
     3349                  {
     3350                        theta = 1-theta;
     3351                        Exchange(tA,tB);
     3352                        Exchange(pA,pB);
     3353                        Exchange(pvA,pvB);
     3354                        Exchange(A,B);
     3355                        e =  vB.onBackgroundEdge->be;
     3356
     3357                  }
     3358                else{ // do the search by walking
     3359                        ISSMERROR("case not supported yet");
     3360                  }
     3361
     3362                // find the direction of walking with sens of edge and pA,PB;
     3363                R2 AB=B-A;
     3364
     3365                double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);
     3366                int kkk=0;
     3367                int sens = (cosE01AB>0) ? 1 : 0;
     3368
     3369                //   double l=0; // length of the edge AB
     3370                double abscisse = -1;
     3371
     3372                for (int cas=0;cas<2;cas++){
     3373                        // 2 times algo:
     3374                        //    1 for computing the length l
     3375                        //    2 for find the vertex
     3376                        int  iii;
     3377                        BamgVertex  *v0=pvA,*v1;
     3378                        Edge *neee,*eee;
     3379                        double lg =0; // length of the curve
     3380                        double te0;
     3381                        // we suppose take the curve's abcisse
     3382                        for ( eee=e,iii=sens,te0=tA;
     3383                                                eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ;
     3384                                                neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) {
     3385
     3386                                kkk=kkk+1;
     3387                                if (kkk>=100){
     3388                                        ISSMERROR("kkk>=100");
     3389                                }
     3390                                if (!eee){
     3391                                        ISSMERROR("!eee");
     3392                                }
     3393                                double lg0 = lg;
     3394                                double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
     3395                                lg += dp;
     3396                                if (cas && abscisse <= lg) { // ok we find the geom edge
     3397                                        double sss  =   (abscisse-lg0)/dp;
     3398                                        double thetab = te0*(1-sss)+ sss*iii;
     3399                                        if (thetab<0 || thetab>1){
     3400                                                ISSMERROR("thetab<0 || thetab>1");
     3401                                        }
     3402                                        BR = VertexOnEdge(&R,eee,thetab);
     3403                                        return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
     3404                                  }
     3405                          }
     3406                        // we find the end
     3407                        if (v1 != pvB){
     3408                                if (( void*) v1 == pB)
     3409                                 tB = iii;
     3410
     3411                                double lg0 = lg;
     3412                                if (!eee){
     3413                                        ISSMERROR("!eee");
     3414                                }
     3415                                v1 = pvB;
     3416                                double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
     3417                                lg += dp;       
     3418                                abscisse = lg*theta;
     3419                                if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end
     3420                                  { // ok we find the geom edge
     3421                                        double sss  =   (abscisse-lg0)/dp;
     3422                                        double thetab = te0*(1-sss)+ sss*tB;
     3423                                        if (thetab<0 || thetab>1){
     3424                                                ISSMERROR("thetab<0 || thetab>1");
     3425                                        }
     3426                                        BR = VertexOnEdge(&R,eee,thetab);
     3427                                        return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
     3428                                  }
     3429                          }
     3430                        abscisse = lg*theta;
     3431
     3432                  }
     3433                ISSMERROR("Big bug...");
     3434                return 0; // just for the compiler
     3435        }                 
     3436        /*}}}1*/
     3437/*FUNCTION Mesh::ReconstructExistingMesh{{{1*/
     3438void Mesh::ReconstructExistingMesh(){
     3439        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/
     3440
     3441        /*This routine reconstruct an existing mesh to make it CONVEX:
     3442         * -all the holes are filled
     3443         * -concave boundaries are filled
     3444         * A convex mesh is required for a lot of operations. This is why every mesh
     3445         * goes through this process.
     3446         * This routine also generates mesh properties such as adjencies,...
     3447         */
     3448
     3449        /*Intermediary*/
     3450        int verbose=0;
     3451
     3452        // generation of the integer coordinate
     3453
     3454        // find extrema coordinates of vertices pmin,pmax
     3455        long i;
     3456        if(verbose>2) printf("      Reconstruct mesh of %i vertices\n",nbv);
     3457
     3458        //initialize ordre
     3459        ISSMASSERT(ordre);
     3460        for (i=0;i<nbv;i++) ordre[i]=0;
     3461
     3462        //Initialize NbSubDomains
     3463        NbSubDomains =0;
     3464
     3465        /* generation of triangles adjacency*/
     3466
     3467        //First add existing edges
     3468        long kk =0;
     3469        SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
     3470        for (i=0;i<nbe;i++){
     3471                kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
     3472        }
     3473        if (kk != nbe){
     3474                ISSMERROR("There are %i double edges in the mesh",kk-nbe);
     3475        }
     3476
     3477        //Add edges of all triangles in existing mesh
     3478        long* st = new long[nbt*3];
     3479        for (i=0;i<nbt*3;i++) st[i]=-1;
     3480        for (i=0;i<nbt;i++){
     3481                for (int j=0;j<3;j++){
     3482
     3483                        //Add current triangle edge to edge4
     3484                        long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
     3485
     3486                        long invisible=triangles[i].Hidden(j);
     3487
     3488                        //If the edge has not been added to st, add it
     3489                        if(st[k]==-1) st[k]=3*i+j;
     3490
     3491                        //If the edge already exists, add adjacency
     3492                        else if(st[k]>=0) {
     3493                                ISSMASSERT(!triangles[i].TriangleAdj(j));
     3494                                ISSMASSERT(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
     3495
     3496                                triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));
     3497                                if (invisible) triangles[i].SetHidden(j);
     3498                                if (k<nbe)     triangles[i].SetLocked(j);
     3499
     3500                                //Make st[k] negative so that it will throw an error message if it is found again
     3501                                st[k]=-2-st[k];
     3502                        }
     3503
     3504                        //An edge belongs to 2 triangles
     3505                        else {
     3506                                ISSMERROR("The edge (%i , %i) belongs to more than 2 triangles",Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
     3507                        }
     3508                }
     3509        }
     3510
     3511        //Display info if required
     3512        if(verbose>5) {
     3513                printf("         info of Mesh:\n");
     3514                printf("            - number of vertices    = %i \n",nbv);
     3515                printf("            - number of triangles   = %i \n",nbt);
     3516                printf("            - number of given edges = %i \n",nbe);
     3517                printf("            - number of all edges   = %i \n"  ,edge4->nb());
     3518                printf("            - Euler number 1 - nb of holes = %i \n"  ,nbt-edge4->nb()+nbv);
     3519        }
     3520
     3521        //check the consistency of edge[].adj and the geometrical required vertex
     3522        long k=0;
     3523        for (i=0;i<edge4->nb();i++){
     3524                if (st[i]>=0){ // edge alone
     3525                        if (i<nbe){
     3526                                long i0=edge4->i(i);
     3527                                ordre[i0] = vertices+i0;
     3528                                long i1=edge4->j(i);
     3529                                ordre[i1] = vertices+i1;
     3530                        }
     3531                        else {
     3532                                k=k+1;
     3533                                if (k<10) {
     3534                                        //print only 10 edges
     3535                                        printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));
     3536                                }
     3537                                else if (k==10){
     3538                                        printf("Other lost boundary edges not shown...\n");
     3539                                }
     3540                        }
     3541                }
     3542        }
     3543        if(k) {
     3544                ISSMERROR("%i boundary edges (from the geometry) are not defined as mesh edges",k);
     3545        }
     3546
     3547        /* mesh generation with boundary points*/
     3548        long nbvb=0;
     3549        for (i=0;i<nbv;i++){
     3550                vertices[i].t=0;
     3551                vertices[i].vint=0;
     3552                if (ordre[i]) ordre[nbvb++]=ordre[i];
     3553        }
     3554
     3555        Triangle* savetriangles=triangles;
     3556        long savenbt=nbt;
     3557        long savenbtx=nbtx;
     3558        SubDomain* savesubdomains=subdomains;
     3559        subdomains=0;
     3560
     3561        long  Nbtriafillhole=2*nbvb;
     3562        Triangle* triafillhole=new Triangle[Nbtriafillhole];
     3563        triangles = triafillhole;
     3564
     3565        nbt=2;
     3566        nbtx= Nbtriafillhole;
     3567
     3568        //Find a vertex that is not aligned with vertices 0 and 1
     3569        for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;)
     3570         if  (++i>=nbvb) {
     3571                 ISSMERROR("ReconstructExistingMesh: All the vertices are aligned");
     3572         }
     3573        //Move this vertex (i) to the 2d position in ordre
     3574        Exchange(ordre[2], ordre[i]);
     3575
     3576        /*Reconstruct mesh beginning with 2 triangles*/
     3577        BamgVertex *  v0=ordre[0], *v1=ordre[1];
     3578
     3579        triangles[0](0) = NULL; // Infinite vertex
     3580        triangles[0](1) = v0;
     3581        triangles[0](2) = v1;
     3582
     3583        triangles[1](0) = NULL;// Infinite vertex
     3584        triangles[1](2) = v0;
     3585        triangles[1](1) = v1;
     3586        const int e0 = OppositeEdge[0];
     3587        const int e1 = NextEdge[e0];
     3588        const int e2 = PreviousEdge[e0];
     3589        triangles[0].SetAdj2(e0, &triangles[1] ,e0);
     3590        triangles[0].SetAdj2(e1, &triangles[1] ,e2);
     3591        triangles[0].SetAdj2(e2, &triangles[1] ,e1);
     3592
     3593        triangles[0].det = -1;  // boundary triangles
     3594        triangles[1].det = -1;  // boundary triangles
     3595
     3596        triangles[0].SetTriangleContainingTheVertex();
     3597        triangles[1].SetTriangleContainingTheVertex();
     3598
     3599        triangles[0].link=&triangles[1];
     3600        triangles[1].link=&triangles[0];
     3601
     3602        if (!quadtree) delete quadtree; //ReInitialise;
     3603        quadtree = new QuadTree(this,0);
     3604        quadtree->Add(*v0);
     3605        quadtree->Add(*v1);
     3606
     3607        // vertices are added one by one
     3608        long NbSwap=0;
     3609        for (int icount=2; icount<nbvb; icount++) {
     3610                BamgVertex *vi  = ordre[icount];
     3611                Icoor2 dete[3];
     3612                Triangle *tcvi = FindTriangleContaining(vi->i,dete);
     3613                quadtree->Add(*vi);
     3614                AddVertex(*vi,tcvi,dete);
     3615                NbSwap += vi->Optim(1,1);
     3616        }
     3617
     3618        //enforce the boundary
     3619        TriangleAdjacent ta(0,0);
     3620        long nbloss = 0,knbe=0;
     3621        for ( i = 0; i < nbe; i++){
     3622                if (st[i] >=0){ //edge alone => on border
     3623                        BamgVertex &a=edges[i][0], &b=edges[i][1];
     3624                        if (a.t && b.t){
     3625                                knbe++;
     3626                                if (ForceEdge(a,b,ta)<0) nbloss++;
     3627                        }
     3628                }
     3629        }
     3630        if(nbloss) {
     3631                ISSMERROR("we lost %i existing edges other %i",nbloss,knbe);
     3632        }
     3633
     3634        FindSubDomain(1);
     3635        // remove all the hole
     3636        // remove all the good sub domain
     3637        long krm =0;
     3638        for (i=0;i<nbt;i++){
     3639                if (triangles[i].link){ // remove triangles
     3640                        krm++;
     3641                        for (int j=0;j<3;j++){
     3642                                TriangleAdjacent ta =  triangles[i].Adj(j);
     3643                                Triangle &tta = *(Triangle*)ta;
     3644                                //if edge between remove and not remove
     3645                                if(! tta.link){
     3646                                        // change the link of ta;
     3647                                        int ja = ta;
     3648                                        BamgVertex *v0= ta.EdgeVertex(0);
     3649                                        BamgVertex *v1= ta.EdgeVertex(1);
     3650                                        long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
     3651
     3652                                        ISSMASSERT(st[k]>=0);
     3653                                        tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
     3654                                        ta.SetLock();
     3655                                        st[k]=-2-st[k];
     3656                                }
     3657                        }
     3658                }
     3659        }
     3660        long NbTfillHoll =0;
     3661        for (i=0;i<nbt;i++){
     3662                if (triangles[i].link) {
     3663                        triangles[i]=Triangle((BamgVertex *) NULL,(BamgVertex *) NULL,(BamgVertex *) NULL);
     3664                        triangles[i].color=-1;
     3665                }
     3666                else{
     3667                        triangles[i].color= savenbt+ NbTfillHoll++;
     3668                }
     3669        }
     3670        ISSMASSERT(savenbt+NbTfillHoll<=savenbtx);
     3671
     3672        // copy of the outside triangles in saveMesh
     3673        for (i=0;i<nbt;i++){
     3674                if(triangles[i].color>=0) {
     3675                        savetriangles[savenbt]=triangles[i];
     3676                        savetriangles[savenbt].link=0;
     3677                        savenbt++;
     3678                }
     3679        }
     3680        // gestion of the adj
     3681        k =0;
     3682        Triangle * tmax = triangles + nbt;
     3683        for (i=0;i<savenbt;i++) {
     3684                Triangle & ti = savetriangles[i];
     3685                for (int j=0;j<3;j++){
     3686                        Triangle * ta = ti.TriangleAdj(j);
     3687                        int aa = ti.NuEdgeTriangleAdj(j);
     3688                        int lck = ti.Locked(j);
     3689                        if (!ta) k++; // bug
     3690                        else if ( ta >= triangles && ta < tmax){
     3691                                ta= savetriangles + ta->color;
     3692                                ti.SetAdj2(j,ta,aa);
     3693                                if(lck) ti.SetLocked(j);
     3694                        }
     3695                }
     3696        }
     3697
     3698        // restore triangles;
     3699        nbt=savenbt;
     3700        nbtx=savenbtx;
     3701        delete [] triangles;
     3702        delete [] subdomains;
     3703        triangles = savetriangles;
     3704        subdomains = savesubdomains;
     3705        if (k) {
     3706                ISSMERROR("number of triangles edges alone = %i",k);
     3707        }
     3708        FindSubDomain();
     3709
     3710        delete edge4;
     3711        delete [] st;
     3712        for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);
     3713
     3714        SetVertexFieldOn();
     3715
     3716        /*Check requirements consistency*/
     3717        for (i=0;i<nbe;i++){
     3718        /*If the current mesh edge is on Geometry*/
     3719                if(edges[i].onGeometry){
     3720                        for(int j=0;j<2;j++){
     3721                                /*Go through the edges adjacent to current edge (if on the same curve)*/
     3722                                if (!edges[i].adj[j]){
     3723                                        /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/
     3724                                        /*Check that the 2 vertices are on geometry AND required*/
     3725                                        if(!edges[i][j].onGeometry->IsRequiredVertex()){
     3726                                                printf("ReconstructExistingMesh error message: problem with the edge number %i: [%i %i]\n",i+1,Number(edges[i][0])+1,Number(edges[i][1])+1);
     3727                                                printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].onGeometry)+1);
     3728                                                if (edges[i][j].onGeometry->OnGeomVertex())
     3729                                                 printf("the vertex number %i of this edge is a geometric BamgVertex number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->gv)+1);
     3730                                                else if (edges[i][j].onGeometry->OnGeomEdge())
     3731                                                 printf("the vertex number %i of this edge is a geometric Edge number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->ge)+1);
     3732                                                else
     3733                                                 printf("Its pointer is %p\n",edges[i][j].onGeometry);
     3734
     3735                                                printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");
     3736                                                ISSMERROR("See above (might be cryptic...)");
     3737                                        }
     3738                                }
     3739                        }
     3740                }
     3741        }
     3742}
     3743/*}}}1*/
     3744        /*FUNCTION Mesh::ReNumberingTheTriangleBySubDomain{{{1*/
     3745        void Mesh::ReNumberingTheTriangleBySubDomain(bool justcompress){
     3746                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
     3747
     3748                long int verbose=0;
     3749                long *renu= new long[nbt];
     3750                register Triangle *t0,*t,*te=triangles+nbt;
     3751                register long k=0,it,i,j;
     3752
     3753                for ( it=0;it<nbt;it++)
     3754                 renu[it]=-1; // outside triangle
     3755                for ( i=0;i<NbSubDomains;i++)
     3756                  {
     3757                        t=t0=subdomains[i].head;
     3758                        if (!t0){ // not empty sub domain
     3759                                ISSMERROR("!t0");
     3760                        }
     3761                        do {
     3762                                long kt = Number(t);
     3763                                if (kt<0 || kt >= nbt ){
     3764                                        ISSMERROR("kt<0 || kt >= nbt");
     3765                                }
     3766                                if (renu[kt]!=-1){
     3767                                        ISSMERROR("renu[kt]!=-1");
     3768                                }
     3769                                renu[kt]=k++;
     3770                        }
     3771                        while (t0 != (t=t->link));
     3772                  }
     3773                // take is same numbering if possible   
     3774                if(justcompress)
     3775                 for ( k=0,it=0;it<nbt;it++)
     3776                  if(renu[it] >=0 )
     3777                        renu[it]=k++;
     3778
     3779                // put the outside triangles at the end
     3780                for ( it=0;it<nbt;it++){
     3781                        if (renu[it]==-1) renu[it]=k++;
     3782                }
     3783                if (k != nbt){
     3784                        ISSMERROR("k != nbt");
     3785                }
     3786                // do the change on all the pointeur
     3787                for ( it=0;it<nbt;it++)
     3788                 triangles[it].ReNumbering(triangles,te,renu);
     3789
     3790                for ( i=0;i<NbSubDomains;i++)
     3791                 subdomains[i].head=triangles+renu[Number(subdomains[i].head)];
     3792
     3793                // move the Triangles  without a copy of the array
     3794                // be carefull not trivial code
     3795                for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu
     3796                 if (renu[it] >= 0) // a new sub cycle
     3797                        {
     3798                         i=it;
     3799                         Triangle ti=triangles[i],tj;
     3800                         while ( (j=renu[i]) >= 0)
     3801                                { // i is old, and j is new
     3802                                 renu[i] = -1; // mark
     3803                                 tj = triangles[j]; // save new
     3804                                 triangles[j]= ti; // new <- old
     3805                                 i=j;     // next
     3806                                 ti = tj;
     3807                                } 
     3808                        }
     3809                delete [] renu;
     3810                nt = nbt - NbOutT;
     3811
     3812        }
     3813        /*}}}1*/
     3814        /*FUNCTION Mesh::ReNumberingVertex{{{1*/
     3815        void Mesh::ReNumberingVertex(long * renu) {
     3816                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
     3817
     3818                // warning be carfull because pointer
     3819                // from on mesh to over mesh
     3820                //  --  so do ReNumbering at the beginning
     3821                BamgVertex * ve = vertices+nbv;
     3822                long it,ie,i;
     3823
     3824                printf("renumbering triangles\n");
     3825                for ( it=0;it<nbt;it++)
     3826                 triangles[it].ReNumbering(vertices,ve,renu);
     3827
     3828                printf("renumbering edges\n");
     3829                for ( ie=0;ie<nbe;ie++)
     3830                 edges[ie].ReNumbering(vertices,ve,renu);
     3831
     3832                printf("renumbering vertices on geom\n");
     3833                for (i=0;i< NbVerticesOnGeomVertex;i++)
     3834                  {
     3835                        BamgVertex *v = VerticesOnGeomVertex[i].mv;
     3836                        if (v>=vertices && v < ve)
     3837                         VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)];
     3838                  }
     3839
     3840                printf("renumbering vertices on edge\n");
     3841                for (i=0;i< NbVerticesOnGeomEdge;i++)
     3842                  {
     3843                        BamgVertex *v =VerticesOnGeomEdge[i].mv;
     3844                        if (v>=vertices && v < ve)
     3845                         VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)];
     3846                  }
     3847
     3848                printf("renumbering vertices on Bth vertex\n");
     3849                for (i=0;i< NbVertexOnBThVertex;i++)
     3850                  {
     3851                        BamgVertex *v=VertexOnBThVertex[i].v;
     3852                        if (v>=vertices && v < ve)
     3853                         VertexOnBThVertex[i].v=vertices+renu[Number(v)];
     3854                  }
     3855
     3856                for (i=0;i< NbVertexOnBThEdge;i++)
     3857                  {
     3858                        BamgVertex *v=VertexOnBThEdge[i].v;
     3859                        if (v>=vertices && v < ve)
     3860                         VertexOnBThEdge[i].v=vertices+renu[Number(v)];
     3861                  }
     3862
     3863                // move the Vertices without a copy of the array
     3864                // be carefull not trivial code
     3865                long j;
     3866                for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
     3867                 if (renu[it] >= 0) // a new sub cycle
     3868                        {
     3869                         i=it;
     3870                         BamgVertex ti=vertices[i],tj;
     3871                         while ( (j=renu[i]) >= 0){
     3872                                 // i is old, and j is new
     3873                                 renu[i] = -1-renu[i]; // mark
     3874                                 tj = vertices[j];     // save new
     3875                                 vertices[j]= ti;      // new <- old
     3876                                 i=j;     // next
     3877                                 ti = tj;
     3878                                } 
     3879                        }
     3880                if (quadtree){
     3881                        delete quadtree;
     3882                        quadtree = new QuadTree(this);
     3883                }
     3884
     3885                for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
     3886        }
     3887        /*}}}1*/
     3888/*FUNCTION Mesh::SetIntCoor{{{1*/
     3889void Mesh::SetIntCoor(const char * strfrom) {
     3890        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)*/
     3891
     3892        /*Set integer coordinate for existing vertices*/
     3893
     3894        //Get extrema coordinates of the existing vertices
     3895        pmin =  vertices[0].r;
     3896        pmax =  vertices[0].r;
     3897        long i;
     3898        for (i=0;i<nbv;i++) {
     3899                pmin.x = Min(pmin.x,vertices[i].r.x);
     3900                pmin.y = Min(pmin.y,vertices[i].r.y);
     3901                pmax.x = Max(pmax.x,vertices[i].r.x);
     3902                pmax.y = Max(pmax.y,vertices[i].r.y);
     3903        }
     3904        R2 DD = (pmax-pmin)*0.05;
     3905        pmin = pmin-DD;
     3906        pmax = pmax+DD;
     3907
     3908        //Compute coefIcoor
     3909        coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
     3910        if (coefIcoor<=0){
     3911                ISSMERROR("coefIcoor should be positive, a problem in the geometry is likely");
     3912        }
     3913
     3914        // generation of integer coord 
     3915        for (i=0;i<nbv;i++) {
     3916                vertices[i].i = toI2(vertices[i].r);   
     3917        }
     3918
     3919        // computation of the det
     3920        int number_of_errors=0;
     3921        for (i=0;i<nbt;i++) {
     3922                BamgVertex & v0 = triangles[i][0];
     3923                BamgVertex & v1 = triangles[i][1];
     3924                BamgVertex & v2 = triangles[i][2];
     3925
     3926                //If this is not a boundary triangle
     3927                if ( &v0 && &v1 &&  &v2 ){
     3928
     3929                        /*Compute determinant*/
     3930                        triangles[i].det= det(v0,v1,v2);
     3931
     3932                        /*Check that determinant is positive*/
     3933                        if (triangles[i].det <=0){
     3934
     3935                                /*increase number_of_errors and print error only for the first 20 triangles*/
     3936                                number_of_errors++;
     3937                                if (number_of_errors<20){
     3938                                        printf("Area of Triangle %i < 0 (det=%i)\n",i+1,triangles[i].det);
     3939                                }
     3940                        }
     3941                }
     3942
     3943                //else, set as -1
     3944                else triangles[i].det=-1;
     3945        }
     3946
     3947        if (number_of_errors) ISSMERROR("Fatal error: some triangles have negative areas, see above");
     3948}
     3949/*}}}1*/
     3950/*FUNCTION Mesh::ShowRegulaty{{{1*/
     3951void  Mesh::ShowRegulaty() const {
     3952        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/
     3953
     3954        const  double  sqrt32=sqrt(3.)*0.5;
     3955        const double  aireKh=sqrt32*0.5;
     3956        D2  Beq(1,0),Heq(0.5,sqrt32);
     3957        D2xD2 Br(D2xD2(Beq,Heq).t());
     3958        D2xD2 B1r(Br.inv());
     3959        double gammamn=1e100,hmin=1e100;
     3960        double gammamx=0,hmax=0;
     3961        double beta=1e100;
     3962        double beta0=0;
     3963        double  alpha2=0;
     3964        double area=0,Marea=0;
     3965        // double cf= double(coefIcoor);
     3966        // double cf2= 6.*cf*cf;
     3967        int nt=0;
     3968        for (int it=0;it<nbt;it++)
     3969         if ( triangles[it].link)
     3970                {
     3971                 nt++;
     3972                 Triangle &K=triangles[it];
     3973                 double  area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;
     3974                 area+= area3;
     3975                 D2xD2 B_Kt(K[0],K[1],K[2]);
     3976                 D2xD2 B_K(B_Kt.t());
     3977                 D2xD2 B1K = Br*B_K.inv();
     3978                 D2xD2 BK =  B_K*B1r;
     3979                 D2xD2 B1B1 = B1K.t()*B1K;
     3980                 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
     3981                 MatVVP2x2 VMK(MK);
     3982                 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
     3983                 double betaK=0;
     3984
     3985                 for (int j=0;j<3;j++)
     3986                        {
     3987                         double he= Norme2(R2(K[j],K[(j+1)%3]));
     3988                         hmin=Min(hmin,he);
     3989                         hmax=Max(hmax,he);
     3990                         BamgVertex & v=K[j];
     3991                         D2xD2 M((Metric)v);
     3992                         betaK += sqrt(M.det());
     3993
     3994                         D2xD2 BMB = BK.t()*M*BK;
     3995                         Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);
     3996                         MatVVP2x2 VM1(M1);
     3997                         gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
     3998                         gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);         
     3999                        }
     4000                 betaK *= area3;//  1/2 (somme sqrt(det))* area(K)
     4001                 Marea+= betaK;
     4002                 beta=min(beta,betaK);
     4003                 beta0=max(beta0,betaK);
     4004                }   
     4005        area*=3;
     4006        gammamn=sqrt(gammamn);
     4007        gammamx=sqrt(gammamx);   
     4008        printf("   Adaptmesh info:\n");
     4009        printf("      number of triangles = %i\n",nt);
     4010        printf("      hmin = %g, hmax=%g\n",hmin,hmax);
     4011        printf("      area = %g, M area = %g, M area/( |Khat| nt) = %g\n",area,Marea, Marea/(aireKh*nt));
     4012        printf("      infinite-regularity(?): min = %g, max = %g\n",gammamn,gammamx);
     4013        printf("      anisomax = %g, beta max = %g, min = %g\n",pow(alpha2,0.5),1./pow(beta/aireKh,0.5), 1./pow(beta0/aireKh,0.5));
     4014}
     4015/*}}}1*/
     4016/*FUNCTION Mesh::ShowHistogram{{{1*/
     4017void  Mesh::ShowHistogram() const {
     4018        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/
     4019
     4020        const long kmax=10;
     4021        const double llmin = 0.5,llmax=2;
     4022        const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);
     4023        long histo[kmax+1];
     4024        long i,it,k, nbedges =0;
     4025        for (i=0;i<=kmax;i++) histo[i]=0;
     4026        for (it=0;it<nbt;it++)
     4027         if ( triangles[it].link)
     4028                {
     4029
     4030                 for (int j=0;j<3;j++)
     4031                        {
     4032                         Triangle *ta = triangles[it].TriangleAdj(j);
     4033                         if ( !ta || !ta->link || Number(ta) >= it)
     4034                                {
     4035                                 BamgVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
     4036                                 BamgVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
     4037                                 if ( !& vP || !&vQ) continue;
     4038                                 R2 PQ = vQ.r - vP.r;
     4039                                 double l = log(LengthInterpole(vP,vQ,PQ));
     4040                                 nbedges++;
     4041                                 k = (int) ((l - lmin)*delta);
     4042                                 k = Min(Max(k,0L),kmax);
     4043                                 histo[k]++;
     4044                                }
     4045                        }
     4046                } 
     4047        printf(" --- Histogram of the unit mesh,  nb of edges = %i\n",nbedges);
     4048        printf("      length of edge in   | %% of edge  | Nb of edges \n");
     4049        printf("      --------------------+-------------+-------------\n");
     4050        for   (i=0;i<=kmax;i++){
     4051                if (i==0) printf("      %10i",0);
     4052                else      printf("      %10g",exp(lmin+i/delta));
     4053                if (i==kmax) printf("          +inf   ");
     4054                else printf("      %10g",exp(lmin+(i+1)/delta));
     4055                printf("|  %10g |\n",((long)  ((10000.0 * histo[i])/ nbedges))/100.0);
     4056                printf("  %i\n",histo[i]);
     4057        }
     4058        printf("      --------------------+-------------+-------------\n");
     4059}
     4060/*}}}1*/
     4061/*FUNCTION Mesh::SmoothingVertex{{{1*/
     4062void Mesh::SmoothingVertex(int nbiter,double omega ) {
     4063        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/
     4064
     4065        long int verbose=0;
     4066        //  if quatree exist remove it end reconstruct
     4067        if (quadtree) delete quadtree;
     4068        quadtree=0;
     4069        ReMakeTriangleContainingTheVertex();
     4070        Triangle vide; // a triangle to mark the boundary vertex
     4071        Triangle   ** tstart= new Triangle* [nbv];
     4072        long i,j,k;
     4073        //   attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide
     4074        if ( this == & BTh)
     4075         for ( i=0;i<nbv;i++)
     4076          tstart[i]=vertices[i].t;     
     4077        else
     4078         for ( i=0;i<nbv;i++)
     4079          tstart[i]=0;
     4080        for ( j=0;j<NbVerticesOnGeomVertex;j++ )
     4081         tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide;
     4082        for ( k=0;k<NbVerticesOnGeomEdge;k++ )
     4083         tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;
     4084        if(verbose>2) printf("   SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega);
     4085        for (k=0;k<nbiter;k++)
     4086          {
     4087                long i,NbSwap =0;
     4088                double delta =0;
     4089                for ( i=0;i<nbv;i++)
     4090                 if (tstart[i] != &vide) // not a boundary vertex
     4091                  delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
     4092                if (!NbOfQuad)
     4093                 for ( i=0;i<nbv;i++)
     4094                  if (tstart[i] != &vide) // not a boundary vertex
     4095                        NbSwap += vertices[i].Optim(1);
     4096                if (verbose>3) printf("      move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap);
     4097          }
     4098
     4099        delete [] tstart;
     4100        if (quadtree) quadtree= new QuadTree(this);
     4101}
     4102/*}}}1*/
     4103/*FUNCTION Mesh::SmoothMetric{{{1*/
     4104void Mesh::SmoothMetric(double raisonmax) {
     4105        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/
     4106
     4107        long int verbose=0;
     4108
     4109        if(raisonmax<1.1) return;
     4110        if(verbose > 1) printf("   Mesh::SmoothMetric raisonmax = %g\n",raisonmax);
     4111        ReMakeTriangleContainingTheVertex();
     4112        long i,j,kch,kk,ip;
     4113        long *first_np_or_next_t0 = new long[nbv];
     4114        long *first_np_or_next_t1 = new long[nbv];
     4115        long Head0 =0,Head1=-1;
     4116        double logseuil= log(raisonmax);
     4117
     4118        for(i=0;i<nbv-1;i++)
     4119         first_np_or_next_t0[i]=i+1;
     4120        first_np_or_next_t0[nbv-1]=-1;// end;
     4121        for(i=0;i<nbv;i++)
     4122         first_np_or_next_t1[i]=-1;
     4123        kk=0;
     4124        while (Head0>=0&& kk++<100) {
     4125                kch=0;
     4126                for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
     4127                        //  pour tous les triangles autour du sommet s
     4128                        register Triangle* t= vertices[i].t;
     4129                        if (!t){
     4130                                ISSMERROR("!t");
     4131                        }
     4132                        BamgVertex & vi = vertices[i];
     4133                        TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);
     4134                        BamgVertex *pvj0 = ta.EdgeVertex(0);
     4135                        while (1) {
     4136                                ta=Previous(Adj(ta));
     4137                                if (vertices+i != ta.EdgeVertex(1)){
     4138                                        ISSMERROR("vertices+i != ta.EdgeVertex(1)");
     4139                                }
     4140                                BamgVertex & vj = *(ta.EdgeVertex(0));
     4141                                if ( &vj ) {
     4142                                        j= &vj-vertices;
     4143                                        if (j<0 || j >= nbv){
     4144                                                ISSMERROR("j<0 || j >= nbv");
     4145                                        }
     4146                                        R2 Aij = (R2) vj - (R2) vi;
     4147                                        double ll =  Norme2(Aij);
     4148                                        if (0) { 
     4149                                                double hi = ll/vi.m(Aij);
     4150                                                double hj = ll/vj.m(Aij);
     4151                                                if(hi < hj)
     4152                                                  {
     4153                                                        double dh=(hj-hi)/ll;
     4154                                                        if (dh>logseuil) {
     4155                                                                vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
     4156                                                                if(first_np_or_next_t1[j]<0)
     4157                                                                 kch++,first_np_or_next_t1[j]=Head1,Head1=j;
     4158                                                        }
     4159                                                  }
     4160                                        }
     4161                                        else
     4162                                          {
     4163                                                double li = vi.m(Aij);
     4164                                                if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
     4165                                                 if(first_np_or_next_t1[j]<0) // if the metrix change
     4166                                                  kch++,first_np_or_next_t1[j]=Head1,Head1=j;
     4167                                          }
     4168                                }
     4169                                if  ( &vj ==  pvj0 ) break;
     4170                        }
     4171                }
     4172                Head0 = Head1;
     4173                Head1 = -1;
     4174                Exchange(first_np_or_next_t0,first_np_or_next_t1);
     4175        }
     4176        if(verbose>2) printf("      number of iterations = %i\n",kch);
     4177        delete [] first_np_or_next_t0;
     4178        delete [] first_np_or_next_t1;
     4179}
     4180/*}}}1*/
     4181        /*FUNCTION Mesh::SplitElement{{{1*/
     4182        int  Mesh::SplitElement(int choice){
     4183                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/
     4184
     4185                long int verbose=0;
     4186
     4187                Direction NoDirOfSearch;
     4188                const  int withBackground = &BTh != this && &BTh;
     4189
     4190                ReNumberingTheTriangleBySubDomain();
     4191                //int nswap =0;
     4192                const long nfortria( choice ? 4 : 6);
     4193                if(withBackground)
     4194                  {
     4195                        BTh.SetVertexFieldOn();
     4196                        SetVertexFieldOnBTh();
     4197                  }
     4198                else
     4199                 BTh.SetVertexFieldOn();
     4200
     4201                long newnbt=0,newnbv=0;
     4202                long * kedge = 0;
     4203                long newNbOfQuad=NbOfQuad;
     4204                long * ksplit= 0, * ksplitarray=0;
     4205                long kkk=0;
     4206                int ret =0;
     4207                if (maxnbv<nbv+nbe) return 1;//   
     4208                // 1) create  the new points by spliting the internal edges
     4209                // set the
     4210                long nbvold = nbv;
     4211                long nbtold = nbt;
     4212                long NbOutTold  = NbOutT;
     4213                long  NbEdgeOnGeom=0;
     4214                long i;
     4215
     4216                nbt = nbt - NbOutT; // remove all the  the ouside triangles
     4217                long nbtsave = nbt;
     4218                Triangle * lastT = triangles + nbt;
     4219                for (i=0;i<nbe;i++)
     4220                 if(edges[i].onGeometry) NbEdgeOnGeom++;
     4221                long newnbe=nbe+nbe;
     4222                //  long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
     4223                long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;
     4224                // long newNbVertexOnBThVertex=NbVertexOnBThVertex;
     4225                long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;
     4226
     4227                // do allocation for pointeur to the geometry and background
     4228                VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];
     4229                VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ?  new VertexOnEdge[newNbVertexOnBThEdge]:0;
     4230                if (NbVerticesOnGeomEdge)
     4231                 memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);
     4232                if (NbVertexOnBThEdge)
     4233                 memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);
     4234                Edge *newedges = new Edge [newnbe];
     4235                //  memcpy(newedges,edges,sizeof(Edge)*nbe);
     4236                SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);
     4237                long k=nbv;
     4238                long kk=0;
     4239                long kvb = NbVertexOnBThEdge;
     4240                long kvg = NbVerticesOnGeomEdge;
     4241                long ie =0;
     4242                Edge ** edgesGtoB=0;
     4243                if (withBackground)
     4244                 edgesGtoB= BTh.MakeGeometricalEdgeToEdge();
     4245                long ferr=0;
     4246                for (i=0;i<nbe;i++)
     4247                 newedges[ie].onGeometry=0;
     4248
     4249                for (i=0;i<nbe;i++)
     4250                  {
     4251                        GeometricalEdge *ong =  edges[i].onGeometry;
     4252
     4253                        newedges[ie]=edges[i];
     4254                        newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;
     4255                        newedges[ie].adj[1]=newedges + ie +1;
     4256                        R2 A = edges[i][0],B = edges[i][1];
     4257
     4258
     4259                        kk += (i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
     4260                        if (ong) // a geometrical edges
     4261                          {
     4262                                if (withBackground){
     4263                                        // walk on back ground mesh
     4264                                        //  newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge);
     4265                                        // a faire -- difficile
     4266                                        // the first PB is to now a background edge between the 2 vertices
     4267                                        if (!edgesGtoB){
     4268                                                ISSMERROR("!edgesGtoB");
     4269                                        }
     4270                                        ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].onGeometry)],
     4271                                                                edges[i][0],edges[i][1],0.5,vertices[k],
     4272                                                                newVertexOnBThEdge[kvb],
     4273                                                                newVerticesOnGeomEdge[kvg++]);
     4274                                        vertices[k].ReferenceNumber= edges[i].ref;
     4275                                        vertices[k].DirOfSearch =   NoDirOfSearch;       
     4276                                        ;
     4277                                        // get the Info on background mesh
     4278                                        double s =        newVertexOnBThEdge[kvb];
     4279                                        BamgVertex &  bv0  = newVertexOnBThEdge[kvb][0];
     4280                                        BamgVertex &  bv1  = newVertexOnBThEdge[kvb][1];
     4281                                        // compute the metrix of the new points
     4282                                        vertices[k].m =  Metric(1-s,bv0,s,bv1);
     4283                                        kvb++;
     4284                                  }
     4285                                else
     4286                                  {
     4287                                        ong=Gh.ProjectOnCurve(edges[i],
     4288                                                                0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);
     4289                                        // vertices[k].i = toI2( vertices[k].r);
     4290                                        vertices[k].ReferenceNumber = edges[i].ref;
     4291                                        vertices[k].DirOfSearch = NoDirOfSearch;
     4292                                        vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);             
     4293                                  } 
     4294                          }
     4295                        else // straigth line edge ---
     4296                          {
     4297                                vertices[k].r = ((R2) edges[i][0] + (R2)  edges[i][1] )*0.5;
     4298                                vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);
     4299                                vertices[k].onGeometry = 0;
     4300                          }
     4301                        //vertices[k].i = toI2( vertices[k].r);
     4302                        R2 AB =  vertices[k].r;
     4303                        R2 AA = (A+AB)*0.5;
     4304                        R2 BB = (AB+B)*0.5;
     4305                        vertices[k].ReferenceNumber = edges[i].ref;
     4306                        vertices[k].DirOfSearch = NoDirOfSearch;
     4307
     4308                        newedges[ie].onGeometry = Gh.Containing(AA,ong);
     4309                        newedges[ie++].v[1]=vertices+k;
     4310
     4311                        newedges[ie]=edges[i];
     4312                        newedges[ie].adj[0]=newedges + ie -1;
     4313                        newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;
     4314                        newedges[ie].onGeometry =  Gh.Containing(BB,ong);
     4315                        newedges[ie++].v[0]=vertices+k;
     4316                        k++;
     4317                  }
     4318                if (edgesGtoB) delete [] edgesGtoB;
     4319                edgesGtoB=0;
     4320
     4321                newnbv=k;
     4322                newNbVerticesOnGeomEdge=kvg;
     4323                if (newnbv> maxnbv) goto Error;// bug
     4324
     4325                nbv = k;
     4326
     4327
     4328                kedge = new long[3*nbt+1];
     4329                ksplitarray = new long[nbt+1];
     4330                ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0]
     4331
     4332                for (i=0;i<3*nbt;i++)
     4333                 kedge[i]=-1;
     4334
     4335                // 
     4336
     4337                for (i=0;i<nbt;i++) {
     4338                        Triangle & t = triangles[i];
     4339                        if (!t.link){
     4340                                ISSMERROR("!t.link");
     4341                        }
     4342                        for(int j=0;j<3;j++)
     4343                          {
     4344                                const TriangleAdjacent ta = t.Adj(j);
     4345                                const Triangle & tt = ta;
     4346                                if (&tt >= lastT)
     4347                                 t.SetAdj2(j,0,0);// unset adj
     4348                                const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
     4349                                const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
     4350                                long  ke =edge4->SortAndFind(Number(v0),Number(v1));
     4351                                if (ke>0)
     4352                                  {
     4353                                        long ii = Number(tt);
     4354                                        int  jj = ta;
     4355                                        long ks = ke + nbvold;
     4356                                        kedge[3*i+j] = ks;
     4357                                        if (ii<nbt) // good triangle
     4358                                         kedge[3*ii+jj] = ks;
     4359                                        BamgVertex &A=vertices[ks];
     4360                                        double aa,bb,cc,dd;
     4361                                        if ((dd=Area2(v0.r,v1.r,A.r)) >=0){
     4362                                                // warning PB roundoff error
     4363                                                if (t.link && ( (aa=Area2( A.r    , t[1].r , t[2].r )) < 0.0
     4364                                                                                ||   (bb=Area2( t[0].r , A.r    , t[2].r )) < 0.0 
     4365                                                                                ||   (cc=Area2( t[0].r , t[1].r , A.r    )) < 0.0)){
     4366                                                        printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd);
     4367                                                        ISSMERROR("Number of triangles with P2 interpolation Problem");
     4368                                                }
     4369                                        }
     4370                                        else {
     4371                                                if (tt.link && ( (aa=Area2( A.r     , tt[1].r , tt[2].r )) < 0
     4372                                                                                ||   (bb=Area2( tt[0].r , A.r     , tt[2].r )) < 0
     4373                                                                                ||   (cc=Area2( tt[0].r , tt[1].r , A.r     )) < 0)){
     4374                                                        printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd);
     4375                                                        ISSMERROR("Number of triangles with P2 interpolation Problem");
     4376                                                }
     4377                                        }
     4378                                  }
     4379                          }
     4380                }
     4381
     4382                for (i=0;i<nbt;i++){
     4383                        ksplit[i]=1; // no split by default
     4384                        const Triangle & t = triangles[ i];
     4385                        int nbsplitedge =0;
     4386                        int nbinvisible =0;
     4387                        int invisibleedge=0;
     4388                        int kkk[3];     
     4389                        for (int j=0;j<3;j++)
     4390                          {
     4391                                if (t.Hidden(j)) invisibleedge=j,nbinvisible++;
     4392
     4393                                const TriangleAdjacent ta = t.Adj(j);
     4394                                const Triangle & tt = ta;
     4395
     4396
     4397                                const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
     4398                                const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
     4399                                if ( kedge[3*i+j] < 0)
     4400                                  {
     4401                                        long  ke =edge4->SortAndFind(Number(v0),Number(v1));
     4402                                        if (ke<0) // new
     4403                                          {
     4404                                                if (&tt) // internal triangles all the boundary
     4405                                                  { // new internal edges
     4406                                                        long ii = Number(tt);
     4407                                                        int  jj = ta;
     4408
     4409                                                        kedge[3*i+j]=k;// save the vertex number
     4410                                                        kedge[3*ii+jj]=k;
     4411                                                        if (k<maxnbv)
     4412                                                          {
     4413                                                                vertices[k].r = ((R2) v0+(R2) v1 )/2;
     4414                                                                //vertices[k].i = toI2( vertices[k].r);
     4415                                                                vertices[k].ReferenceNumber=0;
     4416                                                                vertices[k].DirOfSearch =NoDirOfSearch;
     4417                                                                vertices[k].m =  Metric(0.5,v0,0.5,v1);
     4418                                                          }
     4419                                                        k++;
     4420                                                        kkk[nbsplitedge++]=j;                 
     4421                                                  } // tt
     4422                                                else
     4423                                                 ISSMERROR("Bug...");
     4424                                          } // ke<0           
     4425                                        else
     4426                                          { // ke >=0
     4427                                                kedge[3*i+j]=nbvold+ke;
     4428                                                kkk[nbsplitedge++]=j;// previously splited
     4429                                          }
     4430                                  }
     4431                                else
     4432                                 kkk[nbsplitedge++]=j;// previously splited
     4433
     4434                          }
     4435                        if (nbinvisible>=2){
     4436                                ISSMERROR("nbinvisible>=2");
     4437                        }
     4438                        switch (nbsplitedge) {
     4439                                case 0: ksplit[i]=10; newnbt++; break;   // nosplit
     4440                                case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2
     4441                                case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3
     4442                                case 3:
     4443                                                  if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4;
     4444                                                  else   ksplit[i]=10*nfortria,newnbt+=nfortria;
     4445                                                  break;
     4446                        }
     4447                        if (ksplit[i]<40){
     4448                                ISSMERROR("ksplit[i]<40");
     4449                        }
     4450                  }
     4451                //  now do the element split
     4452                newNbOfQuad = 4*NbOfQuad;
     4453                nbv = k;
     4454                kkk = nbt;
     4455                ksplit[-1] = nbt;
     4456                // look on  old true  triangles
     4457
     4458                for (i=0;i<nbtsave;i++){
     4459                        int  nbmkadj=0;
     4460                        long mkadj [100];
     4461                        mkadj[0]=i;
     4462                        long kk=ksplit[i]/10;
     4463                        int  ke=(int) (ksplit[i]%10);
     4464                        if (kk>=7 || kk<=0){
     4465                                ISSMERROR("kk>=7 || kk<=0");
     4466                        }
     4467
     4468                        // def the numbering   k (edge) i vertex
     4469                        int k0 = ke;
     4470                        int k1 = NextEdge[k0];
     4471                        int k2 = PreviousEdge[k0];
     4472                        int i0 = OppositeVertex[k0];
     4473                        int i1 = OppositeVertex[k1];
     4474                        int i2 = OppositeVertex[k2];
     4475
     4476                        Triangle &t0=triangles[i];
     4477                        BamgVertex * v0=t0(i0);           
     4478                        BamgVertex * v1=t0(i1);           
     4479                        BamgVertex * v2=t0(i2);
     4480
     4481                        if (nbmkadj>=10){
     4482                                ISSMERROR("nbmkadj>=10");
     4483                        }
     4484                        // --------------------------
     4485                        TriangleAdjacent ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2));
     4486                        // save the flag Hidden
     4487                        int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)};
     4488                        // un set all adj -- save Hidden flag --
     4489                        t0.SetAdj2(0,0,hid[0]);
     4490                        t0.SetAdj2(1,0,hid[1]);
     4491                        t0.SetAdj2(2,0,hid[2]);
     4492                        // --  remake
     4493                        switch  (kk) {
     4494                                case 1: break;// nothing
     4495                                case 2: //
     4496                                                  {
     4497                                                        Triangle &t1=triangles[kkk++];
     4498                                                        t1=t0;
     4499                                                        if (kedge[3*i+i0]<0){
     4500                                                                ISSMERROR("kedge[3*i+i0]<0");
     4501                                                        }
     4502                                                        BamgVertex * v3 = vertices + kedge[3*i+k0];
     4503
     4504                                                        t0(i2) = v3;
     4505                                                        t1(i1) = v3;
     4506                                                        t0.SetAllFlag(k2,0);
     4507                                                        t1.SetAllFlag(k1,0);
     4508                                                  }
     4509                                                break;
     4510                                case 3: //
     4511                                                  {
     4512                                                        Triangle &t1=triangles[kkk++];
     4513                                                        Triangle &t2=triangles[kkk++];
     4514                                                        t2=t1=t0;
     4515                                                        if (kedge[3*i+k1]<0){
     4516                                                                ISSMERROR("kedge[3*i+k1]<0");
     4517                                                        }
     4518                                                        if (kedge[3*i+k2]<0){
     4519                                                                ISSMERROR("kedge[3*i+k2]<0");
     4520                                                        }
     4521
     4522                                                        BamgVertex * v01 = vertices + kedge[3*i+k2];
     4523                                                        BamgVertex * v02 = vertices + kedge[3*i+k1];
     4524                                                        t0(i1) = v01;
     4525                                                        t0(i2) = v02;
     4526                                                        t1(i2) = v02;
     4527                                                        t1(i0) = v01;
     4528                                                        t2(i0) = v02;
     4529                                                        t0.SetAllFlag(k0,0);
     4530                                                        t1.SetAllFlag(k1,0);
     4531                                                        t1.SetAllFlag(k0,0);
     4532                                                        t2.SetAllFlag(k2,0);
     4533                                                  }
     4534                                                break;
     4535                                case 4: //
     4536                                case 6: // split in 4
     4537                                                  {
     4538                                                        Triangle &t1=triangles[kkk++];
     4539                                                        Triangle &t2=triangles[kkk++];
     4540                                                        Triangle &t3=triangles[kkk++];
     4541                                                        t3=t2=t1=t0;
     4542                                                        if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){
     4543                                                                ISSMERROR("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0");
     4544                                                        }
     4545                                                        BamgVertex * v12 = vertices + kedge[3*i+k0];
     4546                                                        BamgVertex * v02 = vertices + kedge[3*i+k1];
     4547                                                        BamgVertex * v01 = vertices + kedge[3*i+k2];
     4548                                                        t0(i1) = v01;
     4549                                                        t0(i2) = v02;
     4550                                                        t0.SetAllFlag(k0,hid[k0]);
     4551
     4552                                                        t1(i0) = v01;
     4553                                                        t1(i2) = v12;
     4554                                                        t0.SetAllFlag(k1,hid[k1]);
     4555
     4556                                                        t2(i0) = v02;
     4557                                                        t2(i1) = v12;
     4558                                                        t2.SetAllFlag(k2,hid[k2]);
     4559
     4560                                                        t3(i0) = v12;
     4561                                                        t3(i1) = v02;
     4562                                                        t3(i2) = v01;
     4563
     4564                                                        t3.SetAllFlag(0,hid[0]);           
     4565                                                        t3.SetAllFlag(1,hid[1]);           
     4566                                                        t3.SetAllFlag(2,hid[2]);
     4567
     4568                                                        if ( kk == 6)
     4569                                                          {
     4570
     4571                                                                Triangle &t4=triangles[kkk++];
     4572                                                                Triangle &t5=triangles[kkk++];
     4573
     4574                                                                t4 = t3;
     4575                                                                t5 = t3;
     4576
     4577                                                                t0.SetHidden(k0);
     4578                                                                t1.SetHidden(k1);
     4579                                                                t2.SetHidden(k2);
     4580                                                                t3.SetHidden(0);
     4581                                                                t4.SetHidden(1);
     4582                                                                t5.SetHidden(2);
     4583
     4584                                                                if (nbv < maxnbv )
     4585                                                                  {
     4586                                                                        vertices[nbv].r = ((R2) *v01 + (R2) *v12  + (R2) *v02 ) / 3.0;
     4587                                                                        vertices[nbv].ReferenceNumber =0;
     4588                                                                        vertices[nbv].DirOfSearch =NoDirOfSearch;
     4589                                                                        //vertices[nbv].i = toI2(vertices[nbv].r);
     4590                                                                        double a3[]={1./3.,1./3.,1./3.};
     4591                                                                        vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m);
     4592                                                                        BamgVertex * vc =  vertices +nbv++;
     4593                                                                        t3(i0) = vc;
     4594                                                                        t4(i1) = vc;
     4595                                                                        t5(i2) = vc;
     4596
     4597                                                                  }
     4598                                                                else
     4599                                                                 goto Error;
     4600                                                          }
     4601
     4602                                                  }
     4603                                                break;         
     4604                        }
     4605
     4606                        // save all the new triangles
     4607                        mkadj[nbmkadj++]=i;
     4608                        long jj;
     4609                        if (t0.link)
     4610                         for (jj=nbt;jj<kkk;jj++)
     4611                                {
     4612                                 triangles[jj].link=t0.link;
     4613                                 t0.link= triangles+jj;
     4614                                 mkadj[nbmkadj++]=jj;
     4615                                }
     4616                        if (nbmkadj>13){// 13 = 6 + 4 +
     4617                                ISSMERROR("nbmkadj>13");
     4618                        }
     4619
     4620                        if (kk==6)  newNbOfQuad+=3;
     4621                        for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;
     4622                        ksplit[i]= nbt; // save last adresse of the new triangles
     4623                        kkk = nbt;
     4624                  }
     4625
     4626                for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.;
     4627
     4628                if(withBackground)
     4629                 for (i=0;i<BTh.nbv;i++)
     4630                  BTh.vertices[i].m =  BTh.vertices[i].m*2.;
     4631
     4632
     4633                ret = 2;
     4634                if (nbt>= nbtx) goto Error; // bug
     4635                if (nbv>= maxnbv) goto Error; // bug
     4636                // generation of the new triangles
     4637
     4638                SetIntCoor("In SplitElement");
     4639
     4640                ReMakeTriangleContainingTheVertex();
     4641                if(withBackground)
     4642                 BTh.ReMakeTriangleContainingTheVertex();
     4643
     4644                delete [] edges;
     4645                edges = newedges;
     4646                nbe = newnbe;
     4647                NbOfQuad = newNbOfQuad;
     4648
     4649                for (i=0;i<NbSubDomains;i++)
     4650                  {
     4651                        long k = subdomains[i].edge- edges;
     4652                        subdomains[i].edge =  edges+2*k; // spilt all edge in 2
     4653                  }
     4654
     4655                if (ksplitarray) delete [] ksplitarray;
     4656                if (kedge) delete [] kedge;
     4657                if (edge4) delete edge4;
     4658                if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
     4659                VerticesOnGeomEdge= newVerticesOnGeomEdge;
     4660                if(VertexOnBThEdge) delete []  VertexOnBThEdge;
     4661                VertexOnBThEdge = newVertexOnBThEdge;
     4662                NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
     4663                NbVertexOnBThEdge=newNbVertexOnBThEdge;
     4664                //  ReMakeTriangleContainingTheVertex();
     4665
     4666                ReconstructExistingMesh();
     4667
     4668                if (verbose>2){
     4669                        printf("   number of quadrilaterals    = %i\n",NbOfQuad);
     4670                        printf("   number of triangles         = %i\n",nbt-NbOutT- NbOfQuad*2);
     4671                        printf("   number of outside triangles = %i\n",NbOutT);
     4672                }
     4673
     4674                return 0; //ok
     4675
     4676Error:
     4677                nbv = nbvold;
     4678                nbt = nbtold;
     4679                NbOutT = NbOutTold;
     4680                // cleaning memory ---
     4681                delete newedges;
     4682                if (ksplitarray) delete [] ksplitarray;
     4683                if (kedge) delete [] kedge;
     4684                if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge;
     4685                if (edge4) delete edge4;
     4686                if(newVertexOnBThEdge) delete []  newVertexOnBThEdge;
     4687
     4688                return ret; // ok
     4689        }
     4690        /*}}}1*/
     4691/*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{1*/
     4692long  Mesh::SplitInternalEdgeWithBorderVertices(){
     4693        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)*/
     4694
     4695        long NbSplitEdge=0;
     4696        SetVertexFieldOn(); 
     4697        long it;
     4698        long nbvold=nbv;
     4699        long int verbose=2;
     4700        for (it=0;it<nbt;it++){
     4701                Triangle &t=triangles[it];
     4702                if (t.link)
     4703                 for (int j=0;j<3;j++)
     4704                  if(!t.Locked(j) && !t.Hidden(j)){
     4705                          Triangle &tt = *t.TriangleAdj(j);
     4706                          if ( &tt && tt.link && it < Number(tt))
     4707                                 { // an internal edge
     4708                                  BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
     4709                                  BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
     4710                                  if (v0.onGeometry && v1.onGeometry){
     4711                                          R2 P= ((R2) v0 + (R2) v1)*0.5;
     4712                                          if ( nbv<maxnbv) {
     4713                                                  vertices[nbv].r = P;
     4714                                                  vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);
     4715                                                  vertices[nbv].ReferenceNumber=0;
     4716                                                  vertices[nbv].DirOfSearch = NoDirOfSearch ;
     4717                                          }
     4718                                          NbSplitEdge++;
     4719                                  }
     4720                                 }
     4721                  }
     4722        }
     4723        ReMakeTriangleContainingTheVertex();   
     4724        if (nbvold!=nbv){
     4725                long  iv = nbvold;
     4726                long NbSwap = 0;
     4727                Icoor2 dete[3]; 
     4728                for (int i=nbvold;i<nbv;i++) {// for all the new point
     4729                        BamgVertex & vi = vertices[i];
     4730                        vi.i = toI2(vi.r);
     4731                        vi.r = toR2(vi.i);
     4732
     4733                        // a good new point
     4734                        vi.ReferenceNumber=0;
     4735                        vi.DirOfSearch =NoDirOfSearch;
     4736                        Triangle *tcvi = FindTriangleContaining(vi.i,dete);
     4737                        if (tcvi && !tcvi->link) {
     4738                                printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
     4739                        }
     4740
     4741                        quadtree->Add(vi);
     4742                        if (!tcvi || tcvi->det<0){// internal
     4743                                ISSMERROR("!tcvi || tcvi->det < 0");
     4744                        }
     4745                        AddVertex(vi,tcvi,dete);
     4746                        NbSwap += vi.Optim(1);         
     4747                        iv++;
     4748                        //      }
     4749        }
     4750        if (verbose>3) {
     4751                printf("   number of points: %i\n",iv);
     4752                printf("   number of swap to  split internal edges with border vertices: %i\n",NbSwap);
     4753                nbv = iv;
     4754        }
     4755}
     4756if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices  to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold));
     4757if (verbose>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge);
     4758
     4759return  NbSplitEdge;
     4760}
     4761/*}}}1*/
     4762        /*FUNCTION Mesh::TriangulateFromGeom0{{{1*/
     4763        void Mesh::TriangulateFromGeom0(long imaxnbv,BamgOpts* bamgopts){
    27264764                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles0)*/
    27274765
    … …  
    27374775                R2 AB;
    27384776                GeometricalVertex *a,*b;
    2739                 MeshVertex *va,*vb;
     4777                BamgVertex *va,*vb;
    27404778                GeometricalEdge *e;
    27414779
    … …  
    27454783                //initialize this
    27464784                if (verbose>3) printf("      Generating Boundary vertices\n");
    2747                 PreInit(inbvx);
     4785                Init(imaxnbv);
    27484786                nbv=0;
    27494787                NbVerticesOnGeomVertex=0;
    … …  
    27604798                //initialize VerticesOnGeomVertex
    27614799                VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 
    2762                 if( NbVerticesOnGeomVertex >= nbvx) {
    2763                         ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx);
     4800                if( NbVerticesOnGeomVertex >= maxnbv) {
     4801                        ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);
    27644802                }
    27654803
    … …  
    27734811                        if (Gh[i].Required() && Gh[i].IsThe()) {//Gh  vertices Required
    27744812
    2775                                 //Add the vertex (provided that nbv<nbvx)
    2776                                 if (nbv<nbvx){
     4813                                //Add the vertex (provided that nbv<maxnbv)
     4814                                if (nbv<maxnbv){
    27774815                                        vertices[nbv]=Gh[i];
    27784816                                }
    27794817                                else{
    2780                                         ISSMERROR("Maximum number of vertices (nbvx = %i) too small",nbvx);
     4818                                        ISSMERROR("Maximum number of vertices (maxnbv = %i) too small",maxnbv);
    27814819                                }
    27824820                               
    … …  
    28604898                                                                NbNewPoints=0;
    28614899                                                                NbEdgeCurve=0;
    2862                                                                 if (nbvend>=nbvx) ISSMERROR("maximum number of vertices too low! Check the domain outline or increase nbvx");
     4900                                                                if (nbvend>=maxnbv) ISSMERROR("maximum number of vertices too low! Check the domain outline or increase maxnbv");
    28634901                                                                lcurve =0;
    28644902                                                                s = lstep;
    … …  
    30335071        }
    30345072        /*}}}1*/
    3035         /*FUNCTION Mesh::GeomToTriangles1{{{1*/
    3036         void Mesh::GeomToTriangles1(long inbvx,BamgOpts* bamgopts,int KeepVertices){
     5073        /*FUNCTION Mesh::TriangulateFromGeom1{{{1*/
     5074        void Mesh::TriangulateFromGeom1(long imaxnbv,BamgOpts* bamgopts,int KeepVertices){
    30375075                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles1)*/
    30385076
    … …  
    30645102
    30655103                //Initialize new mesh
    3066                 this->PreInit(inbvx);
     5104                this->Init(imaxnbv);
    30675105                BTh.SetVertexFieldOn();
    30685106                int* bcurve = new int[Gh.NbOfCurves]; //
    … …  
    30825120                int i;
    30835121                for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
    3084                 if(NbVerticesOnGeomVertex >= nbvx) { ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx);}
     5122                if(NbVerticesOnGeomVertex >= maxnbv) { ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);}
    30855123
    30865124                VerticesOnGeomVertex = new VertexOnGeom[  NbVerticesOnGeomVertex];
    30875125                VertexOnBThVertex    = new VertexOnVertex[NbVerticesOnGeomVertex];
    30885126
    3089                 //At this point there is NO vertex but vertices should have been allocated by PreInit
     5127                //At this point there is NO vertex but vertices should have been allocated by Init
    30905128                ISSMASSERT(vertices);
    30915129                for (i=0;i<Gh.nbv;i++){
    … …  
    31035141                        if (vog.IsRequiredVertex()){
    31045142                                GeometricalVertex* gv=vog;
    3105                                 MeshVertex *bv = vog;
     5143                                BamgVertex *bv = vog;
    31065144                                ISSMASSERT(gv->to); // use of Geom -> Th
    31075145                                VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->to,bv);
    … …  
    32085246                                                long i=0;// index of new points on the curve
    32095247                                                register GeometricalVertex * GA0 = *(*peequi)[k0equi].onGeometry;
    3210                                                 MeshVertex *A0;
     5248                                                BamgVertex *A0;
    32115249                                                A0 = GA0->to;  // the vertex in new mesh
    3212                                                 MeshVertex *A1;
     5250                                                BamgVertex *A1;
    32135251                                                VertexOnGeom *GA1;
    32145252                                                Edge* PreviousNewEdge = 0;
    … …  
    32285266                                                                ISSMASSERT(pe && ee.onGeometry);
    32295267                                                                ee.onGeometry->SetMark();
    3230                                                                 MeshVertex & v0=ee[0], & v1=ee[1];
     5268                                                                BamgVertex & v0=ee[0], & v1=ee[1];
    32315269                                                                R2 AB=(R2)v1-(R2)v0;
    32325270                                                                double L0=L,LAB;
    … …  
    32415279                                                                                ISSMASSERT(sNew>=L0);
    32425280                                                                                ISSMASSERT(LAB);
    3243                                                                                 ISSMASSERT(vertices && nbv<nbvx);
     5281                                                                                ISSMASSERT(vertices && nbv<maxnbv);
    32445282                                                                                ISSMASSERT(edges && nbe<nbex);
    32455283                                                                                ISSMASSERT(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
    … …  
    33285366                        //Allocate memory
    33295367                        if(step==0){
    3330                                 if(nbv+NbOfNewPoints > nbvx) {
    3331                                         ISSMERROR("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,nbvx);
     5368                                if(nbv+NbOfNewPoints > maxnbv) {
     5369                                        ISSMERROR("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,maxnbv);
    33325370                                }
    33335371                                edges = new Edge[NbOfNewEdge];
    … …  
    33655403        }
    33665404        /*}}}1*/
    3367         /*FUNCTION Mesh::Insert{{{1*/
    3368         void Mesh::Insert() {
    3369                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/
    3370 
    3371                 /*Insert points in the existing Geometry*/
    3372 
    3373                 //Intermediary
    3374                 int i;
    3375 
    3376                 /*Get options*/
    3377                 long int verbose=2;
    3378 
    3379                 //Display info
    3380                 if (verbose>2) printf("   Insert initial %i vertices\n",nbv);
    3381 
    3382                 //Compute integer coordinates and determinants for the existing vertices (from Geometry)
    3383                 SetIntCoor();
    3384 
    3385                 /*Now we want to build a list (ordre) of the vertices in a random
    3386                  * order. To do so, we use the following method:
    3387                  *
    3388                  * From an initial k0 in [0 nbv[ random (vertex number)
    3389                  * the next k (vertex number) is computed using a big
    3390                  * prime number (PN>>nbv) following:
    3391                  *
    3392                  * k_{i+1} = k_i + PN  [nbv]
    3393                  *
    3394                  * let's show that:
    3395                  *
    3396                  *   for all j in [0 nbv[, ∃! i in [0 nbv[ such that k_i=j
    3397                  *
    3398                  * Let's assume that there are 2 distinct j1 and j2 such that
    3399                  * k_j1 = k_j2
    3400                  *
    3401                  * This means that
    3402                  * 
    3403                  *  k0+j1*PN = k0+j2*PN [nbv]
    3404                  *  (j1-j2)*PN =0       [nbv]
    3405                  * since PN is a prime number larger than nbv, and nbv!=1
    3406                  *  j1-j2=0             [nbv]
    3407                  * BUT
    3408                  *  j1 and j2 are in [0 nbv[ which is impossible.
    3409                  *
    3410                  *  We hence have built a random list of nbv elements of
    3411                  *  [0 nbv[ all distincts*/
    3412                 for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ;
    3413                 const long PrimeNumber= BigPrimeNumber(nbv) ;
    3414                 int   k0=rand()%nbv;
    3415                 for (int is3=0; is3<nbv; is3++){
    3416                         ordre[is3]= &vertices[k0=(k0+PrimeNumber)%nbv];
    3417                 }
    3418 
    3419                 /*Modify ordre such that the first 3 vertices form a triangle*/
    3420 
    3421                 //get first vertex i such that [0,1,i] are not aligned
    3422                 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){
    3423                         //if i is higher than nbv, it means that all the determinants are 0,
    3424                         //all vertices are aligned!
    3425                         if  ( ++i >= nbv) {
    3426                                 ISSMERROR("all the vertices are aligned");
    3427                         }
    3428                 }
    3429                 // exchange i et 2 in "ordre" so that
    3430                 // the first 3 vertices are not aligned (real triangle)
    3431                 Exchange(ordre[2], ordre[i]);
    3432 
    3433                 /*Take the first edge formed by the first two vertices and build
    3434                  * the initial simple mesh from this edge and 2 boundary triangles*/
    3435 
    3436                 MeshVertex *  v0=ordre[0], *v1=ordre[1];
    3437 
    3438                 nbt = 2;
    3439                 triangles[0](0) = NULL; //infinite vertex
    3440                 triangles[0](1) = v0;
    3441                 triangles[0](2) = v1;
    3442                 triangles[1](0) = NULL;//infinite vertex
    3443                 triangles[1](2) = v0;
    3444                 triangles[1](1) = v1;
    3445 
    3446                 //Build adjacence
    3447                 const int e0 = OppositeEdge[0];
    3448                 const int e1 = NextEdge[e0];
    3449                 const int e2 = PreviousEdge[e0];
    3450                 triangles[0].SetAdj2(e0, &triangles[1] ,e0);
    3451                 triangles[0].SetAdj2(e1, &triangles[1] ,e2);
    3452                 triangles[0].SetAdj2(e2, &triangles[1] ,e1);
    3453 
    3454                 triangles[0].det = -1;  //boundary triangle: det = -1
    3455                 triangles[1].det = -1;  //boundary triangle: det = -1
    3456 
    3457                 triangles[0].SetTriangleContainingTheVertex();
    3458                 triangles[1].SetTriangleContainingTheVertex();
    3459 
    3460                 triangles[0].link=&triangles[1];
    3461                 triangles[1].link=&triangles[0];
    3462 
    3463                 //build quadtree
    3464                 if (!quadtree)  quadtree = new QuadTree(this,0);
    3465                 quadtree->Add(*v0);
    3466                 quadtree->Add(*v1);
    3467 
    3468                 /*Now, add the vertices One by One*/
    3469 
    3470                 long NbSwap=0;
    3471                 if (verbose>3) printf("   Begining of insertion process...\n");
    3472 
    3473                 for (int icount=2; icount<nbv; icount++) {
    3474 
    3475                         //Get new vertex
    3476                         MeshVertex *newvertex=ordre[icount];
    3477 
    3478                         //Find the triangle in which newvertex is located
    3479                         Icoor2 dete[3];
    3480                         Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates)
    3481 
    3482                         //Add newvertex to the quadtree
    3483                         quadtree->Add(*newvertex);
    3484 
    3485                         //Add newvertex to the existing mesh
    3486                         AddVertex(*newvertex,tcvi,dete);
    3487 
    3488                         //Make the mesh Delaunay around newvertex by swaping the edges
    3489                         NbSwap += newvertex->Optim(1,0);
    3490                 }
    3491 
    3492                 //Display info
    3493                 if (verbose>3) {
    3494                         printf("      NbSwap of insertion: %i\n",NbSwap);
    3495                         printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
    3496                 }
    3497 
    3498 #ifdef NBLOOPOPTIM
    3499 
    3500                 k0 = rand()%nbv ;
    3501                 for (int is4=0; is4<nbv; is4++)
    3502                  ordre[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv];
    3503 
    3504                 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){
    3505                         long  NbSwap = 0;
    3506                         for (int is1=0; is1<nbv; is1++)
    3507                          NbSwap += ordre[is1]->Optim(0,0);
    3508                         if (verbose>3) {
    3509                                 printf("      Optim Loop: %i\n",Nbloop);
    3510                                 printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
    3511                         }
    3512                         if(!NbSwap) break;
    3513                 }
    3514                 ReMakeTriangleContainingTheVertex();
    3515                 // because we break the TriangleContainingTheVertex
    3516 #endif
    3517         }
    3518         /*}}}1*/
    3519         /*FUNCTION Mesh::InsertNewPoints{{{1*/
    3520         long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {
    3521                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/
    3522 
    3523                 long int verbose=0;
    3524                 double seuil= 1.414/2 ;// for two close point
    3525                 long i;
    3526                 long NbSwap=0;
    3527                 Icoor2 dete[3];
    3528 
    3529                 //number of new points
    3530                 const long nbvnew=nbv-nbvold;
    3531 
    3532                 //display info if required
    3533                 if (verbose>5) printf("      Try to Insert %i new points\n",nbvnew);
    3534 
    3535                 //return if no new points
    3536                 if (!nbvnew) return 0;
    3537 
    3538                 /*construction of a random order*/
    3539                 const long PrimeNumber= BigPrimeNumber(nbv)  ;
    3540                 //remainder of the division of rand() by nbvnew
    3541                 long k3 = rand()%nbvnew;
    3542                 //loop over the new points
    3543                 for (int is3=0; is3<nbvnew; is3++){
    3544                         register long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
    3545                         register long i=nbvold+is3;
    3546                         ordre[i]= vertices + j;
    3547                         ordre[i]->ReferenceNumber=i;
    3548                 }
    3549 
    3550                 // for all the new point
    3551                 long iv=nbvold;
    3552                 for (i=nbvold;i<nbv;i++){
    3553                         MeshVertex &vi=*ordre[i];
    3554                         vi.i=toI2(vi.r);
    3555                         vi.r=toR2(vi.i);
    3556                         double hx,hy;
    3557                         vi.m.Box(hx,hy);
    3558                         Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);
    3559                         if (!quadtree->ToClose(vi,seuil,hi,hj)){
    3560                                 // a good new point
    3561                                 MeshVertex &vj = vertices[iv];
    3562                                 long  j=vj.ReferenceNumber;
    3563                                 if (&vj!=ordre[j]){
    3564                                         ISSMERROR("&vj!= ordre[j]");
    3565                                 }
    3566                                 if(i!=j){
    3567                                         Exchange(vi,vj);
    3568                                         Exchange(ordre[j],ordre[i]);
    3569                                 }
    3570                                 vj.ReferenceNumber=0;
    3571                                 Triangle *tcvj=FindTriangleContaining(vj.i,dete);
    3572                                 if (tcvj && !tcvj->link){
    3573                                         tcvj->Echo();
    3574                                         ISSMERROR("problem inserting point in InsertNewPoints (tcvj=%p and tcvj->link=%i)",tcvj,tcvj->link);
    3575                                 }
    3576                                 quadtree->Add(vj);
    3577                                 AddVertex(vj,tcvj,dete);
    3578                                 NbSwap += vj.Optim(1);         
    3579                                 iv++;
    3580                         }
    3581                 }
    3582                 if (verbose>3) {
    3583                         printf("         number of new points: %i\n",iv);
    3584                         printf("         number of to close (?) points: %i\n",nbv-iv);
    3585                         printf("         number of swap after: %i\n",NbSwap);
    3586                 }
    3587                 nbv = iv;
    3588 
    3589                 for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1); 
    3590                 if (verbose>3) printf("   NbSwap=%i\n",NbSwap);
    3591 
    3592                 NbTSwap +=  NbSwap ;
    3593                 return nbv-nbvold;
    3594         }
    3595         /*}}}1*/
    3596         /*FUNCTION Mesh::MakeGeometricalEdgeToEdge{{{1*/
    3597         Edge** Mesh::MakeGeometricalEdgeToEdge() {
    3598                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeometricalEdgeToEdge)*/
    3599 
    3600                 if (!Gh.nbe){
    3601                         ISSMERROR("!Gh.nbe");
    3602                 }
    3603                 Edge **e= new (Edge* [Gh.nbe]);
    3604 
    3605                 long i;
    3606                 for ( i=0;i<Gh.nbe ; i++)
    3607                  e[i]=NULL;
    3608                 for ( i=0;i<nbe ; i++)
    3609                   {
    3610                         Edge * ei = edges+i;
    3611                         GeometricalEdge *onGeometry = ei->onGeometry;
    3612                         e[Gh.Number(onGeometry)] = ei;   
    3613                   }
    3614                 for ( i=0;i<nbe ; i++)
    3615                  for (int ii=0;ii<2;ii++) {
    3616                          Edge * ei = edges+i;
    3617                          GeometricalEdge *onGeometry = ei->onGeometry;
    3618                          int j= ii;
    3619                          while (!(*onGeometry)[j].Required()) {
    3620                                  Adj(onGeometry,j); // next geom edge
    3621                                  j=1-j;
    3622                                  if (e[Gh.Number(onGeometry)])  break; // optimisation
    3623                                  e[Gh.Number(onGeometry)] = ei;
    3624                          }
    3625                  }
    3626 
    3627                 int kk=0;
    3628                 for ( i=0;i<Gh.nbe ; i++){
    3629                         if (!e[i]){
    3630                                 kk++;
    3631                                 if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i);
    3632                         }
    3633                 }
    3634                 if(kk) ISSMERROR("See above");
    3635 
    3636                 return e;
    3637         }
    3638         /*}}}1*/
    3639         /*FUNCTION Mesh::MakeQuadrangles{{{1*/
    3640         void Mesh::MakeQuadrangles(double costheta){
    3641                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadrangles)*/
    3642 
    3643                 long int verbose=0;
    3644 
    3645                 if (verbose>2) printf("MakeQuadrangles costheta = %g\n",costheta);
    3646 
    3647                 if (costheta >1) {
    3648                         if (verbose>5) printf("   do nothing: costheta > 1\n");
    3649                 }
    3650 
    3651 
    3652                         long nbqq = (nbt*3)/2;
    3653                         DoubleAndInt *qq = new DoubleAndInt[nbqq];
    3654 
    3655                         long i,ij;
    3656                         int j;
    3657                         long k=0;
    3658                         for (i=0;i<nbt;i++)
    3659                          for (j=0;j<3;j++)
    3660                           if ((qq[k].q=triangles[i].QualityQuad(j))>=costheta)
    3661                                 qq[k++].i3j=i*3+j;
    3662                         //  sort  qq
    3663                         HeapSort(qq,k);
    3664 
    3665                         long kk=0;
    3666                         for (ij=0;ij<k;ij++) {
    3667                                 i=qq[ij].i3j/3;
    3668                                 j=(int) (qq[ij].i3j%3);
    3669                                 // optisamition no float computation 
    3670                                 if (triangles[i].QualityQuad(j,0) >=costheta)
    3671                                  triangles[i].SetHidden(j),kk++;
    3672                           }
    3673                         NbOfQuad = kk;
    3674                         if (verbose>2){
    3675                                 printf("   number of quadrilaterals    = %i\n",NbOfQuad);
    3676                                 printf("   number of triangles         = %i\n",nbt-NbOutT- NbOfQuad*2);
    3677                                 printf("   number of outside triangles = %i\n",NbOutT);
    3678                         }
    3679                         delete [] qq;
    3680         }
    3681         /*}}}1*/
    3682         /*FUNCTION Mesh::MakeQuadTree{{{1*/
    3683         void Mesh::MakeQuadTree() { 
    3684                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadTree)*/
    3685 
    3686                 long int verbose=0;
    3687                 if (  !quadtree )  quadtree = new QuadTree(this);
    3688 
    3689         }
    3690         /*}}}1*/
    3691         /*FUNCTION Mesh::MaxSubDivision{{{1*/
    3692         void  Mesh::MaxSubDivision(double maxsubdiv) {
    3693                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/
    3694 
    3695                 long int verbose=0;
    3696 
    3697                 const  double maxsubdiv2 = maxsubdiv*maxsubdiv;
    3698                 if(verbose>1) printf("   Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv);
    3699                 // for all the edges
    3700                 // if the len of the edge is to long
    3701                 long it,nbchange=0;   
    3702                 double lmax=0;
    3703                 for (it=0;it<nbt;it++){
    3704                         Triangle &t=triangles[it];
    3705                         for (int j=0;j<3;j++){
    3706                                 Triangle &tt = *t.TriangleAdj(j);
    3707                                 if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)){
    3708                                         MeshVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
    3709                                         MeshVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
    3710                                         R2 AB= (R2) v1-(R2) v0;
    3711                                         Metric M = v0;
    3712                                         double l = M(AB,AB);
    3713                                         lmax = Max(lmax,l);
    3714                                         if(l> maxsubdiv2){
    3715                                           R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
    3716                                                 double lc = M(AC,AC);
    3717                                                 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
    3718                                                 D2xD2 Rt1(Rt.inv());
    3719                                                 D2xD2 D(maxsubdiv2,0,0,lc);
    3720                                                 D2xD2 MM = Rt1*D*Rt1.t();
    3721                                                 v0.m =  M = Metric(MM.x.x,MM.y.x,MM.y.y);
    3722                                                 nbchange++;
    3723                                         }
    3724                                         M = v1;
    3725                                         l = M(AB,AB);
    3726                                         lmax = Max(lmax,l);
    3727                                         if(l> maxsubdiv2){
    3728                                           R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
    3729                                                 double lc = M(AC,AC);
    3730                                                 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
    3731                                                 D2xD2 Rt1(Rt.inv());
    3732                                                 D2xD2 D(maxsubdiv2,0,0,lc);
    3733                                                 D2xD2  MM = Rt1*D*Rt1.t();
    3734                                                 v1.m =  M = Metric(MM.x.x,MM.y.x,MM.y.y);
    3735                                                 nbchange++;
    3736                                         }
    3737                                 }
    3738                         }
    3739                 }
    3740                 if(verbose>3){
    3741                         printf("      number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5));
    3742                 }
    3743         }
    3744         /*}}}1*/
    3745         /*FUNCTION Mesh::MetricAt{{{1*/
    3746         Metric Mesh::MetricAt(const R2 & A) const {
    3747                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/
    3748 
    3749                 I2 a = toI2(A);
    3750                 Icoor2 deta[3];
    3751                 Triangle * t =FindTriangleContaining(a,deta);
    3752                 if (t->det <0) { // outside
    3753                         double ba,bb;
    3754                         TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ;
    3755                         return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}
    3756                 else { // inside
    3757                         double   aa[3];
    3758                         double s = deta[0]+deta[1]+deta[2];
    3759                         aa[0]=deta[0]/s;
    3760                         aa[1]=deta[1]/s;
    3761                         aa[2]=deta[2]/s;
    3762                         return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);
    3763                 }
    3764         }
    3765         /*}}}1*/
    3766 /*FUNCTION Mesh::NearestVertex{{{1*/
    3767 MeshVertex* Mesh::NearestVertex(Icoor1 i,Icoor1 j) {
    3768         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NearestVertex)*/
    3769         return  quadtree->NearestVertex(i,j);
    3770 }
    3771 /*}}}1*/
    3772         /*FUNCTION Mesh::NewPoints{{{1*/
    3773         void  Mesh::NewPoints(Mesh & Bh,BamgOpts* bamgopts,int KeepVertices){
    3774                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
    3775 
    3776                 int i,j,k;
    3777                 long NbTSwap=0;
    3778                 long nbtold=nbt;
    3779                 long nbvold=nbv;
    3780                 long Headt=0;
    3781                 long next_t;
    3782                 long* first_np_or_next_t=new long[nbtx];
    3783                 Triangle* t=NULL;
    3784 
    3785                 /*Recover options*/
    3786                 int verbose=bamgopts->verbose;
    3787 
    3788                 /*First, insert old points if requested*/
    3789                 if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< nbvx)){
    3790                         if (verbose>5) printf("         Inserting initial mesh points\n");
    3791                         for (i=0;i<Bh.nbv;i++){
    3792                                 MeshVertex &bv=Bh[i];
    3793                                 if (!bv.onGeometry){
    3794                                         vertices[nbv].r   = bv.r;
    3795                                         vertices[nbv++].m = bv.m;
    3796                                 }
    3797                         }
    3798                         Bh.ReMakeTriangleContainingTheVertex();     
    3799                         InsertNewPoints(nbvold,NbTSwap)   ;           
    3800                 } 
    3801                 else Bh.ReMakeTriangleContainingTheVertex();     
    3802 
    3803                 // generation of the list of next Triangle
    3804                 for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);
    3805                 // the next traingle of i is -first_np_or_next_t[i]
    3806 
    3807                 // Big loop (most time consuming)
    3808                 int iter=0;
    3809                 if (verbose>5) printf("         Big loop\n");
    3810                 do {
    3811                         /*Update variables*/
    3812                         iter++;
    3813                         nbtold=nbt;
    3814                         nbvold=nbv;
    3815 
    3816                         /*We test all triangles*/
    3817                         i=Headt;
    3818                         next_t=-first_np_or_next_t[i];
    3819                         for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){
    3820 
    3821                                 //check i
    3822                                 if (i<0 || i>=nbt ){
    3823                                         ISSMERROR("Index problem in NewPoints (i=%i not in [0 %i])",i,nbt-1);
    3824                                 }
    3825                                 //change first_np_or_next_t[i]
    3826                                 first_np_or_next_t[i] = iter;
    3827 
    3828                                 //Loop over the edges of t
    3829                                 for(j=0;j<3;j++){
    3830                                         TriangleAdjacent tj(t,j);
    3831                                         MeshVertex &vA = *tj.EdgeVertex(0);
    3832                                         MeshVertex &vB = *tj.EdgeVertex(1);
    3833 
    3834                                         //if t is a boundary triangle, or tj locked, continue
    3835                                         if (!t->link)     continue;
    3836                                         if (t->det <0)    continue;
    3837                                         if (t->Locked(j)) continue;
    3838 
    3839                                         TriangleAdjacent tadjj = t->Adj(j);       
    3840                                         Triangle* ta=tadjj;
    3841 
    3842                                         //if the adjacent triangle is a boundary triangle, continur
    3843                                         if (ta->det<0) continue;         
    3844 
    3845                                         R2 A=vA;
    3846                                         R2 B=vB;
    3847                                         k=Number(ta);
    3848 
    3849                                         //if this edge has already been done, go to next edge of triangle
    3850                                         if(first_np_or_next_t[k]==iter) continue;
    3851 
    3852                                         lIntTria.SplitEdge(Bh,A,B);
    3853                                         lIntTria.NewPoints(vertices,nbv,nbvx);
    3854                                   } // end loop for each edge
    3855                           }// for triangle   
    3856 
    3857                         if (!InsertNewPoints(nbvold,NbTSwap)) break;
    3858                         for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
    3859                         Headt = nbt; // empty list
    3860 
    3861                         // for all the triangle containing the vertex i
    3862                         for (i=nbvold;i<nbv;i++){
    3863                                 MeshVertex*          s  = vertices + i;
    3864                                 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
    3865                                 Triangle*        tbegin= (Triangle*) ta;
    3866                                 long kt;
    3867                                 do {
    3868                                         kt = Number((Triangle*) ta);
    3869                                         if (first_np_or_next_t[kt]>0){
    3870                                                 first_np_or_next_t[kt]=-Headt;
    3871                                                 Headt=kt;
    3872                                         }
    3873                                         if (ta.EdgeVertex(0)!=s){
    3874                                                 ISSMERROR("ta.EdgeVertex(0)!=s");
    3875                                         }
    3876                                         ta = Next(Adj(ta));
    3877                                 } while ( (tbegin != (Triangle*) ta));
    3878                         }
    3879 
    3880                 } while (nbv!=nbvold);
    3881 
    3882                 delete []  first_np_or_next_t;
    3883 
    3884                 long NbSwapf =0,NbSwp;
    3885 
    3886                 NbSwp = NbSwapf;
    3887                 for (i=0;i<nbv;i++)
    3888                  NbSwapf += vertices[i].Optim(0);
    3889                 NbTSwap +=  NbSwapf ;
    3890         }
    3891         /*}}}1*/
    3892 /*FUNCTION Mesh::PreInit{{{1*/
    3893 void Mesh::PreInit(long inbvx) {
    3894         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/
    3895 
    3896         long int verbose=0;
    3897 
    3898         srand(19999999);
    3899         NbRef=0;
    3900         //  allocGeometry=0;
    3901         NbOfTriangleSearchFind =0;
    3902         NbOfSwapTriangle =0;
    3903         nbiv=0;
    3904         nbv=0;
    3905         nbvx=inbvx;
    3906         nbt=0;
    3907         NbOfQuad = 0;
    3908         nbtx=2*inbvx-2;
    3909         NbSubDomains=0;
    3910         NbVertexOnBThVertex=0;
    3911         NbVertexOnBThEdge=0;
    3912         VertexOnBThVertex=NULL;
    3913         VertexOnBThEdge=NULL;
    3914 
    3915         NbCrackedVertices=0;
    3916         CrackedVertices  =NULL; 
    3917         NbCrackedEdges   =0;
    3918         CrackedEdges     =NULL; 
    3919         nbe = 0;
    3920 
    3921         if (inbvx) {
    3922                 vertices=new MeshVertex[nbvx];
    3923                 if (!vertices){
    3924                         ISSMERROR("!vertices");
    3925                 }
    3926                 ordre=new (MeshVertex* [nbvx]);
    3927                 if (!ordre){
    3928                         ISSMERROR("!ordre");
    3929                 }
    3930                 triangles=new Triangle[nbtx];
    3931                 if (!triangles){
    3932                         ISSMERROR("!triangles");
    3933                 }
    3934         }
    3935         else {
    3936                 vertices=NULL;
    3937                 ordre=NULL;
    3938                 triangles=NULL;
    3939                 nbtx=0;
    3940         }
    3941 
    3942         quadtree=NULL;
    3943         edges=NULL;
    3944         VerticesOnGeomVertex=NULL;
    3945         VerticesOnGeomEdge=NULL;
    3946         NbVerticesOnGeomVertex=0;
    3947         NbVerticesOnGeomEdge=0;
    3948         subdomains=NULL;
    3949         NbSubDomains=0;
    3950 }
    3951 /*}}}1*/
    3952         /*FUNCTION Mesh::ProjectOnCurve{{{1*/
    3953         GeometricalEdge*   Mesh::ProjectOnCurve( Edge & BhAB, MeshVertex &  vA, MeshVertex & vB,
    3954                                 double theta,MeshVertex & R,VertexOnEdge &  BR,VertexOnGeom & GR) {
    3955                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/ProjectOnCurve)*/
    3956 
    3957                 void *pA=0,*pB=0;
    3958                 double tA=0,tB=0;
    3959                 R2 A=vA,B=vB;
    3960                 MeshVertex * pvA=&vA, * pvB=&vB;
    3961                 if (vA.vint == IsVertexOnVertex){
    3962                         pA=vA.onBackgroundVertex;
    3963                 }
    3964                 else if (vA.vint == IsVertexOnEdge){
    3965                         pA=vA.onBackgroundEdge->be;
    3966                         tA=vA.onBackgroundEdge->abcisse;
    3967                 }
    3968                 else {
    3969                         ISSMERROR("ProjectOnCurve On MeshVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA));
    3970                 }
    3971 
    3972                 if (vB.vint == IsVertexOnVertex){
    3973                         pB=vB.onBackgroundVertex;
    3974                 }
    3975                 else if(vB.vint == IsVertexOnEdge){
    3976                         pB=vB.onBackgroundEdge->be;
    3977                         tB=vB.onBackgroundEdge->abcisse;
    3978                 }
    3979                 else {
    3980                         ISSMERROR("ProjectOnCurve On MeshVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB));
    3981                 }
    3982                 Edge * e = &BhAB;
    3983                 if (!pA || !pB || !e){
    3984                         ISSMERROR("!pA || !pB || !e");
    3985                 }
    3986                 // be carefull the back ground edge e is on same geom edge
    3987                 // of the initiale edge def by the 2 vertex A B;
    3988                 //check Is a background Mesh;   
    3989                 if (e<BTh.edges || e>=BTh.edges+BTh.nbe){
    3990                         ISSMERROR("e<BTh.edges || e>=BTh.edges+BTh.nbe");
    3991                 }
    3992                 // walk on BTh edge
    3993                 //not finish ProjectOnCurve with BackGround Mesh);
    3994                 // 1 first find a back ground edge contening the vertex A
    3995                 // 2 walk n back gound boundary to find the final vertex B
    3996 
    3997                 if( vA.vint == IsVertexOnEdge)
    3998                   { // find the start edge
    3999                         e = vA.onBackgroundEdge->be;     
    4000 
    4001                   }
    4002                 else if (vB.vint == IsVertexOnEdge)
    4003                   {
    4004                         theta = 1-theta;
    4005                         Exchange(tA,tB);
    4006                         Exchange(pA,pB);
    4007                         Exchange(pvA,pvB);
    4008                         Exchange(A,B);
    4009                         e =  vB.onBackgroundEdge->be;
    4010 
    4011                   }
    4012                 else{ // do the search by walking
    4013                         ISSMERROR("case not supported yet");
    4014                   }
    4015 
    4016                 // find the direction of walking with sens of edge and pA,PB;
    4017                 R2 AB=B-A;
    4018 
    4019                 double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);
    4020                 int kkk=0;
    4021                 int sens = (cosE01AB>0) ? 1 : 0;
    4022 
    4023                 //   double l=0; // length of the edge AB
    4024                 double abscisse = -1;
    4025 
    4026                 for (int cas=0;cas<2;cas++){
    4027                         // 2 times algo:
    4028                         //    1 for computing the length l
    4029                         //    2 for find the vertex
    4030                         int  iii;
    4031                         MeshVertex  *v0=pvA,*v1;
    4032                         Edge *neee,*eee;
    4033                         double lg =0; // length of the curve
    4034                         double te0;
    4035                         // we suppose take the curve's abcisse
    4036                         for ( eee=e,iii=sens,te0=tA;
    4037                                                 eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ;
    4038                                                 neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) {
    4039 
    4040                                 kkk=kkk+1;
    4041                                 if (kkk>=100){
    4042                                         ISSMERROR("kkk>=100");
    4043                                 }
    4044                                 if (!eee){
    4045                                         ISSMERROR("!eee");
    4046                                 }
    4047                                 double lg0 = lg;
    4048                                 double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
    4049                                 lg += dp;
    4050                                 if (cas && abscisse <= lg) { // ok we find the geom edge
    4051                                         double sss  =   (abscisse-lg0)/dp;
    4052                                         double thetab = te0*(1-sss)+ sss*iii;
    4053                                         if (thetab<0 || thetab>1){
    4054                                                 ISSMERROR("thetab<0 || thetab>1");
    4055                                         }
    4056                                         BR = VertexOnEdge(&R,eee,thetab);
    4057                                         return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
    4058                                   }
    4059                           }
    4060                         // we find the end
    4061                         if (v1 != pvB){
    4062                                 if (( void*) v1 == pB)
    4063                                  tB = iii;
    4064 
    4065                                 double lg0 = lg;
    4066                                 if (!eee){
    4067                                         ISSMERROR("!eee");
    4068                                 }
    4069                                 v1 = pvB;
    4070                                 double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
    4071                                 lg += dp;       
    4072                                 abscisse = lg*theta;
    4073                                 if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end
    4074                                   { // ok we find the geom edge
    4075                                         double sss  =   (abscisse-lg0)/dp;
    4076                                         double thetab = te0*(1-sss)+ sss*tB;
    4077                                         if (thetab<0 || thetab>1){
    4078                                                 ISSMERROR("thetab<0 || thetab>1");
    4079                                         }
    4080                                         BR = VertexOnEdge(&R,eee,thetab);
    4081                                         return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
    4082                                   }
    4083                           }
    4084                         abscisse = lg*theta;
    4085 
    4086                   }
    4087                 ISSMERROR("Big bug...");
    4088                 return 0; // just for the compiler
    4089         }                 
    4090         /*}}}1*/
    4091 /*FUNCTION Mesh::ReconstructExistingMesh{{{1*/
    4092 void Mesh::ReconstructExistingMesh(){
    4093         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/
    4094 
    4095         /*This routine reconstruct an existing mesh to make it CONVEX:
    4096          * -all the holes are filled
    4097          * -concave boundaries are filled
    4098          * A convex mesh is required for a lot of operations. This is why every mesh
    4099          * goes through this process.
    4100          * This routine also generates mesh properties such as adjencies,...
    4101          */
    4102 
    4103         /*Intermediary*/
    4104         int verbose=0;
    4105 
    4106         // generation of the integer coordinate
    4107 
    4108         // find extrema coordinates of vertices pmin,pmax
    4109         long i;
    4110         if(verbose>2) printf("      Reconstruct mesh of %i vertices\n",nbv);
    4111 
    4112         //initialize ordre
    4113         ISSMASSERT(ordre);
    4114         for (i=0;i<nbv;i++) ordre[i]=0;
    4115 
    4116         //Initialize NbSubDomains
    4117         NbSubDomains =0;
    4118 
    4119         /* generation of triangles adjacency*/
    4120 
    4121         //First add existing edges
    4122         long kk =0;
    4123         SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
    4124         for (i=0;i<nbe;i++){
    4125                 kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
    4126         }
    4127         if (kk != nbe){
    4128                 ISSMERROR("There are %i double edges in the mesh",kk-nbe);
    4129         }
    4130 
    4131         //Add edges of all triangles in existing mesh
    4132         long* st = new long[nbt*3];
    4133         for (i=0;i<nbt*3;i++) st[i]=-1;
    4134         for (i=0;i<nbt;i++){
    4135                 for (int j=0;j<3;j++){
    4136 
    4137                         //Add current triangle edge to edge4
    4138                         long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
    4139 
    4140                         long invisible=triangles[i].Hidden(j);
    4141 
    4142                         //If the edge has not been added to st, add it
    4143                         if(st[k]==-1) st[k]=3*i+j;
    4144 
    4145                         //If the edge already exists, add adjacency
    4146                         else if(st[k]>=0) {
    4147                                 ISSMASSERT(!triangles[i].TriangleAdj(j));
    4148                                 ISSMASSERT(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
    4149 
    4150                                 triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));
    4151                                 if (invisible) triangles[i].SetHidden(j);
    4152                                 if (k<nbe)     triangles[i].SetLocked(j);
    4153 
    4154                                 //Make st[k] negative so that it will throw an error message if it is found again
    4155                                 st[k]=-2-st[k];
    4156                         }
    4157 
    4158                         //An edge belongs to 2 triangles
    4159                         else {
    4160                                 ISSMERROR("The edge (%i , %i) belongs to more than 2 triangles",Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
    4161                         }
    4162                 }
    4163         }
    4164 
    4165         //Display info if required
    4166         if(verbose>5) {
    4167                 printf("         info of Mesh:\n");
    4168                 printf("            - number of vertices    = %i \n",nbv);
    4169                 printf("            - number of triangles   = %i \n",nbt);
    4170                 printf("            - number of given edges = %i \n",nbe);
    4171                 printf("            - number of all edges   = %i \n"  ,edge4->nb());
    4172                 printf("            - Euler number 1 - nb of holes = %i \n"  ,nbt-edge4->nb()+nbv);
    4173         }
    4174 
    4175         //check the consistency of edge[].adj and the geometrical required vertex
    4176         long k=0;
    4177         for (i=0;i<edge4->nb();i++){
    4178                 if (st[i]>=0){ // edge alone
    4179                         if (i<nbe){
    4180                                 long i0=edge4->i(i);
    4181                                 ordre[i0] = vertices+i0;
    4182                                 long i1=edge4->j(i);
    4183                                 ordre[i1] = vertices+i1;
    4184                         }
    4185                         else {
    4186                                 k=k+1;
    4187                                 if (k<10) {
    4188                                         //print only 10 edges
    4189                                         printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));
    4190                                 }
    4191                                 else if (k==10){
    4192                                         printf("Other lost boundary edges not shown...\n");
    4193                                 }
    4194                         }
    4195                 }
    4196         }
    4197         if(k) {
    4198                 ISSMERROR("%i boundary edges (from the geometry) are not defined as mesh edges",k);
    4199         }
    4200 
    4201         /* mesh generation with boundary points*/
    4202         long nbvb=0;
    4203         for (i=0;i<nbv;i++){
    4204                 vertices[i].t=0;
    4205                 vertices[i].vint=0;
    4206                 if (ordre[i]) ordre[nbvb++]=ordre[i];
    4207         }
    4208 
    4209         Triangle* savetriangles=triangles;
    4210         long savenbt=nbt;
    4211         long savenbtx=nbtx;
    4212         SubDomain* savesubdomains=subdomains;
    4213         subdomains=0;
    4214 
    4215         long  Nbtriafillhole=2*nbvb;
    4216         Triangle* triafillhole=new Triangle[Nbtriafillhole];
    4217         triangles = triafillhole;
    4218 
    4219         nbt=2;
    4220         nbtx= Nbtriafillhole;
    4221 
    4222         //Find a vertex that is not aligned with vertices 0 and 1
    4223         for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;)
    4224          if  (++i>=nbvb) {
    4225                  ISSMERROR("ReconstructExistingMesh: All the vertices are aligned");
    4226          }
    4227         //Move this vertex (i) to the 2d position in ordre
    4228         Exchange(ordre[2], ordre[i]);
    4229 
    4230         /*Reconstruct mesh beginning with 2 triangles*/
    4231         MeshVertex *  v0=ordre[0], *v1=ordre[1];
    4232 
    4233         triangles[0](0) = NULL; // Infinite vertex
    4234         triangles[0](1) = v0;
    4235         triangles[0](2) = v1;
    4236 
    4237         triangles[1](0) = NULL;// Infinite vertex
    4238         triangles[1](2) = v0;
    4239         triangles[1](1) = v1;
    4240         const int e0 = OppositeEdge[0];
    4241         const int e1 = NextEdge[e0];
    4242         const int e2 = PreviousEdge[e0];
    4243         triangles[0].SetAdj2(e0, &triangles[1] ,e0);
    4244         triangles[0].SetAdj2(e1, &triangles[1] ,e2);
    4245         triangles[0].SetAdj2(e2, &triangles[1] ,e1);
    4246 
    4247         triangles[0].det = -1;  // boundary triangles
    4248         triangles[1].det = -1;  // boundary triangles
    4249 
    4250         triangles[0].SetTriangleContainingTheVertex();
    4251         triangles[1].SetTriangleContainingTheVertex();
    4252 
    4253         triangles[0].link=&triangles[1];
    4254         triangles[1].link=&triangles[0];
    4255 
    4256         if (!quadtree) delete quadtree; //ReInitialise;
    4257         quadtree = new QuadTree(this,0);
    4258         quadtree->Add(*v0);
    4259         quadtree->Add(*v1);
    4260 
    4261         // vertices are added one by one
    4262         long NbSwap=0;
    4263         for (int icount=2; icount<nbvb; icount++) {
    4264                 MeshVertex *vi  = ordre[icount];
    4265                 Icoor2 dete[3];
    4266                 Triangle *tcvi = FindTriangleContaining(vi->i,dete);
    4267                 quadtree->Add(*vi);
    4268                 AddVertex(*vi,tcvi,dete);
    4269                 NbSwap += vi->Optim(1,1);
    4270         }
    4271 
    4272         //enforce the boundary
    4273         TriangleAdjacent ta(0,0);
    4274         long nbloss = 0,knbe=0;
    4275         for ( i = 0; i < nbe; i++){
    4276                 if (st[i] >=0){ //edge alone => on border
    4277                         MeshVertex &a=edges[i][0], &b=edges[i][1];
    4278                         if (a.t && b.t){
    4279                                 knbe++;
    4280                                 if (ForceEdge(a,b,ta)<0) nbloss++;
    4281                         }
    4282                 }
    4283         }
    4284         if(nbloss) {
    4285                 ISSMERROR("we lost %i existing edges other %i",nbloss,knbe);
    4286         }
    4287 
    4288         FindSubDomain(1);
    4289         // remove all the hole
    4290         // remove all the good sub domain
    4291         long krm =0;
    4292         for (i=0;i<nbt;i++){
    4293                 if (triangles[i].link){ // remove triangles
    4294                         krm++;
    4295                         for (int j=0;j<3;j++){
    4296                                 TriangleAdjacent ta =  triangles[i].Adj(j);
    4297                                 Triangle &tta = *(Triangle*)ta;
    4298                                 //if edge between remove and not remove
    4299                                 if(! tta.link){
    4300                                         // change the link of ta;
    4301                                         int ja = ta;
    4302                                         MeshVertex *v0= ta.EdgeVertex(0);
    4303                                         MeshVertex *v1= ta.EdgeVertex(1);
    4304                                         long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
    4305 
    4306                                         ISSMASSERT(st[k]>=0);
    4307                                         tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
    4308                                         ta.SetLock();
    4309                                         st[k]=-2-st[k];
    4310                                 }
    4311                         }
    4312                 }
    4313         }
    4314         long NbTfillHoll =0;
    4315         for (i=0;i<nbt;i++){
    4316                 if (triangles[i].link) {
    4317                         triangles[i]=Triangle((MeshVertex *) NULL,(MeshVertex *) NULL,(MeshVertex *) NULL);
    4318                         triangles[i].color=-1;
    4319                 }
    4320                 else{
    4321                         triangles[i].color= savenbt+ NbTfillHoll++;
    4322                 }
    4323         }
    4324         ISSMASSERT(savenbt+NbTfillHoll<=savenbtx);
    4325 
    4326         // copy of the outside triangles in saveMesh
    4327         for (i=0;i<nbt;i++){
    4328                 if(triangles[i].color>=0) {
    4329                         savetriangles[savenbt]=triangles[i];
    4330                         savetriangles[savenbt].link=0;
    4331                         savenbt++;
    4332                 }
    4333         }
    4334         // gestion of the adj
    4335         k =0;
    4336         Triangle * tmax = triangles + nbt;
    4337         for (i=0;i<savenbt;i++) {
    4338                 Triangle & ti = savetriangles[i];
    4339                 for (int j=0;j<3;j++){
    4340                         Triangle * ta = ti.TriangleAdj(j);
    4341                         int aa = ti.NuEdgeTriangleAdj(j);
    4342                         int lck = ti.Locked(j);
    4343                         if (!ta) k++; // bug
    4344                         else if ( ta >= triangles && ta < tmax){
    4345                                 ta= savetriangles + ta->color;
    4346                                 ti.SetAdj2(j,ta,aa);
    4347                                 if(lck) ti.SetLocked(j);
    4348                         }
    4349                 }
    4350         }
    4351 
    4352         // restore triangles;
    4353         nbt=savenbt;
    4354         nbtx=savenbtx;
    4355         delete [] triangles;
    4356         delete [] subdomains;
    4357         triangles = savetriangles;
    4358         subdomains = savesubdomains;
    4359         if (k) {
    4360                 ISSMERROR("number of triangles edges alone = %i",k);
    4361         }
    4362         FindSubDomain();
    4363 
    4364         delete edge4;
    4365         delete [] st;
    4366         for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);
    4367 
    4368         SetVertexFieldOn();
    4369 
    4370         /*Check requirements consistency*/
    4371         for (i=0;i<nbe;i++){
    4372         /*If the current mesh edge is on Geometry*/
    4373                 if(edges[i].onGeometry){
    4374                         for(int j=0;j<2;j++){
    4375                                 /*Go through the edges adjacent to current edge (if on the same curve)*/
    4376                                 if (!edges[i].adj[j]){
    4377                                         /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/
    4378                                         /*Check that the 2 vertices are on geometry AND required*/
    4379                                         if(!edges[i][j].onGeometry->IsRequiredVertex()){
    4380                                                 printf("ReconstructExistingMesh error message: problem with the edge number %i: [%i %i]\n",i+1,Number(edges[i][0])+1,Number(edges[i][1])+1);
    4381                                                 printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].onGeometry)+1);
    4382                                                 if (edges[i][j].onGeometry->OnGeomVertex())
    4383                                                  printf("the vertex number %i of this edge is a geometric MeshVertex number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->gv)+1);
    4384                                                 else if (edges[i][j].onGeometry->OnGeomEdge())
    4385                                                  printf("the vertex number %i of this edge is a geometric Edge number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->ge)+1);
    4386                                                 else
    4387                                                  printf("Its pointer is %p\n",edges[i][j].onGeometry);
    4388 
    4389                                                 printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");
    4390                                                 ISSMERROR("See above (might be cryptic...)");
    4391                                         }
    4392                                 }
    4393                         }
    4394                 }
    4395         }
    4396 }
    4397 /*}}}1*/
    4398         /*FUNCTION Mesh::ReNumberingTheTriangleBySubDomain{{{1*/
    4399         void Mesh::ReNumberingTheTriangleBySubDomain(bool justcompress){
    4400                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
    4401 
    4402                 long int verbose=0;
    4403                 long *renu= new long[nbt];
    4404                 register Triangle *t0,*t,*te=triangles+nbt;
    4405                 register long k=0,it,i,j;
    4406 
    4407                 for ( it=0;it<nbt;it++)
    4408                  renu[it]=-1; // outside triangle
    4409                 for ( i=0;i<NbSubDomains;i++)
    4410                   {
    4411                         t=t0=subdomains[i].head;
    4412                         if (!t0){ // not empty sub domain
    4413                                 ISSMERROR("!t0");
    4414                         }
    4415                         do {
    4416                                 long kt = Number(t);
    4417                                 if (kt<0 || kt >= nbt ){
    4418                                         ISSMERROR("kt<0 || kt >= nbt");
    4419                                 }
    4420                                 if (renu[kt]!=-1){
    4421                                         ISSMERROR("renu[kt]!=-1");
    4422                                 }
    4423                                 renu[kt]=k++;
    4424                         }
    4425                         while (t0 != (t=t->link));
    4426                   }
    4427                 // take is same numbering if possible   
    4428                 if(justcompress)
    4429                  for ( k=0,it=0;it<nbt;it++)
    4430                   if(renu[it] >=0 )
    4431                         renu[it]=k++;
    4432 
    4433                 // put the outside triangles at the end
    4434                 for ( it=0;it<nbt;it++){
    4435                         if (renu[it]==-1) renu[it]=k++;
    4436                 }
    4437                 if (k != nbt){
    4438                         ISSMERROR("k != nbt");
    4439                 }
    4440                 // do the change on all the pointeur
    4441                 for ( it=0;it<nbt;it++)
    4442                  triangles[it].ReNumbering(triangles,te,renu);
    4443 
    4444                 for ( i=0;i<NbSubDomains;i++)
    4445                  subdomains[i].head=triangles+renu[Number(subdomains[i].head)];
    4446 
    4447                 // move the Triangles  without a copy of the array
    4448                 // be carefull not trivial code
    4449                 for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu
    4450                  if (renu[it] >= 0) // a new sub cycle
    4451                         {
    4452                          i=it;
    4453                          Triangle ti=triangles[i],tj;
    4454                          while ( (j=renu[i]) >= 0)
    4455                                 { // i is old, and j is new
    4456                                  renu[i] = -1; // mark
    4457                                  tj = triangles[j]; // save new
    4458                                  triangles[j]= ti; // new <- old
    4459                                  i=j;     // next
    4460                                  ti = tj;
    4461                                 } 
    4462                         }
    4463                 delete [] renu;
    4464                 nt = nbt - NbOutT;
    4465 
    4466         }
    4467         /*}}}1*/
    4468         /*FUNCTION Mesh::ReNumberingVertex{{{1*/
    4469         void Mesh::ReNumberingVertex(long * renu) {
    4470                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
    4471 
    4472                 // warning be carfull because pointer
    4473                 // from on mesh to over mesh
    4474                 //  --  so do ReNumbering at the beginning
    4475                 MeshVertex * ve = vertices+nbv;
    4476                 long it,ie,i;
    4477 
    4478                 printf("renumbering triangles\n");
    4479                 for ( it=0;it<nbt;it++)
    4480                  triangles[it].ReNumbering(vertices,ve,renu);
    4481 
    4482                 printf("renumbering edges\n");
    4483                 for ( ie=0;ie<nbe;ie++)
    4484                  edges[ie].ReNumbering(vertices,ve,renu);
    4485 
    4486                 printf("renumbering vertices on geom\n");
    4487                 for (i=0;i< NbVerticesOnGeomVertex;i++)
    4488                   {
    4489                         MeshVertex *v = VerticesOnGeomVertex[i].mv;
    4490                         if (v>=vertices && v < ve)
    4491                          VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)];
    4492                   }
    4493 
    4494                 printf("renumbering vertices on edge\n");
    4495                 for (i=0;i< NbVerticesOnGeomEdge;i++)
    4496                   {
    4497                         MeshVertex *v =VerticesOnGeomEdge[i].mv;
    4498                         if (v>=vertices && v < ve)
    4499                          VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)];
    4500                   }
    4501 
    4502                 printf("renumbering vertices on Bth vertex\n");
    4503                 for (i=0;i< NbVertexOnBThVertex;i++)
    4504                   {
    4505                         MeshVertex *v=VertexOnBThVertex[i].v;
    4506                         if (v>=vertices && v < ve)
    4507                          VertexOnBThVertex[i].v=vertices+renu[Number(v)];
    4508                   }
    4509 
    4510                 for (i=0;i< NbVertexOnBThEdge;i++)
    4511                   {
    4512                         MeshVertex *v=VertexOnBThEdge[i].v;
    4513                         if (v>=vertices && v < ve)
    4514                          VertexOnBThEdge[i].v=vertices+renu[Number(v)];
    4515                   }
    4516 
    4517                 // move the Vertices without a copy of the array
    4518                 // be carefull not trivial code
    4519                 long j;
    4520                 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
    4521                  if (renu[it] >= 0) // a new sub cycle
    4522                         {
    4523                          i=it;
    4524                          MeshVertex ti=vertices[i],tj;
    4525                          while ( (j=renu[i]) >= 0){
    4526                                  // i is old, and j is new
    4527                                  renu[i] = -1-renu[i]; // mark
    4528                                  tj = vertices[j];     // save new
    4529                                  vertices[j]= ti;      // new <- old
    4530                                  i=j;     // next
    4531                                  ti = tj;
    4532                                 } 
    4533                         }
    4534                 if (quadtree){
    4535                         delete quadtree;
    4536                         quadtree = new QuadTree(this);
    4537                 }
    4538 
    4539                 for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
    4540         }
    4541         /*}}}1*/
    4542 /*FUNCTION Mesh::SetIntCoor{{{1*/
    4543 void Mesh::SetIntCoor(const char * strfrom) {
    4544         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)*/
    4545 
    4546         /*Set integer coordinate for existing vertices*/
    4547 
    4548         //Get extrema coordinates of the existing vertices
    4549         pmin =  vertices[0].r;
    4550         pmax =  vertices[0].r;
    4551         long i;
    4552         for (i=0;i<nbv;i++) {
    4553                 pmin.x = Min(pmin.x,vertices[i].r.x);
    4554                 pmin.y = Min(pmin.y,vertices[i].r.y);
    4555                 pmax.x = Max(pmax.x,vertices[i].r.x);
    4556                 pmax.y = Max(pmax.y,vertices[i].r.y);
    4557         }
    4558         R2 DD = (pmax-pmin)*0.05;
    4559         pmin = pmin-DD;
    4560         pmax = pmax+DD;
    4561 
    4562         //Compute coefIcoor
    4563         coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
    4564         if (coefIcoor<=0){
    4565                 ISSMERROR("coefIcoor should be positive, a problem in the geometry is likely");
    4566         }
    4567 
    4568         // generation of integer coord 
    4569         for (i=0;i<nbv;i++) {
    4570                 vertices[i].i = toI2(vertices[i].r);   
    4571         }
    4572 
    4573         // computation of the det
    4574         int number_of_errors=0;
    4575         for (i=0;i<nbt;i++) {
    4576                 MeshVertex & v0 = triangles[i][0];
    4577                 MeshVertex & v1 = triangles[i][1];
    4578                 MeshVertex & v2 = triangles[i][2];
    4579 
    4580                 //If this is not a boundary triangle
    4581                 if ( &v0 && &v1 &&  &v2 ){
    4582 
    4583                         /*Compute determinant*/
    4584                         triangles[i].det= det(v0,v1,v2);
    4585 
    4586                         /*Check that determinant is positive*/
    4587                         if (triangles[i].det <=0){
    4588 
    4589                                 /*increase number_of_errors and print error only for the first 20 triangles*/
    4590                                 number_of_errors++;
    4591                                 if (number_of_errors<20){
    4592                                         printf("Area of Triangle %i < 0 (det=%i)\n",i+1,triangles[i].det);
    4593                                 }
    4594                         }
    4595                 }
    4596 
    4597                 //else, set as -1
    4598                 else triangles[i].det=-1;
    4599         }
    4600 
    4601         if (number_of_errors) ISSMERROR("Fatal error: some triangles have negative areas, see above");
    4602 }
    4603 /*}}}1*/
    4604 /*FUNCTION Mesh::ShowRegulaty{{{1*/
    4605 void  Mesh::ShowRegulaty() const {
    4606         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/
    4607 
    4608         const  double  sqrt32=sqrt(3.)*0.5;
    4609         const double  aireKh=sqrt32*0.5;
    4610         D2  Beq(1,0),Heq(0.5,sqrt32);
    4611         D2xD2 Br(D2xD2(Beq,Heq).t());
    4612         D2xD2 B1r(Br.inv());
    4613         double gammamn=1e100,hmin=1e100;
    4614         double gammamx=0,hmax=0;
    4615         double beta=1e100;
    4616         double beta0=0;
    4617         double  alpha2=0;
    4618         double area=0,Marea=0;
    4619         // double cf= double(coefIcoor);
    4620         // double cf2= 6.*cf*cf;
    4621         int nt=0;
    4622         for (int it=0;it<nbt;it++)
    4623          if ( triangles[it].link)
    4624                 {
    4625                  nt++;
    4626                  Triangle &K=triangles[it];
    4627                  double  area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;
    4628                  area+= area3;
    4629                  D2xD2 B_Kt(K[0],K[1],K[2]);
    4630                  D2xD2 B_K(B_Kt.t());
    4631                  D2xD2 B1K = Br*B_K.inv();
    4632                  D2xD2 BK =  B_K*B1r;
    4633                  D2xD2 B1B1 = B1K.t()*B1K;
    4634                  Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
    4635                  MatVVP2x2 VMK(MK);
    4636                  alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
    4637                  double betaK=0;
    4638 
    4639                  for (int j=0;j<3;j++)
    4640                         {
    4641                          double he= Norme2(R2(K[j],K[(j+1)%3]));
    4642                          hmin=Min(hmin,he);
    4643                          hmax=Max(hmax,he);
    4644                          MeshVertex & v=K[j];
    4645                          D2xD2 M((Metric)v);
    4646                          betaK += sqrt(M.det());
    4647 
    4648                          D2xD2 BMB = BK.t()*M*BK;
    4649                          Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);
    4650                          MatVVP2x2 VM1(M1);
    4651                          gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
    4652                          gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);         
    4653                         }
    4654                  betaK *= area3;//  1/2 (somme sqrt(det))* area(K)
    4655                  Marea+= betaK;
    4656                  beta=min(beta,betaK);
    4657                  beta0=max(beta0,betaK);
    4658                 }   
    4659         area*=3;
    4660         gammamn=sqrt(gammamn);
    4661         gammamx=sqrt(gammamx);   
    4662         printf("   Adaptmesh info:\n");
    4663         printf("      number of triangles = %i\n",nt);
    4664         printf("      hmin = %g, hmax=%g\n",hmin,hmax);
    4665         printf("      area = %g, M area = %g, M area/( |Khat| nt) = %g\n",area,Marea, Marea/(aireKh*nt));
    4666         printf("      infinite-regularity(?): min = %g, max = %g\n",gammamn,gammamx);
    4667         printf("      anisomax = %g, beta max = %g, min = %g\n",pow(alpha2,0.5),1./pow(beta/aireKh,0.5), 1./pow(beta0/aireKh,0.5));
    4668 }
    4669 /*}}}1*/
    4670 /*FUNCTION Mesh::ShowHistogram{{{1*/
    4671 void  Mesh::ShowHistogram() const {
    4672         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/
    4673 
    4674         const long kmax=10;
    4675         const double llmin = 0.5,llmax=2;
    4676         const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);
    4677         long histo[kmax+1];
    4678         long i,it,k, nbedges =0;
    4679         for (i=0;i<=kmax;i++) histo[i]=0;
    4680         for (it=0;it<nbt;it++)
    4681          if ( triangles[it].link)
    4682                 {
    4683 
    4684                  for (int j=0;j<3;j++)
    4685                         {
    4686                          Triangle *ta = triangles[it].TriangleAdj(j);
    4687                          if ( !ta || !ta->link || Number(ta) >= it)
    4688                                 {
    4689                                  MeshVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
    4690                                  MeshVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
    4691                                  if ( !& vP || !&vQ) continue;
    4692                                  R2 PQ = vQ.r - vP.r;
    4693                                  double l = log(LengthInterpole(vP,vQ,PQ));
    4694                                  nbedges++;
    4695                                  k = (int) ((l - lmin)*delta);
    4696                                  k = Min(Max(k,0L),kmax);
    4697                                  histo[k]++;
    4698                                 }
    4699                         }
    4700                 } 
    4701         printf(" --- Histogram of the unit mesh,  nb of edges = %i\n",nbedges);
    4702         printf("      length of edge in   | %% of edge  | Nb of edges \n");
    4703         printf("      --------------------+-------------+-------------\n");
    4704         for   (i=0;i<=kmax;i++){
    4705                 if (i==0) printf("      %10i",0);
    4706                 else      printf("      %10g",exp(lmin+i/delta));
    4707                 if (i==kmax) printf("          +inf   ");
    4708                 else printf("      %10g",exp(lmin+(i+1)/delta));
    4709                 printf("|  %10g |\n",((long)  ((10000.0 * histo[i])/ nbedges))/100.0);
    4710                 printf("  %i\n",histo[i]);
    4711         }
    4712         printf("      --------------------+-------------+-------------\n");
    4713 }
    4714 /*}}}1*/
    4715 /*FUNCTION Mesh::SmoothingVertex{{{1*/
    4716 void Mesh::SmoothingVertex(int nbiter,double omega ) {
    4717         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/
    4718 
    4719         long int verbose=0;
    4720         //  if quatree exist remove it end reconstruct
    4721         if (quadtree) delete quadtree;
    4722         quadtree=0;
    4723         ReMakeTriangleContainingTheVertex();
    4724         Triangle vide; // a triangle to mark the boundary vertex
    4725         Triangle   ** tstart= new Triangle* [nbv];
    4726         long i,j,k;
    4727         //   attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide
    4728         if ( this == & BTh)
    4729          for ( i=0;i<nbv;i++)
    4730           tstart[i]=vertices[i].t;     
    4731         else
    4732          for ( i=0;i<nbv;i++)
    4733           tstart[i]=0;
    4734         for ( j=0;j<NbVerticesOnGeomVertex;j++ )
    4735          tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide;
    4736         for ( k=0;k<NbVerticesOnGeomEdge;k++ )
    4737          tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;
    4738         if(verbose>2) printf("   SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega);
    4739         for (k=0;k<nbiter;k++)
    4740           {
    4741                 long i,NbSwap =0;
    4742                 double delta =0;
    4743                 for ( i=0;i<nbv;i++)
    4744                  if (tstart[i] != &vide) // not a boundary vertex
    4745                   delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
    4746                 if (!NbOfQuad)
    4747                  for ( i=0;i<nbv;i++)
    4748                   if (tstart[i] != &vide) // not a boundary vertex
    4749                         NbSwap += vertices[i].Optim(1);
    4750                 if (verbose>3) printf("      move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap);
    4751           }
    4752 
    4753         delete [] tstart;
    4754         if (quadtree) quadtree= new QuadTree(this);
    4755 }
    4756 /*}}}1*/
    4757 /*FUNCTION Mesh::SmoothMetric{{{1*/
    4758 void Mesh::SmoothMetric(double raisonmax) {
    4759         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/
    4760 
    4761         long int verbose=0;
    4762 
    4763         if(raisonmax<1.1) return;
    4764         if(verbose > 1) printf("   Mesh::SmoothMetric raisonmax = %g\n",raisonmax);
    4765         ReMakeTriangleContainingTheVertex();
    4766         long i,j,kch,kk,ip;
    4767         long *first_np_or_next_t0 = new long[nbv];
    4768         long *first_np_or_next_t1 = new long[nbv];
    4769         long Head0 =0,Head1=-1;
    4770         double logseuil= log(raisonmax);
    4771 
    4772         for(i=0;i<nbv-1;i++)
    4773          first_np_or_next_t0[i]=i+1;
    4774         first_np_or_next_t0[nbv-1]=-1;// end;
    4775         for(i=0;i<nbv;i++)
    4776          first_np_or_next_t1[i]=-1;
    4777         kk=0;
    4778         while (Head0>=0&& kk++<100) {
    4779                 kch=0;
    4780                 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
    4781                         //  pour tous les triangles autour du sommet s
    4782                         register Triangle* t= vertices[i].t;
    4783                         if (!t){
    4784                                 ISSMERROR("!t");
    4785                         }
    4786                         MeshVertex & vi = vertices[i];
    4787                         TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);
    4788                         MeshVertex *pvj0 = ta.EdgeVertex(0);
    4789                         while (1) {
    4790                                 ta=Previous(Adj(ta));
    4791                                 if (vertices+i != ta.EdgeVertex(1)){
    4792                                         ISSMERROR("vertices+i != ta.EdgeVertex(1)");
    4793                                 }
    4794                                 MeshVertex & vj = *(ta.EdgeVertex(0));
    4795                                 if ( &vj ) {
    4796                                         j= &vj-vertices;
    4797                                         if (j<0 || j >= nbv){
    4798                                                 ISSMERROR("j<0 || j >= nbv");
    4799                                         }
    4800                                         R2 Aij = (R2) vj - (R2) vi;
    4801                                         double ll =  Norme2(Aij);
    4802                                         if (0) { 
    4803                                                 double hi = ll/vi.m(Aij);
    4804                                                 double hj = ll/vj.m(Aij);
    4805                                                 if(hi < hj)
    4806                                                   {
    4807                                                         double dh=(hj-hi)/ll;
    4808                                                         if (dh>logseuil) {
    4809                                                                 vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
    4810                                                                 if(first_np_or_next_t1[j]<0)
    4811                                                                  kch++,first_np_or_next_t1[j]=Head1,Head1=j;
    4812                                                         }
    4813                                                   }
    4814                                         }
    4815                                         else
    4816                                           {
    4817                                                 double li = vi.m(Aij);
    4818                                                 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
    4819                                                  if(first_np_or_next_t1[j]<0) // if the metrix change
    4820                                                   kch++,first_np_or_next_t1[j]=Head1,Head1=j;
    4821                                           }
    4822                                 }
    4823                                 if  ( &vj ==  pvj0 ) break;
    4824                         }
    4825                 }
    4826                 Head0 = Head1;
    4827                 Head1 = -1;
    4828                 Exchange(first_np_or_next_t0,first_np_or_next_t1);
    4829         }
    4830         if(verbose>2) printf("      number of iterations = %i\n",kch);
    4831         delete [] first_np_or_next_t0;
    4832         delete [] first_np_or_next_t1;
    4833 }
    4834 /*}}}1*/
    4835         /*FUNCTION Mesh::SplitElement{{{1*/
    4836         int  Mesh::SplitElement(int choice){
    4837                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/
    4838 
    4839                 long int verbose=0;
    4840 
    4841                 Direction NoDirOfSearch;
    4842                 const  int withBackground = &BTh != this && &BTh;
    4843 
    4844                 ReNumberingTheTriangleBySubDomain();
    4845                 //int nswap =0;
    4846                 const long nfortria( choice ? 4 : 6);
    4847                 if(withBackground)
    4848                   {
    4849                         BTh.SetVertexFieldOn();
    4850                         SetVertexFieldOnBTh();
    4851                   }
    4852                 else
    4853                  BTh.SetVertexFieldOn();
    4854 
    4855                 long newnbt=0,newnbv=0;
    4856                 long * kedge = 0;
    4857                 long newNbOfQuad=NbOfQuad;
    4858                 long * ksplit= 0, * ksplitarray=0;
    4859                 long kkk=0;
    4860                 int ret =0;
    4861                 if (nbvx<nbv+nbe) return 1;//   
    4862                 // 1) create  the new points by spliting the internal edges
    4863                 // set the
    4864                 long nbvold = nbv;
    4865                 long nbtold = nbt;
    4866                 long NbOutTold  = NbOutT;
    4867                 long  NbEdgeOnGeom=0;
    4868                 long i;
    4869 
    4870                 nbt = nbt - NbOutT; // remove all the  the ouside triangles
    4871                 long nbtsave = nbt;
    4872                 Triangle * lastT = triangles + nbt;
    4873                 for (i=0;i<nbe;i++)
    4874                  if(edges[i].onGeometry) NbEdgeOnGeom++;
    4875                 long newnbe=nbe+nbe;
    4876                 //  long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
    4877                 long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;
    4878                 // long newNbVertexOnBThVertex=NbVertexOnBThVertex;
    4879                 long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;
    4880 
    4881                 // do allocation for pointeur to the geometry and background
    4882                 VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];
    4883                 VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ?  new VertexOnEdge[newNbVertexOnBThEdge]:0;
    4884                 if (NbVerticesOnGeomEdge)
    4885                  memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);
    4886                 if (NbVertexOnBThEdge)
    4887                  memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);
    4888                 Edge *newedges = new Edge [newnbe];
    4889                 //  memcpy(newedges,edges,sizeof(Edge)*nbe);
    4890                 SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);
    4891                 long k=nbv;
    4892                 long kk=0;
    4893                 long kvb = NbVertexOnBThEdge;
    4894                 long kvg = NbVerticesOnGeomEdge;
    4895                 long ie =0;
    4896                 Edge ** edgesGtoB=0;
    4897                 if (withBackground)
    4898                  edgesGtoB= BTh.MakeGeometricalEdgeToEdge();
    4899                 long ferr=0;
    4900                 for (i=0;i<nbe;i++)
    4901                  newedges[ie].onGeometry=0;
    4902 
    4903                 for (i=0;i<nbe;i++)
    4904                   {
    4905                         GeometricalEdge *ong =  edges[i].onGeometry;
    4906 
    4907                         newedges[ie]=edges[i];
    4908                         newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;
    4909                         newedges[ie].adj[1]=newedges + ie +1;
    4910                         R2 A = edges[i][0],B = edges[i][1];
    4911 
    4912 
    4913                         kk += (i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
    4914                         if (ong) // a geometrical edges
    4915                           {
    4916                                 if (withBackground){
    4917                                         // walk on back ground mesh
    4918                                         //  newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge);
    4919                                         // a faire -- difficile
    4920                                         // the first PB is to now a background edge between the 2 vertices
    4921                                         if (!edgesGtoB){
    4922                                                 ISSMERROR("!edgesGtoB");
    4923                                         }
    4924                                         ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].onGeometry)],
    4925                                                                 edges[i][0],edges[i][1],0.5,vertices[k],
    4926                                                                 newVertexOnBThEdge[kvb],
    4927                                                                 newVerticesOnGeomEdge[kvg++]);
    4928                                         vertices[k].ReferenceNumber= edges[i].ref;
    4929                                         vertices[k].DirOfSearch =   NoDirOfSearch;       
    4930                                         ;
    4931                                         // get the Info on background mesh
    4932                                         double s =        newVertexOnBThEdge[kvb];
    4933                                         MeshVertex &  bv0  = newVertexOnBThEdge[kvb][0];
    4934                                         MeshVertex &  bv1  = newVertexOnBThEdge[kvb][1];
    4935                                         // compute the metrix of the new points
    4936                                         vertices[k].m =  Metric(1-s,bv0,s,bv1);
    4937                                         kvb++;
    4938                                   }
    4939                                 else
    4940                                   {
    4941                                         ong=Gh.ProjectOnCurve(edges[i],
    4942                                                                 0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);
    4943                                         // vertices[k].i = toI2( vertices[k].r);
    4944                                         vertices[k].ReferenceNumber = edges[i].ref;
    4945                                         vertices[k].DirOfSearch = NoDirOfSearch;
    4946                                         vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);             
    4947                                   } 
    4948                           }
    4949                         else // straigth line edge ---
    4950                           {
    4951                                 vertices[k].r = ((R2) edges[i][0] + (R2)  edges[i][1] )*0.5;
    4952                                 vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);
    4953                                 vertices[k].onGeometry = 0;
    4954                           }
    4955                         //vertices[k].i = toI2( vertices[k].r);
    4956                         R2 AB =  vertices[k].r;
    4957                         R2 AA = (A+AB)*0.5;
    4958                         R2 BB = (AB+B)*0.5;
    4959                         vertices[k].ReferenceNumber = edges[i].ref;
    4960                         vertices[k].DirOfSearch = NoDirOfSearch;
    4961 
    4962                         newedges[ie].onGeometry = Gh.Containing(AA,ong);
    4963                         newedges[ie++].v[1]=vertices+k;
    4964 
    4965                         newedges[ie]=edges[i];
    4966                         newedges[ie].adj[0]=newedges + ie -1;
    4967                         newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;
    4968                         newedges[ie].onGeometry =  Gh.Containing(BB,ong);
    4969                         newedges[ie++].v[0]=vertices+k;
    4970                         k++;
    4971                   }
    4972                 if (edgesGtoB) delete [] edgesGtoB;
    4973                 edgesGtoB=0;
    4974 
    4975                 newnbv=k;
    4976                 newNbVerticesOnGeomEdge=kvg;
    4977                 if (newnbv> nbvx) goto Error;// bug
    4978 
    4979                 nbv = k;
    4980 
    4981 
    4982                 kedge = new long[3*nbt+1];
    4983                 ksplitarray = new long[nbt+1];
    4984                 ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0]
    4985 
    4986                 for (i=0;i<3*nbt;i++)
    4987                  kedge[i]=-1;
    4988 
    4989                 // 
    4990 
    4991                 for (i=0;i<nbt;i++) {
    4992                         Triangle & t = triangles[i];
    4993                         if (!t.link){
    4994                                 ISSMERROR("!t.link");
    4995                         }
    4996                         for(int j=0;j<3;j++)
    4997                           {
    4998                                 const TriangleAdjacent ta = t.Adj(j);
    4999                                 const Triangle & tt = ta;
    5000                                 if (&tt >= lastT)
    5001                                  t.SetAdj2(j,0,0);// unset adj
    5002                                 const MeshVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
    5003                                 const MeshVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
    5004                                 long  ke =edge4->SortAndFind(Number(v0),Number(v1));
    5005                                 if (ke>0)
    5006                                   {
    5007                                         long ii = Number(tt);
    5008                                         int  jj = ta;
    5009                                         long ks = ke + nbvold;
    5010                                         kedge[3*i+j] = ks;
    5011                                         if (ii<nbt) // good triangle
    5012                                          kedge[3*ii+jj] = ks;
    5013                                         MeshVertex &A=vertices[ks];
    5014                                         double aa,bb,cc,dd;
    5015                                         if ((dd=Area2(v0.r,v1.r,A.r)) >=0){
    5016                                                 // warning PB roundoff error
    5017                                                 if (t.link && ( (aa=Area2( A.r    , t[1].r , t[2].r )) < 0.0
    5018                                                                                 ||   (bb=Area2( t[0].r , A.r    , t[2].r )) < 0.0 
    5019                                                                                 ||   (cc=Area2( t[0].r , t[1].r , A.r    )) < 0.0)){
    5020                                                         printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd);
    5021                                                         ISSMERROR("Number of triangles with P2 interpolation Problem");
    5022                                                 }
    5023                                         }
    5024                                         else {
    5025                                                 if (tt.link && ( (aa=Area2( A.r     , tt[1].r , tt[2].r )) < 0
    5026                                                                                 ||   (bb=Area2( tt[0].r , A.r     , tt[2].r )) < 0
    5027                                                                                 ||   (cc=Area2( tt[0].r , tt[1].r , A.r     )) < 0)){
    5028                                                         printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd);
    5029                                                         ISSMERROR("Number of triangles with P2 interpolation Problem");
    5030                                                 }
    5031                                         }
    5032                                   }
    5033                           }
    5034                 }
    5035 
    5036                 for (i=0;i<nbt;i++){
    5037                         ksplit[i]=1; // no split by default
    5038                         const Triangle & t = triangles[ i];
    5039                         int nbsplitedge =0;
    5040                         int nbinvisible =0;
    5041                         int invisibleedge=0;
    5042                         int kkk[3];     
    5043                         for (int j=0;j<3;j++)
    5044                           {
    5045                                 if (t.Hidden(j)) invisibleedge=j,nbinvisible++;
    5046 
    5047                                 const TriangleAdjacent ta = t.Adj(j);
    5048                                 const Triangle & tt = ta;
    5049 
    5050 
    5051                                 const MeshVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
    5052                                 const MeshVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
    5053                                 if ( kedge[3*i+j] < 0)
    5054                                   {
    5055                                         long  ke =edge4->SortAndFind(Number(v0),Number(v1));
    5056                                         if (ke<0) // new
    5057                                           {
    5058                                                 if (&tt) // internal triangles all the boundary
    5059                                                   { // new internal edges
    5060                                                         long ii = Number(tt);
    5061                                                         int  jj = ta;
    5062 
    5063                                                         kedge[3*i+j]=k;// save the vertex number
    5064                                                         kedge[3*ii+jj]=k;
    5065                                                         if (k<nbvx)
    5066                                                           {
    5067                                                                 vertices[k].r = ((R2) v0+(R2) v1 )/2;
    5068                                                                 //vertices[k].i = toI2( vertices[k].r);
    5069                                                                 vertices[k].ReferenceNumber=0;
    5070                                                                 vertices[k].DirOfSearch =NoDirOfSearch;
    5071                                                                 vertices[k].m =  Metric(0.5,v0,0.5,v1);
    5072                                                           }
    5073                                                         k++;
    5074                                                         kkk[nbsplitedge++]=j;                 
    5075                                                   } // tt
    5076                                                 else
    5077                                                  ISSMERROR("Bug...");
    5078                                           } // ke<0           
    5079                                         else
    5080                                           { // ke >=0
    5081                                                 kedge[3*i+j]=nbvold+ke;
    5082                                                 kkk[nbsplitedge++]=j;// previously splited
    5083                                           }
    5084                                   }
    5085                                 else
    5086                                  kkk[nbsplitedge++]=j;// previously splited
    5087 
    5088                           }
    5089                         if (nbinvisible>=2){
    5090                                 ISSMERROR("nbinvisible>=2");
    5091                         }
    5092                         switch (nbsplitedge) {
    5093                                 case 0: ksplit[i]=10; newnbt++; break;   // nosplit
    5094                                 case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2
    5095                                 case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3
    5096                                 case 3:
    5097                                                   if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4;
    5098                                                   else   ksplit[i]=10*nfortria,newnbt+=nfortria;
    5099                                                   break;
    5100                         }
    5101                         if (ksplit[i]<40){
    5102                                 ISSMERROR("ksplit[i]<40");
    5103                         }
    5104                   }
    5105                 //  now do the element split
    5106                 newNbOfQuad = 4*NbOfQuad;
    5107                 nbv = k;
    5108                 kkk = nbt;
    5109                 ksplit[-1] = nbt;
    5110                 // look on  old true  triangles
    5111 
    5112                 for (i=0;i<nbtsave;i++){
    5113                         int  nbmkadj=0;
    5114                         long mkadj [100];
    5115                         mkadj[0]=i;
    5116                         long kk=ksplit[i]/10;
    5117                         int  ke=(int) (ksplit[i]%10);
    5118                         if (kk>=7 || kk<=0){
    5119                                 ISSMERROR("kk>=7 || kk<=0");
    5120                         }
    5121 
    5122                         // def the numbering   k (edge) i vertex
    5123                         int k0 = ke;
    5124                         int k1 = NextEdge[k0];
    5125                         int k2 = PreviousEdge[k0];
    5126                         int i0 = OppositeVertex[k0];
    5127                         int i1 = OppositeVertex[k1];
    5128                         int i2 = OppositeVertex[k2];
    5129 
    5130                         Triangle &t0=triangles[i];
    5131                         MeshVertex * v0=t0(i0);           
    5132                         MeshVertex * v1=t0(i1);           
    5133                         MeshVertex * v2=t0(i2);
    5134 
    5135                         if (nbmkadj>=10){
    5136                                 ISSMERROR("nbmkadj>=10");
    5137                         }
    5138                         // --------------------------
    5139                         TriangleAdjacent ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2));
    5140                         // save the flag Hidden
    5141                         int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)};
    5142                         // un set all adj -- save Hidden flag --
    5143                         t0.SetAdj2(0,0,hid[0]);
    5144                         t0.SetAdj2(1,0,hid[1]);
    5145                         t0.SetAdj2(2,0,hid[2]);
    5146                         // --  remake
    5147                         switch  (kk) {
    5148                                 case 1: break;// nothing
    5149                                 case 2: //
    5150                                                   {
    5151                                                         Triangle &t1=triangles[kkk++];
    5152                                                         t1=t0;
    5153                                                         if (kedge[3*i+i0]<0){
    5154                                                                 ISSMERROR("kedge[3*i+i0]<0");
    5155                                                         }
    5156                                                         MeshVertex * v3 = vertices + kedge[3*i+k0];
    5157 
    5158                                                         t0(i2) = v3;
    5159                                                         t1(i1) = v3;
    5160                                                         t0.SetAllFlag(k2,0);
    5161                                                         t1.SetAllFlag(k1,0);
    5162                                                   }
    5163                                                 break;
    5164                                 case 3: //
    5165                                                   {
    5166                                                         Triangle &t1=triangles[kkk++];
    5167                                                         Triangle &t2=triangles[kkk++];
    5168                                                         t2=t1=t0;
    5169                                                         if (kedge[3*i+k1]<0){
    5170                                                                 ISSMERROR("kedge[3*i+k1]<0");
    5171                                                         }
    5172                                                         if (kedge[3*i+k2]<0){
    5173                                                                 ISSMERROR("kedge[3*i+k2]<0");
    5174                                                         }
    5175 
    5176                                                         MeshVertex * v01 = vertices + kedge[3*i+k2];
    5177                                                         MeshVertex * v02 = vertices + kedge[3*i+k1];
    5178                                                         t0(i1) = v01;
    5179                                                         t0(i2) = v02;
    5180                                                         t1(i2) = v02;
    5181                                                         t1(i0) = v01;
    5182                                                         t2(i0) = v02;
    5183                                                         t0.SetAllFlag(k0,0);
    5184                                                         t1.SetAllFlag(k1,0);
    5185                                                         t1.SetAllFlag(k0,0);
    5186                                                         t2.SetAllFlag(k2,0);
    5187                                                   }
    5188                                                 break;
    5189                                 case 4: //
    5190                                 case 6: // split in 4
    5191                                                   {
    5192                                                         Triangle &t1=triangles[kkk++];
    5193                                                         Triangle &t2=triangles[kkk++];
    5194                                                         Triangle &t3=triangles[kkk++];
    5195                                                         t3=t2=t1=t0;
    5196                                                         if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){
    5197                                                                 ISSMERROR("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0");
    5198                                                         }
    5199                                                         MeshVertex * v12 = vertices + kedge[3*i+k0];
    5200                                                         MeshVertex * v02 = vertices + kedge[3*i+k1];
    5201                                                         MeshVertex * v01 = vertices + kedge[3*i+k2];
    5202                                                         t0(i1) = v01;
    5203                                                         t0(i2) = v02;
    5204                                                         t0.SetAllFlag(k0,hid[k0]);
    5205 
    5206                                                         t1(i0) = v01;
    5207                                                         t1(i2) = v12;
    5208                                                         t0.SetAllFlag(k1,hid[k1]);
    5209 
    5210                                                         t2(i0) = v02;
    5211                                                         t2(i1) = v12;
    5212                                                         t2.SetAllFlag(k2,hid[k2]);
    5213 
    5214                                                         t3(i0) = v12;
    5215                                                         t3(i1) = v02;
    5216                                                         t3(i2) = v01;
    5217 
    5218                                                         t3.SetAllFlag(0,hid[0]);           
    5219                                                         t3.SetAllFlag(1,hid[1]);           
    5220                                                         t3.SetAllFlag(2,hid[2]);
    5221 
    5222                                                         if ( kk == 6)
    5223                                                           {
    5224 
    5225                                                                 Triangle &t4=triangles[kkk++];
    5226                                                                 Triangle &t5=triangles[kkk++];
    5227 
    5228                                                                 t4 = t3;
    5229                                                                 t5 = t3;
    5230 
    5231                                                                 t0.SetHidden(k0);
    5232                                                                 t1.SetHidden(k1);
    5233                                                                 t2.SetHidden(k2);
    5234                                                                 t3.SetHidden(0);
    5235                                                                 t4.SetHidden(1);
    5236                                                                 t5.SetHidden(2);
    5237 
    5238                                                                 if (nbv < nbvx )
    5239                                                                   {
    5240                                                                         vertices[nbv].r = ((R2) *v01 + (R2) *v12  + (R2) *v02 ) / 3.0;
    5241                                                                         vertices[nbv].ReferenceNumber =0;
    5242                                                                         vertices[nbv].DirOfSearch =NoDirOfSearch;
    5243                                                                         //vertices[nbv].i = toI2(vertices[nbv].r);
    5244                                                                         double a3[]={1./3.,1./3.,1./3.};
    5245                                                                         vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m);
    5246                                                                         MeshVertex * vc =  vertices +nbv++;
    5247                                                                         t3(i0) = vc;
    5248                                                                         t4(i1) = vc;
    5249                                                                         t5(i2) = vc;
    5250 
    5251                                                                   }
    5252                                                                 else
    5253                                                                  goto Error;
    5254                                                           }
    5255 
    5256                                                   }
    5257                                                 break;         
    5258                         }
    5259 
    5260                         // save all the new triangles
    5261                         mkadj[nbmkadj++]=i;
    5262                         long jj;
    5263                         if (t0.link)
    5264                          for (jj=nbt;jj<kkk;jj++)
    5265                                 {
    5266                                  triangles[jj].link=t0.link;
    5267                                  t0.link= triangles+jj;
    5268                                  mkadj[nbmkadj++]=jj;
    5269                                 }
    5270                         if (nbmkadj>13){// 13 = 6 + 4 +
    5271                                 ISSMERROR("nbmkadj>13");
    5272                         }
    5273 
    5274                         if (kk==6)  newNbOfQuad+=3;
    5275                         for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;
    5276                         ksplit[i]= nbt; // save last adresse of the new triangles
    5277                         kkk = nbt;
    5278                   }
    5279 
    5280                 for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.;
    5281 
    5282                 if(withBackground)
    5283                  for (i=0;i<BTh.nbv;i++)
    5284                   BTh.vertices[i].m =  BTh.vertices[i].m*2.;
    5285 
    5286 
    5287                 ret = 2;
    5288                 if (nbt>= nbtx) goto Error; // bug
    5289                 if (nbv>= nbvx) goto Error; // bug
    5290                 // generation of the new triangles
    5291 
    5292                 SetIntCoor("In SplitElement");
    5293 
    5294                 ReMakeTriangleContainingTheVertex();
    5295                 if(withBackground)
    5296                  BTh.ReMakeTriangleContainingTheVertex();
    5297 
    5298                 delete [] edges;
    5299                 edges = newedges;
    5300                 nbe = newnbe;
    5301                 NbOfQuad = newNbOfQuad;
    5302 
    5303                 for (i=0;i<NbSubDomains;i++)
    5304                   {
    5305                         long k = subdomains[i].edge- edges;
    5306                         subdomains[i].edge =  edges+2*k; // spilt all edge in 2
    5307                   }
    5308 
    5309                 if (ksplitarray) delete [] ksplitarray;
    5310                 if (kedge) delete [] kedge;
    5311                 if (edge4) delete edge4;
    5312                 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
    5313                 VerticesOnGeomEdge= newVerticesOnGeomEdge;
    5314                 if(VertexOnBThEdge) delete []  VertexOnBThEdge;
    5315                 VertexOnBThEdge = newVertexOnBThEdge;
    5316                 NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
    5317                 NbVertexOnBThEdge=newNbVertexOnBThEdge;
    5318                 //  ReMakeTriangleContainingTheVertex();
    5319 
    5320                 ReconstructExistingMesh();
    5321 
    5322                 if (verbose>2){
    5323                         printf("   number of quadrilaterals    = %i\n",NbOfQuad);
    5324                         printf("   number of triangles         = %i\n",nbt-NbOutT- NbOfQuad*2);
    5325                         printf("   number of outside triangles = %i\n",NbOutT);
    5326                 }
    5327 
    5328                 return 0; //ok
    5329 
    5330 Error:
    5331                 nbv = nbvold;
    5332                 nbt = nbtold;
    5333                 NbOutT = NbOutTold;
    5334                 // cleaning memory ---
    5335                 delete newedges;
    5336                 if (ksplitarray) delete [] ksplitarray;
    5337                 if (kedge) delete [] kedge;
    5338                 if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge;
    5339                 if (edge4) delete edge4;
    5340                 if(newVertexOnBThEdge) delete []  newVertexOnBThEdge;
    5341 
    5342                 return ret; // ok
    5343         }
    5344         /*}}}1*/
    5345 /*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{1*/
    5346 long  Mesh::SplitInternalEdgeWithBorderVertices(){
    5347         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)*/
    5348 
    5349         long NbSplitEdge=0;
    5350         SetVertexFieldOn(); 
    5351         long it;
    5352         long nbvold=nbv;
    5353         long int verbose=2;
    5354         for (it=0;it<nbt;it++){
    5355                 Triangle &t=triangles[it];
    5356                 if (t.link)
    5357                  for (int j=0;j<3;j++)
    5358                   if(!t.Locked(j) && !t.Hidden(j)){
    5359                           Triangle &tt = *t.TriangleAdj(j);
    5360                           if ( &tt && tt.link && it < Number(tt))
    5361                                  { // an internal edge
    5362                                   MeshVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
    5363                                   MeshVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
    5364                                   if (v0.onGeometry && v1.onGeometry){
    5365                                           R2 P= ((R2) v0 + (R2) v1)*0.5;
    5366                                           if ( nbv<nbvx) {
    5367                                                   vertices[nbv].r = P;
    5368                                                   vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);
    5369                                                   vertices[nbv].ReferenceNumber=0;
    5370                                                   vertices[nbv].DirOfSearch = NoDirOfSearch ;
    5371                                           }
    5372                                           NbSplitEdge++;
    5373                                   }
    5374                                  }
    5375                   }
    5376         }
    5377         ReMakeTriangleContainingTheVertex();   
    5378         if (nbvold!=nbv){
    5379                 long  iv = nbvold;
    5380                 long NbSwap = 0;
    5381                 Icoor2 dete[3]; 
    5382                 for (int i=nbvold;i<nbv;i++) {// for all the new point
    5383                         MeshVertex & vi = vertices[i];
    5384                         vi.i = toI2(vi.r);
    5385                         vi.r = toR2(vi.i);
    5386 
    5387                         // a good new point
    5388                         vi.ReferenceNumber=0;
    5389                         vi.DirOfSearch =NoDirOfSearch;
    5390                         Triangle *tcvi = FindTriangleContaining(vi.i,dete);
    5391                         if (tcvi && !tcvi->link) {
    5392                                 printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
    5393                         }
    5394 
    5395                         quadtree->Add(vi);
    5396                         if (!tcvi || tcvi->det<0){// internal
    5397                                 ISSMERROR("!tcvi || tcvi->det < 0");
    5398                         }
    5399                         AddVertex(vi,tcvi,dete);
    5400                         NbSwap += vi.Optim(1);         
    5401                         iv++;
    5402                         //      }
    5403         }
    5404         if (verbose>3) {
    5405                 printf("   number of points: %i\n",iv);
    5406                 printf("   number of swap to  split internal edges with border vertices: %i\n",NbSwap);
    5407                 nbv = iv;
    5408         }
    5409 }
    5410 if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices  to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold));
    5411 if (verbose>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge);
    5412 
    5413 return  NbSplitEdge;
    5414 }
    5415 /*}}}1*/
    54165405/*FUNCTION Mesh::TriangleReferenceList{{{1*/
    54175406long  Mesh::TriangleReferenceList(long* reft) const {
    … …  
    54845473                                ISSMERROR("kkk>=1000");
    54855474                        }
    5486                         MeshVertex  &vI =  *edge.EdgeVertex(0);
    5487                         MeshVertex  &vJ =  *edge.EdgeVertex(1);
     5475                        BamgVertex  &vI =  *edge.EdgeVertex(0);
     5476                        BamgVertex  &vJ =  *edge.EdgeVertex(1);
    54885477                        I2 I=vI, J=vJ, IJ= J-I;
    54895478                        IJ_IA = (IJ ,(A-I));
    … …  
    55255514                // the probleme is in case of  the fine and long internal hole
    55265515                // for exemple neart the training edge of a wing
    5527                 MeshVertex * s=0,*s1=0, *s0=0;
     5516                BamgVertex * s=0,*s1=0, *s0=0;
    55285517                Icoor2 imax = MaxICoor22;
    55295518                Icoor2 l0 = imax,l1 = imax;
    … …  
    56125601                                if ((link + linkp) == 1)
    56135602                                  { // a boundary edge
    5614                                         MeshVertex * st = edge.EdgeVertex(1);
     5603                                        BamgVertex * st = edge.EdgeVertex(1);
    56155604                                        I2 I=*st;
    56165605                                        Icoor2  ll = Norme2_2 (C-I);
    … …  
    56575646        /*}}}1*/
    56585647/*FUNCTION ForceEdge{{{1*/
    5659 int ForceEdge(MeshVertex &a, MeshVertex & b,TriangleAdjacent & taret)  {
     5648int ForceEdge(BamgVertex &a, BamgVertex & b,TriangleAdjacent & taret)  {
    56605649        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/
    56615650
    … …  
    56685657
    56695658        TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]);
    5670         MeshVertex   *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
     5659        BamgVertex   *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
    56715660        // we turn around a in the  direct sens 
    56725661
    … …  
    56935682                if((det1 < 0) && (det2 >0)) {
    56945683                        // try to force the edge
    5695                         MeshVertex * va = &a, *vb = &b;
     5684                        BamgVertex * va = &a, *vb = &b;
    56965685                        tc = Previous(tc);
    56975686                        if (!v1 || !v2){
    … …  
    57035692                                 ISSMERROR("Loop in forcing Egde, nb de swap=%i, nb of try swap (%i) too big",NbSwap,l);
    57045693                         }
    5705                         MeshVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
     5694                        BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
    57065695                        if (( aa == &a ) && (bb == &b) ||  (bb ==  &a ) && (aa == &b)) {
    57075696                                tc.SetLock();
    … …  
    57335722/*}}}1*/
    57345723/*FUNCTION swap{{{1*/
    5735 void  swap(Triangle *t1,short a1, Triangle *t2,short a2, MeshVertex *s1,MeshVertex *s2,Icoor2 det1,Icoor2 det2){
     5724void  swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2){
    57365725        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)*/
    57375726        // --------------------------------------------------------------
    … …  
    57785767/*}}}1*/
    57795768        /*FUNCTION SwapForForcingEdge{{{1*/
    5780         int SwapForForcingEdge(MeshVertex   *  & pva ,MeshVertex  * &   pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) {
     5769        int SwapForForcingEdge(BamgVertex   *  & pva ,BamgVertex  * &   pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) {
    57815770                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/
    57825771                // l'arete ta coupe l'arete pva pvb
    … …  
    57955784                }
    57965785
    5797                 MeshVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
    5798                 MeshVertex & s1= (*t1)[OppositeVertex[a1]];
    5799                 MeshVertex & s2= (*t2)[OppositeVertex[a2]];
     5786                BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
     5787                BamgVertex & s1= (*t1)[OppositeVertex[a1]];
     5788                BamgVertex & s2= (*t2)[OppositeVertex[a2]];
    58005789
    58015790

  • issm/trunk/src/c/objects/Bamg/Mesh.h

    r5095 r5120  
    2222
    2323        class Mesh {
     24
    2425                public:
    2526
    26                         Geometry & Gh;   // Geometry
    27                         Mesh & BTh; // Background Mesh Bth==*this =>no  background
    28                         long NbRef;      // counter of ref on the this class if 0 we can delete
    29                         long nbvx,nbtx;  // nombre max  de sommets , de triangles
    30                         long nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference,
    31                         long NbOfQuad;  // nb of quadrangle
    32                         long NbSubDomains;
    33                         long NbOutT;    // Nb of oudeside triangle
    34                         long NbOfTriangleSearchFind;
    35                         long NbOfSwapTriangle;
    36                         MeshVertex* vertices;
    37                         long NbVerticesOnGeomVertex;
    38                         VertexOnGeom * VerticesOnGeomVertex;
    39                         long NbVerticesOnGeomEdge;
    40                         VertexOnGeom * VerticesOnGeomEdge;
    41                         long NbVertexOnBThVertex;
    42                         VertexOnVertex *VertexOnBThVertex;
    43                         long NbVertexOnBThEdge;
    44                         VertexOnEdge *VertexOnBThEdge;
    45                         long NbCrackedVertices;
    46                         long* CrackedVertices;
    47                         long NbCrackedEdges;
    48                         CrackedEdge* CrackedEdges;
    49                         R2 pmin,pmax;    // extrema
    50                         double coefIcoor; // coef to integer Icoor1;
    51                         Triangle* triangles;
    52                         Edge* edges;
    53                         QuadTree* quadtree;
    54                         MeshVertex** ordre;
    55                         SubDomain* subdomains;
    56                         ListofIntersectionTriangles lIntTria;
     27                        Geometry                    & Gh;                    // Geometry
     28                        Mesh                        & BTh;                   // Background Mesh Bth== *this =>no background
     29                        long                          NbRef;                 // counter of ref on the this class if 0 we can delete
     30                        long                          maxnbv,nbtx;           // nombre max de sommets , de triangles
     31                        long                          nt,nbv,nbt,nbe;        // nb of legal triangles, nb of vertex, of triangles, of edges with reference,
     32                        long                          NbOfQuad;              // nb of quadrangle
     33                        long                          NbSubDomains;
     34                        long                          NbOutT;                // Nb of oudeside triangle
     35                        long                          NbOfTriangleSearchFind;
     36                        long                          NbOfSwapTriangle;
     37                        BamgVertex                   *vertices;
     38                        long                          NbVerticesOnGeomVertex;
     39                        VertexOnGeom                 *VerticesOnGeomVertex;
     40                        long                          NbVerticesOnGeomEdge;
     41                        VertexOnGeom                 *VerticesOnGeomEdge;
     42                        long                          NbVertexOnBThVertex;
     43                        VertexOnVertex               *VertexOnBThVertex;
     44                        long                          NbVertexOnBThEdge;
     45                        VertexOnEdge                 *VertexOnBThEdge;
     46                        long                          NbCrackedVertices;
     47                        long                         *CrackedVertices;
     48                        long                          NbCrackedEdges;
     49                        CrackedEdge                  *CrackedEdges;
     50                        R2                            pmin,pmax;             // extrema
     51                        double                        coefIcoor;            // coef to integer Icoor1;
     52                        Triangle                     *triangles;
     53                        Edge                         *edges;
     54                        QuadTree                     *quadtree;
     55                        BamgVertex                  **ordre;
     56                        SubDomain                    *subdomains;
     57                        ListofIntersectionTriangles   lIntTria;
    5758
    5859                        //Constructors/Destructors
    5960                        Mesh(BamgGeom* bamggeom,BamgMesh* bamgmesh,BamgOpts* bamgopts);
    6061                        Mesh(double* index,double* x,double* y,int nods,int nels);
    61                         Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long nbvxx=0 ); //copy operator
     62                        Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long maxnbv_in=0 ); //copy operator
    6263                        Mesh(const Mesh &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
    63                         Mesh(long nbvx,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
    64                                 try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);}
     64                        Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
     65                                try {TriangulateFromGeom1(maxnbv,bamgopts,keepBackVertices);}
    6566                                catch(...) { this->~Mesh(); throw; }
    6667                        }
    67                         Mesh(long nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
    68                                 try { GeomToTriangles0(nbvx,bamgopts);}
     68                        Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
     69                                try { TriangulateFromGeom0(maxnbv,bamgopts);}
    6970                                catch(...) { this->~Mesh(); throw; }
    7071                        }
    … …  
    7273
    7374                        //Operators
    74                         const MeshVertex & operator[]  (long i) const { return vertices[i];};
    75                         MeshVertex & operator[](long i) { return vertices[i];};
     75                        const BamgVertex & operator[]  (long i) const { return vertices[i];};
     76                        BamgVertex & operator[](long i) { return vertices[i];};
    7677                        const Triangle & operator()  (long i) const { return triangles[i];};
    7778                        Triangle & operator()(long i) { return triangles[i];};
    … …  
    8788                                return  R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
    8889                        }
    89                         void AddVertex(MeshVertex & s,Triangle * t,Icoor2 *  =0) ;
     90                        void AddVertex(BamgVertex & s,Triangle * t,Icoor2 *  =0) ;
    9091                        void Insert();
    9192                        void ForceBoundary();
    … …  
    110111                        void SmoothingVertex(int =3,double=0.3);
    111112                        Metric MetricAt (const R2 &) const;
    112                         GeometricalEdge* ProjectOnCurve( Edge & AB, MeshVertex &  A, MeshVertex & B,double theta, MeshVertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
     113                        GeometricalEdge* ProjectOnCurve( Edge & AB, BamgVertex &  A, BamgVertex & B,double theta, BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
    113114                        long Number(const Triangle & t) const  { return &t - triangles;}
    114115                        long Number(const Triangle * t) const  { return t - triangles;}
    115                         long Number(const MeshVertex & t) const  { return &t - vertices;}
    116                         long Number(const MeshVertex * t) const  { return t - vertices;}
     116                        long Number(const BamgVertex & t) const  { return &t - vertices;}
     117                        long Number(const BamgVertex * t) const  { return t - vertices;}
    117118                        long Number(const Edge & t) const  { return &t - edges;}
    118119                        long Number(const Edge * t) const  { return t - edges;}
    119120                        long Number2(const Triangle * t) const  { return t - triangles; }
    120                         MeshVertex* NearestVertex(Icoor1 i,Icoor1 j) ;
     121                        BamgVertex* NearestVertex(Icoor1 i,Icoor1 j) ;
    121122                        Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
    122123                        void ReadMesh(double* index,double* x,double* y,int nods,int nels);
    … …  
    155156
    156157                private:
    157                         void GeomToTriangles1(long nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption
    158                         void GeomToTriangles0(long nbvx,BamgOpts* bamgopts);// the real constructor mesh generator
    159                         void PreInit(long);
     158                        void TriangulateFromGeom1(long maxnbv,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption
     159                        void TriangulateFromGeom0(long maxnbv,BamgOpts* bamgopts);// the real constructor mesh generator
     160                        void Init(long);
    160161
    161162        };
    … …  
    166167        void  swap(Triangle *t1,short a1,
    167168                                Triangle *t2,short a2,
    168                                 MeshVertex *s1,MeshVertex *s2,Icoor2 det1,Icoor2 det2);
    169         int SwapForForcingEdge(MeshVertex   *  & pva ,MeshVertex  * &   pvb ,
     169                                BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2);
     170        int SwapForForcingEdge(BamgVertex   *  & pva ,BamgVertex  * &   pvb ,
    170171                                TriangleAdjacent & tt1,Icoor2 & dets1,
    171172                                Icoor2 & detsa,Icoor2 & detsb, int & nbswap);
    172         int ForceEdge(MeshVertex &a, MeshVertex & b,TriangleAdjacent & taret) ;
     173        int ForceEdge(BamgVertex &a, BamgVertex & b,TriangleAdjacent & taret) ;
    173174        inline TriangleAdjacent Previous(const TriangleAdjacent & ta){
    174175                return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);
    … …  
    183184                int j=i;i=on->DirAdj[i];on=on->Adj[j];
    184185        }
    185         inline double qualite(const MeshVertex &va,const MeshVertex &vb,const MeshVertex &vc){
     186        inline double qualite(const BamgVertex &va,const BamgVertex &vb,const BamgVertex &vc){
    186187                double ret;
    187188                I2 ia=va,ib=vb,ic=vc;

  • issm/trunk/src/c/objects/Bamg/QuadTree.cpp

    r5095 r5120  
    2424        /*}}}*/
    2525        /*DOCUMENTATION What is a QuadTree? {{{1
    26          * A Quadtree is a very simple way to group the vertices according
    27          * to their location. A square that holds all the points of the mesh
    28          * (or the geometry) is divided into 4 boxes. As soom as one box
     26         * A Quadtree is a very simple way to group vertices according
     27         * to their locations. A square that holds all the points of the mesh
     28         * (or the geometry) is divided into 4 boxes. As soon as one box
    2929         * hold more than 4 vertices, it is divided into 4 new boxes, etc...
    3030         * There cannot be more than MAXDEEP (=30) subdivision.
    … …  
    6767         * 0 1 1 1 .... 1    0 1 1 1 .... 1 in binary
    6868         *  \--   29  --/     \--   29  --/
    69          * Using the binaries is therefor very easy to locate a vertex in a box:
     69         * Using binaries is therefore very easy to locate a vertex in a box:
    7070         * we just need to look at the bits from the left to the right (See ::Add)
    7171         }}}1*/
    … …  
    7676                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/QuadTree)*/
    7777
    78                 lenStorageQuadTreeBox(t->nbvx/8+10),
     78                lenStorageQuadTreeBox(t->maxnbv/8+10),
    7979                th(t),
    8080                NbQuadTreeBox(0),
    … …  
    120120        /*Methods*/
    121121        /*FUNCTION QuadTree::Add{{{1*/
    122         void  QuadTree::Add(MeshVertex &w){
     122        void  QuadTree::Add(BamgVertex &w){
    123123                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/Add)*/
    124124
    … …  
    163163
    164164                        //Copy the 4 vertices in the current QuadTreebox
    165                         MeshVertex* v4[4];
     165                        BamgVertex* v4[4];
    166166                        v4[0]= b->v[0];
    167167                        v4[1]= b->v[1];
    … …  
    205205        /*}}}1*/
    206206        /*FUNCTION QuadTree::NearestVertex{{{1*/
    207         MeshVertex*  QuadTree::NearestVertex(Icoor1 i,Icoor1 j) {
     207        BamgVertex*  QuadTree::NearestVertex(Icoor1 i,Icoor1 j) {
    208208                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/NearestVertex)*/
    209209
    … …  
    224224                Icoor1   jplus( j<MaxISize?(j<0?0:j):MaxISize-1);
    225225                //initial nearest vertex pointer
    226                 MeshVertex*  vn=NULL;
     226                BamgVertex*  vn=NULL;
    227227
    228228                //Get initial Quadtree box (largest)
    … …  
    293293
    294294                                //if the current subbox is holding vertices,
    295                                 if (b->n>0){ // MeshVertex QuadTreeBox not empty
     295                                if (b->n>0){ // BamgVertex QuadTreeBox not empty
    296296                                        NbVerticesSearch++;
    297297                                        I2 i2 =  b->v[k]->i;
    … …  
    344344        /*}}}1*/
    345345        /*FUNCTION QuadTree::NearestVertexWithNormal{{{1*/
    346         MeshVertex*  QuadTree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) {
     346        BamgVertex*  QuadTree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) {
    347347                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/NearestVertexWithNormal)*/
    348348
    … …  
    358358                Icoor1  jplus( j<MaxISize?(j<0?0:j):MaxISize-1);
    359359
    360                 MeshVertex *vn=0;
     360                BamgVertex *vn=0;
    361361
    362362                // init for optimisation ---
    … …  
    411411                                int k = pi[l];
    412412
    413                                 if (b->n>0) // MeshVertex QuadTreeBox none empty
     413                                if (b->n>0) // BamgVertex QuadTreeBox none empty
    414414                                  {
    415415                                        NbVerticesSearch++;
    … …  
    472472        /*}}}1*/
    473473        /*FUNCTION QuadTree::ToClose {{{1*/
    474         MeshVertex *   QuadTree::ToClose(MeshVertex & v,double seuil,Icoor1 hx,Icoor1 hy){
     474        BamgVertex *   QuadTree::ToClose(BamgVertex & v,double seuil,Icoor1 hx,Icoor1 hy){
    475475                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/ToClose)*/
    476476
    … …  
    489489                Icoor1 i0=0,j0=0;
    490490
    491                 //  MeshVertex *vn=0;
     491                //  BamgVertex *vn=0;
    492492
    493493                if (!root->n)
    … …  
    506506                                register int k = pi[l];
    507507
    508                                 if (b->n>0) // MeshVertex QuadTreeBox none empty
     508                                if (b->n>0) // BamgVertex QuadTreeBox none empty
    509509                                  {
    510510                                        NbVerticesSearch++;

  • issm/trunk/src/c/objects/Bamg/QuadTree.h

    r5095 r5120  
    1111
    1212        class Mesh;
    13         class MeshVertex;
     13        class BamgVertex;
    1414
    1515        class QuadTree{
    … …  
    1818                                public:
    1919                                        long n;
    20                                         //contains only one object form the list (either MeshVertex or QuadTreeBox)
    21                                         // if n < 4 => MeshVertex else =>  QuadTreeBox;
     20                                        //contains only one object form the list (either BamgVertex or QuadTreeBox)
     21                                        // if n < 4 => BamgVertex else =>  QuadTreeBox;
    2222                                        union{
    2323                                                QuadTreeBox* b[4];
    24                                                 MeshVertex* v[4];
     24                                                BamgVertex* v[4];
    2525                                        };
    2626                        };
    … …  
    5050                        long    NbQuadTreeBox,NbVertices;
    5151                        long    NbQuadTreeBoxSearch,NbVerticesSearch;
    52                         MeshVertex* NearestVertex(Icoor1 i,Icoor1 j);
    53                         MeshVertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j);
    54                         MeshVertex* ToClose(MeshVertex & ,double ,Icoor1,Icoor1);
     52                        BamgVertex* NearestVertex(Icoor1 i,Icoor1 j);
     53                        BamgVertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j);
     54                        BamgVertex* ToClose(BamgVertex & ,double ,Icoor1,Icoor1);
    5555                        long    SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();}
    56                         void    Add( MeshVertex & w);
     56                        void    Add( BamgVertex & w);
    5757                        QuadTreeBox* NewQuadTreeBox(){
    5858                                if(! (sb->bc<sb->be)) sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);

  • issm/trunk/src/c/objects/Bamg/R2.h

    r3913 r5120  
    77namespace bamg {
    88
    9         template <class R,class RR> class P2xP2;
    10         template <class R,class RR>
    11           class P2 {
     9        template <class R,class RR> class P2 {
     10
    1211                  public: 
     12
    1313                          //fields
    1414                          R x,y;
     15
    1516                          //functions
    1617                          P2 () :x(0),y(0) {};
    … …  
    1920                          void Echo(){
    2021                                  printf("Member of P2:\n");
    21                                   printf("   x: %g\n",x);
    22                                   printf("   y: %g\n",y);
     22                                  printf("   x: %g or %i\n",x,x);
     23                                  printf("   y: %g or %i\n",y,y);
    2324                          }
    2425                          //operators
    … …  
    3334                          P2<R,RR> operator*=(const  R r) {x *= r;y *= r;return *this;}
    3435                          P2<R,RR> operator-=(const  P2<R,RR> & cc) {x -= cc.x;y -= cc.y;return *this;}
     36
    3537          };
    3638
    37         template <class R,class RR>
    38           class P2xP2 {
     39        template <class R,class RR> class P2xP2 {
     40
    3941                  private:
     42
    4043                          friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc){
    4144                                  return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y);
    4245                          }
    4346                  public:
     47
    4448                          //fields
    4549                          P2<R,RR> x,y;
     50
    4651                          //functions
    4752                          P2xP2 (): x(),y()  {}

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

    r5095 r5120  
    1111        /*FUNCTION Triangle(Mesh *Th,long i,long j,long k) {{{1*/
    1212        Triangle::Triangle(Mesh *Th,long i,long j,long k) {
    13                 MeshVertex *v=Th->vertices;
     13                BamgVertex *v=Th->vertices;
    1414                long nbv = Th->nbv;
    1515                if (i<0 || j<0 || k<0){
    … …  
    2727        }
    2828        /*}}}*/
    29         /*FUNCTION Triangle(MeshVertex *v0,MeshVertex *v1,MeshVertex *v2) {{{1*/
    30         Triangle::Triangle(MeshVertex *v0,MeshVertex *v1,MeshVertex *v2){
     29        /*FUNCTION Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2) {{{1*/
     30        Triangle::Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2){
    3131                TriaVertices[0]=v0;
    3232                TriaVertices[1]=v1;
    … …  
    154154        /*}}}1*/
    155155        /*FUNCTION Triangle::Quadrangle {{{1*/
    156         Triangle* Triangle::Quadrangle(MeshVertex * & v0,MeshVertex * & v1,MeshVertex * & v2,MeshVertex * & v3) const{
     156        Triangle* Triangle::Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const{
    157157                // return the other triangle of the quad if a quad or 0 if not a quat
    158158                Triangle * t =0;
    … …  
    186186                         q= -1;
    187187                        else if(option){
    188                                 const MeshVertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
    189                                 const MeshVertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
    190                                 const MeshVertex & v1 = *TriaVertices[OppositeEdge[a]];
    191                                 const MeshVertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
     188                                const BamgVertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
     189                                const BamgVertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
     190                                const BamgVertex & v1 = *TriaVertices[OppositeEdge[a]];
     191                                const BamgVertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
    192192                                q =  QuadQuality(v0,v1,v2,v3); // do the float part
    193193                        }
    … …  
    220220                if(a2/4 !=0) return 0; // arete lock or MarkUnSwap
    221221
    222                 register MeshVertex  *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]];
    223                 register MeshVertex  *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]];
    224                 register MeshVertex  *s1=t1->TriaVertices[OppositeVertex[a1]];
    225                 register MeshVertex  *s2=t2->TriaVertices[OppositeVertex[a2]];
     222                register BamgVertex  *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]];
     223                register BamgVertex  *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]];
     224                register BamgVertex  *s1=t1->TriaVertices[OppositeVertex[a1]];
     225                register BamgVertex  *s2=t2->TriaVertices[OppositeVertex[a2]];
    226226
    227227                Icoor2 det1=t1->det , det2=t2->det ;

  • issm/trunk/src/c/objects/Bamg/Triangle.h

    r5095 r5120  
    99        //classes
    1010        class Mesh;
    11         class MeshVertex;
     11        class BamgVertex;
    1212        class Triangle;
    1313
    … …  
    1717
    1818                private:
    19                         MeshVertex*   TriaVertices[3];      // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
     19                        BamgVertex*   TriaVertices[3];      // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
    2020                        Triangle* TriaAdjTriangles[3];  // 3 pointers toward the adjacent triangles
    2121                        short     TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1
    … …  
    3131                        Triangle() {}
    3232                        Triangle(Mesh *Th,long i,long j,long k);
    33                         Triangle(MeshVertex *v0,MeshVertex *v1,MeshVertex *v2);
     33                        Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2);
    3434
    3535                        //Operators
    36                         const MeshVertex & operator[](int i) const {return *TriaVertices[i];};
    37                         MeshVertex & operator[](int i)  {return *TriaVertices[i];};
    38                         const MeshVertex * operator()(int i) const {return TriaVertices[i];};
    39                         MeshVertex * & operator()(int i)  {return TriaVertices[i];};
     36                        const BamgVertex & operator[](int i) const {return *TriaVertices[i];};
     37                        BamgVertex & operator[](int i)  {return *TriaVertices[i];};
     38                        const BamgVertex * operator()(int i) const {return TriaVertices[i];};
     39                        BamgVertex * & operator()(int i)  {return TriaVertices[i];};
    4040
    4141                        //Methods
    … …  
    5252                        TriangleAdjacent Adj(int i)  const {return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);};
    5353                        Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}
    54                         Triangle* Quadrangle(MeshVertex * & v0,MeshVertex * & v1,MeshVertex * & v2,MeshVertex * & v3) const ;
     54                        Triangle* Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const ;
    5555                        void  ReNumbering(Triangle *tb,Triangle *te, long *renu){
    5656                                if (link  >=tb && link  <te) link  = tb + renu[link -tb];
    … …  
    5959                                if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];   
    6060                        }
    61                         void ReNumbering(MeshVertex *vb,MeshVertex *ve, long *renu){
     61                        void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
    6262                                if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
    6363                                if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
    … …  
    118118                        double qualite() ;
    119119                        void  Set(const Triangle &,const Mesh &,Mesh &);
    120                         int   In(MeshVertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
     120                        int   In(BamgVertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
    121121
    122122        };

  • issm/trunk/src/c/objects/Bamg/TriangleAdjacent.cpp

    r5095 r5120  
    4141        /*}}}*/
    4242        /*FUNCTION TriangleAdjacent::EdgeVertex {{{1*/
    43         MeshVertex* TriangleAdjacent::EdgeVertex(const int & i) const {
     43        BamgVertex* TriangleAdjacent::EdgeVertex(const int & i) const {
    4444                return t->TriaVertices[VerticesOfTriangularEdge[a][i]];
    4545        }
    4646        /*}}}*/
    4747        /*FUNCTION TriangleAdjacent::OppositeVertex {{{1*/
    48         MeshVertex* TriangleAdjacent::OppositeVertex() const {
     48        BamgVertex* TriangleAdjacent::OppositeVertex() const {
    4949                return t->TriaVertices[bamg::OppositeVertex[a]];
    5050        }

  • issm/trunk/src/c/objects/Bamg/TriangleAdjacent.h

    r3913 r5120  
    33
    44#include "./include.h"
    5 #include "./MeshVertex.h"
     5#include "./BamgVertex.h"
    66
    77namespace bamg {
    … …  
    3838                        int  swap();
    3939                        TriangleAdjacent Adj() const;
    40                         MeshVertex* EdgeVertex(const int & i) const;
    41                         MeshVertex* OppositeVertex() const;
     40                        BamgVertex* EdgeVertex(const int & i) const;
     41                        BamgVertex* OppositeVertex() const;
    4242                        Icoor2& det() const;
    4343        };

  • issm/trunk/src/c/objects/Bamg/VertexOnEdge.h

    r5095 r5120  
    99        //classes
    1010        class Mesh;
    11         class MeshVertex;
     11        class BamgVertex;
    1212
    1313        class VertexOnEdge {
    1414
    1515                public:
    16                         MeshVertex* v;
     16                        BamgVertex* v;
    1717                        Edge*   be;
    1818                        double abcisse;
    1919
    2020                        //Constructors
    21                         VertexOnEdge( MeshVertex * w, Edge *bw,double s) :v(w),be(bw),abcisse(s) {}
     21                        VertexOnEdge( BamgVertex * w, Edge *bw,double s) :v(w),be(bw),abcisse(s) {}
    2222                        VertexOnEdge(){}
    2323
    2424                        //Operators
    2525                        operator double () const { return abcisse;}
    26                         operator MeshVertex* () const { return v;} 
     26                        operator BamgVertex* () const { return v;} 
    2727                        operator Edge* () const { return be;} 
    28                         MeshVertex & operator[](int i) const { return (*be)[i];}
     28                        BamgVertex & operator[](int i) const { return (*be)[i];}
    2929
    3030                        //Methods

  • issm/trunk/src/c/objects/Bamg/VertexOnGeom.h

    r5095 r5120  
    99        //classes
    1010        class Mesh;
    11         class MeshVertex;
     11        class BamgVertex;
    1212        class GeometricalEdge;
    1313
    … …  
    1616                public:
    1717
    18                         MeshVertex* mv;
     18                        BamgVertex* mv;
    1919                        double abscisse; 
    2020                        union{
    … …  
    2525                        //Constructors/Destructors
    2626                        VertexOnGeom(): mv(0),abscisse(0){gv=0;}
    27                         VertexOnGeom(MeshVertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}
    28                         VertexOnGeom(MeshVertex & m,GeometricalEdge &g,double s) : mv(&m),abscisse(s){ge=&g;}
     27                        VertexOnGeom(BamgVertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}
     28                        VertexOnGeom(BamgVertex & m,GeometricalEdge &g,double s) : mv(&m),abscisse(s){ge=&g;}
    2929
    3030                        //Operators
    31                         operator MeshVertex*() const  {return mv;}
     31                        operator BamgVertex*() const  {return mv;}
    3232                        operator GeometricalVertex * () const  {return gv;}
    3333                        operator GeometricalEdge * () const  {return ge;}

  • issm/trunk/src/c/objects/Bamg/VertexOnVertex.h

    r5095 r5120  
    33
    44#include "./include.h"
    5 #include "./MeshVertex.h"
     5#include "./BamgVertex.h"
    66
    77namespace bamg {
    … …  
    1313
    1414                public:
    15                         MeshVertex* v;
    16                         MeshVertex* bv;
     15                        BamgVertex* v;
     16                        BamgVertex* bv;
    1717
    1818                        //Constructors
    19                         VertexOnVertex(MeshVertex * w,MeshVertex *bw) :v(w),bv(bw){}
     19                        VertexOnVertex(BamgVertex * w,BamgVertex *bw) :v(w),bv(bw){}
    2020                        VertexOnVertex() {};
    2121

  • issm/trunk/src/c/objects/Bamg/typedefs.h

    r3913 r5120  
    66namespace bamg {
    77
    8         typedef int  Icoor1; 
    9 #if LONG_BIT > 63 //64 bits or more
     8        /*Integer coordinates types*/
     9        typedef int  Icoor1;
     10        #if LONG_BIT > 63 //64 bits or more
    1011        typedef long Icoor2;
    11 #else //32 bits
     12        #else //32 bits
    1213        typedef double Icoor2;
    13 #endif
     14        #endif
     15
     16        /*I2 and R2*/
    1417        typedef P2<Icoor1,Icoor2>  I2;
    1518        typedef P2xP2<short,long>  I2xI2;

  • issm/trunk/src/c/objects/objects.h

    r5095 r5120  
    9999#include "./Bamg/DoubleAndInt.h"
    100100#include "./Bamg/Direction.h"
    101 #include "./Bamg/MeshVertex.h"
     101#include "./Bamg/BamgVertex.h"
    102102#include "./Bamg/TriangleAdjacent.h"
    103103#include "./Bamg/Edge.h"
Note: See TracChangeset for help on using the changeset viewer.